This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / d / cft.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1996/04/01 15:02:24  mclareni
6 * Mathlib gen
7 *
8 *
9 #include "gen/pilot.h"
10       SUBROUTINE CFT(A,B,NTOT,N,NSPAN,ISN)
11 C
12 C     MULTIVARIATE COMPLEX FOURIER TRANSFORM, COMPUTED IN PLACE
13 C     USING MIXED-RADIX FAST FOURIER TRANSFORM ALGORITHM.
14 C     BY R. C. SINGLETON, STANFORD RESEARCH INSTITUTE, OCT. 1968
15 C     ARRAYS A AND B ORIGINALLY HOLD THE REAL AND IMAGINARY
16 C     COMPONENTS OF THE DATA, AND RETURN THE REAL AND
17 C     IMAGINARY COMPONENTS OF THE RESULTING FOURIER COEFFICIENTS.
18 C     MULTIVARIATE DATA IS INDEXED ACCORDING TO THE FORTRAN
19 C     ARRAY ELEMENT SUCCESSOR FUNCTION, WITHOUT LIMIT
20 C     ON THE NUMBER OF IMPLIED MULTIPLE SUBSCRIPTS.
21 C     THE SUBROUTINE IS CALLED ONCE FOR EACH VARIATE.
22 C     THE CALLS FOR A MULTIVARIATE TRANSFORM MAY BE IN ANY ORDER.
23 C     NTOT IS THE TOTAL NUMBER OF COMPLEX DATA VALUES.
24 C     N IS THE DIMENSION OF THE CURRENT VARIABLE.
25 C     NSPAN/N IS THE SPACING OF CONSUCUTIVE DATA VALUES
26 C     WHILE INDEXING THE CURRENT VARIABLE.
27 C     THE SIGN OF ISN DETERMINES THE SIGN OF THE COMPLEX
28 C     EXPONENTIAL, AND THE MAGNITUDE OF ISN IS NORMALLY ONE.
29 C
30 C     FOR A SINGLE-VARIATE TRANSFORM,
31 C     NTOT = N = NSPAN = (NUMBER OF COMPLEX DATA VALUES), F.G.
32 C     CALL CFT(A,B,N,N,N,1)
33 C
34 C     A TRI-VARIATE TRANSFORM WITH A(N1,N2,N3), B(N1,N2,N3)
35 C     IS COMPUTED BY
36 C     CALL CFT(A,B,N1*N2*N3,N1,N1,1)
37 C     CALL CFT(A,B,N1*N2*N3,N2,N1*N2,1)
38 C     CALL CFT(A,B,N1*N2*N3,N3,N1*N2*N3,1)
39 C
40 C     THE DATA MAY ALTERNATIVELY BE STORED IN A SINGLE COMPLEX
41 C     ARRAY A, THEN THE MAGNITUDE OF ISN CHANGED TO TWO TO
42 C     GIVE THE CORRECT INDEXING INCREMENT AND THE SECOND PARAMETER
43 C     USED TO PASS THE INITIAL ADDRESS FOR THE SEQUENCE OF
44 C     IMAGINARY VALUES, E.G.
45 C
46 C        REAL S(2)
47 C        EQUIVALENCE (A,S)
48 C        ....
49 C        ....
50 C        CALL CFT(A,S(2),NTOT,N,NSPAN,2)
51 C
52 C     ARRAYS AT(MAXF), CK(MAXF), BT(MAXF), SK(MAXF), AND NP(MAXP)
53 C     ARE USED FOR TEMPORARY STORAGE. IF THE AVAILABLE STORAGE
54 C     IS INSUFFICIENT, THE PROGRAM IS TERMINATED BY A STOP.
55 C     MAXF MUST BE .GE. THE MAXIMUM PRIME FACTOR OF N.
56 C     MAXP MUST BE .GT. THE NUMBER OF PRIME FACTORS OF N.
57 C     IN ADDITION, IF THE SQUARE-FREE PORTION K CF N HAS TWO OR
58 C     MORE PRIME FACTORS, THEN MAXP MUST BE .GE. K-1.
59 C     ARRAY STORAGE IN NFAC FOR A MAXIMUM OF 11 FACTORS OF N.
60 C     IF N HAS MORE THAN ONE SQUARE-FREE FACTOR, THE PRODUCT OF THE
61 C     SQUARE-FREE FACTORS MUST BE .LE. 210
62 C
63       DIMENSION A(1),B(1)
64       DIMENSION NFAC(11),NP(209)
65 C     ARRAY STORAGE FOR MAXIMUM PRIME FACTOR OF 23
66       DIMENSION AT(23),CK(23),BT(23),SK(23)
67       EQUIVALENCE (I,II)
68 C     THE FOLLOWING TWO CONSTANTS SHOULD AGREE WITH THE ARRAY DIMENSIONS
69       MAXF=23
70       MAXP=209
71       IF(N .LT. 2) RETURN
72       INC=ISN
73 C     THE FOLLOWING CONSTANTS ARE RAD = 2.*PI , S72 = SIN(0.4*PI) ,
74 C     C72 = COS(0.4*PI) AND S120 = SQRT(0.75)
75       RAD = 6.2831853071796
76       S72 = 0.95105651629515
77       C72 = 0.30901699437495
78       S120 = 0.86602540378444
79       IF(ISN .GE. 0) GO TO 10
80       S72=-S72
81       S120=-S120
82       RAD=-RAD
83       INC=-INC
84    10 NT=INC*NTOT
85       KS=INC*NSPAN
86       KSPAN=KS
87       NN=NT-INC
88       JC=KS/N
89       RADF=RAD*JC*0.5
90       I=0
91       JF=0
92 C     DETERMINE THE FACTORS OF N
93       M=0
94       K=N
95       GO TO 20
96    15 M=M+1
97       NFAC(M)=4
98       K=K/16
99    20 IF(K-(K/16)*16 .EQ. 0) GO TO 15
100       J=3
101       JJ=9
102       GO TO 30
103    25 M=M+1
104       NFAC(M)=J
105       K=K/JJ
106    30 IF(MOD(K,JJ) .EQ. 0) GO TO 25
107       J=J+2
108       JJ=J**2
109       IF(JJ .LE. K) GO TO 30
110       IF(K .GT. 4) GO TO 40
111       KT=M
112       NFAC(M+1)=K
113       IF(K .NE. 1) M=M+1
114       GO TO 80
115    40 IF(K-(K/4)*4 .NE. 0) GO TO 50
116       M=M+1
117       NFAC(M)=2
118       K=K/4
119    50 KT=M
120       J=2
121    60 IF(MOD(K,J) .NE. 0) GO TO 70
122       M=M+1
123       NFAC(M)=J
124       K=K/J
125    70 J=((J+1)/2)*2+1
126       IF(J .LE. K) GO TO 60
127    80 IF(KT .EQ. 0) GO TO 100
128       J=KT
129    90 M=M+1
130       NFAC(M)=NFAC(J)
131       J=J-1
132       IF(J .NE. 0) GO TO 90
133 C     COMPUTE FOURIER TRANSFORM
134   100 SD=RADF/KSPAN
135       CD=2.0*SIN(SD)**2
136       SD=SIN(SD+SD)
137       KK=1
138       I=I+1
139       IF(NFAC(I) .NE. 2) GO TO 400
140 C     TRANSFORM FOR FACTOR OF 2 (INCLUDING ROTATION FACTOR)
141       KSPAN=KSPAN/2
142       K1=KSPAN+2
143   210 K2=KK+KSPAN
144       AK=A(K2)
145       BK=B(K2)
146       A(K2)=A(KK)-AK
147       B(K2)=B(KK)-BK
148       A(KK)=A(KK)+AK
149       B(KK)=B(KK)+BK
150       KK=K2+KSPAN
151       IF(KK .LE. NN) GO TO 210
152       KK=KK-NN
153       IF(KK .LE. JC) GO TO 210
154       IF(KK .GT. KSPAN) GO TO 800
155   220 C1=1.0-CD
156       S1=SD
157   230 K2=KK+KSPAN
158       AK=A(KK)-A(K2)
159       BK=B(KK)-B(K2)
160       A(KK)=A(KK)+A(K2)
161       B(KK)=B(KK)+B(K2)
162       A(K2)=C1*AK-S1*BK
163       B(K2)=S1*AK+C1*BK
164       KK=K2+KSPAN
165       IF(KK .LT. NT) GO TO 230
166       K2=KK-NT
167       C1=-C1
168       KK=K1-K2
169       IF(KK .GT. K2) GO TO 230
170       AK=C1-(CD*C1+SD*S1)
171       S1=(SD*C1-CD*S1)+S1
172 C     THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
173 C     ERROR. IF ROUNDED ARITHMETIC IS USED, THEY MAY BE DELETED.
174 C     C1=0.5/(AK**2+S1**2)+0.5
175 C     S1=C1*S1
176 C     C1=C1*AK
177 C     NEXT STATEMENT SHOULD BE DELETED IF NON-ROUNDED ARITHMETIC IS USED
178       C1=AK
179       KK=KK+JC
180       IF(KK .LT. K2) GO TO 230
181       K1=K1+INC+INC
182       KK=(K1-KSPAN)/2+JC
183       IF(KK .LE. JC+JC) GO TO 220
184       GO TO 100
185 C     TRANSFORM FOR FACTOR OF 3 (OPTIONAL CODE)
186   320 K1=KK+KSPAN
187       K2=K1+KSPAN
188       AK=A(KK)
189       BK=B(KK)
190       AJ=A(K1)+A(K2)
191       BJ=B(K1)+B(K2)
192       A(KK)=AK+AJ
193       B(KK)=BK+BJ
194       AK=-0.5*AJ+AK
195       BK=-0.5*BJ+BK
196       AJ=(A(K1)-A(K2))*S120
197       BJ=(B(K1)-B(K2))*S120
198       A(K1)=AK-BJ
199       B(K1)=BK+AJ
200       A(K2)=AK+BJ
201       B(K2)=BK-AJ
202       KK=K2+KSPAN
203       IF(KK .LT. NN) GO TO 320
204       KK=KK-NN
205       IF(KK .LE. KSPAN) GO TO 320
206       GO TO 700
207 C     TRANSFORM FOR FACTOR OF 4
208   400 IF(NFAC(I) .NE. 4) GO TO 600
209       KSPNN=KSPAN
210       KSPAN=KSPAN/4
211   410 C1=1.0
212       S1=0
213   420 K1=KK+KSPAN
214       K2=K1+KSPAN
215       K3=K2+KSPAN
216       AKP=A(KK)+A(K2)
217       AKM=A(KK)-A(K2)
218       AJP=A(K1)+A(K3)
219       AJM=A(K1)-A(K3)
220       A(KK)=AKP+AJP
221       AJP=AKP-AJP
222       BKP=B(KK)+B(K2)
223       BKM=B(KK)-B(K2)
224       BJP=B(K1)+B(K3)
225       BJM=B(K1)-B(K3)
226       B(KK)=BKP+BJP
227       BJP=BKP-BJP
228       IF(ISN .LT. 0) GO TO 450
229       AKP=AKM-BJM
230       AKM=AKM+BJM
231       BKP=BKM+AJM
232       BKM=BKM-AJM
233       IF(S1 .EQ. 0.0) GO TO 460
234   430 A(K1)=AKP*C1-BKP*S1
235       B(K1)=AKP*S1+BKP*C1
236       A(K2)=AJP*C2-BJP*S2
237       B(K2)=AJP*S2+BJP*C2
238       A(K3)=AKM*C3-BKM*S3
239       B(K3)=AKM*S3+BKM*C3
240       KK=K3+KSPAN
241       IF(KK .LE. NT) GO TO 420
242   440 C2=C1-(CD*C1+SD*S1)
243       S1=(SD*C1-CD*S1)+S1
244 C     THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
245 C     ERROR. IF ROUNDED ARITHMETIC IS USED, THEY MAY BE DELETED.
246 C     C1=0.5/(C2**2+S1**2)+0.5
247 C     S1=C1*S1
248 C     C1=C1*C2
249 C     NEXT STATEMENT SHOULD BE DELETED IF NON-ROUNDED ARITHMETIC IS USED
250       C1=C2
251       C2=C1**2-S1**2
252       S2=2.0*C1*S1
253       C3=C2*C1-S2*S1
254       S3=C2*S1+S2*C1
255       KK=KK-NT+JC
256       IF(KK .LE. KSPAN) GO TO 420
257       KK=KK-KSPAN+INC
258       IF(KK .LE. JC) GO TO 410
259       IF(KSPAN .EQ. JC) GO TO 800
260       GO TO 100
261   450 AKP=AKM+BJM
262       AKM=AKM-BJM
263       BKP=BKM-AJM
264       BKM=BKM+AJM
265       IF(S1 .NE. 0.0) GO TO 430
266   460 A(K1)=AKP
267       B(K1)=BKP
268       A(K2)=AJP
269       B(K2)=BJP
270       A(K3)=AKM
271       B(K3)=BKM
272       KK=K3+KSPAN
273       IF(KK .LE. NT) GO TO 420
274       GO TO 440
275 C     TRANSFORM FOR FACTOR OF 5 (OPTIONAL CODE)
276   510 C2=C72**2-S72**2
277       S2=2.0*C72*S72
278   520 K1=KK+KSPAN
279       K2=K1+KSPAN
280       K3=K2+KSPAN
281       K4=K3+KSPAN
282       AKP=A(K1)+A(K4)
283       AKM=A(K1)-A(K4)
284       BKP=B(K1)+B(K4)
285       BKM=B(K1)-B(K4)
286       AJP=A(K2)+A(K3)
287       AJM=A(K2)-A(K3)
288       BJP=B(K2)+B(K3)
289       BJM=B(K2)-B(K3)
290       AA=A(KK)
291       BB=B(KK)
292       A(KK)=AA+AKP+AJP
293       B(KK)=BB+BKP+BJP
294       AK=AKP*C72+AJP*C2+AA
295       BK=BKP*C72+BJP*C2+BB
296       AJ=AKM*S72+AJM*S2
297       BJ=BKM*S72+BJM*S2
298       A(K1)=AK-BJ
299       A(K4)=AK+BJ
300       B(K1)=BK+AJ
301       B(K4)=BK-AJ
302       AK=AKP*C2+AJP*C72+AA
303       BK=BKP*C2+BJP*C72+BB
304       AJ=AKM*S2-AJM*S72
305       BJ=BKM*S2-BJM*S72
306       A(K2)=AK-BJ
307       A(K3)=AK+BJ
308       B(K2)=BK+AJ
309       B(K3)=BK-AJ
310       KK=K4+KSPAN
311       IF(KK .LT. NN) GO TO 520
312       KK=KK-NN
313       IF(KK .LE. KSPAN) GO TO 520
314       GO TO 700
315 C     TRANSFORM FOR ODD FACTORS
316   600 K=NFAC(I)
317       KSPNN=KSPAN
318       KSPAN=KSPAN/K
319       IF(K .EQ. 3) GO TO 320
320       IF(K .EQ. 5) GO TO 510
321       IF(K .EQ. JF) GO TO 640
322       JF=K
323       S1=RAD/K
324       C1=COS(S1)
325       S1=SIN(S1)
326       IF(JF .GT. MAXF) GO TO 998
327       CK(JF)=1.0
328       SK(JF)=0.0
329       J=1
330   630 CK(J)=CK(K)*C1+SK(K)*S1
331       SK(J)=CK(K)*S1-SK(K)*C1
332       K=K-1
333       CK(K)=CK(J)
334       SK(K)=-SK(J)
335       J=J+1
336       IF(J .LT. K) GO TO 630
337   640 K1=KK
338       K2=KK+KSPNN
339       AA=A(KK)
340       BB=B(KK)
341       AK=AA
342       BK=BB
343       J=1
344       K1=K1+KSPAN
345   650 K2=K2-KSPAN
346       J=J+1
347       AT(J)=A(K1)+A(K2)
348       AK=AT(J)+AK
349       BT(J)=B(K1)+B(K2)
350       BK=BT(J)+BK
351       J=J+1
352       AT(J)=A(K1)-A(K2)
353       BT(J)=B(K1)-B(K2)
354       K1=K1+KSPAN
355       IF(K1 .LT. K2) GO TO 650
356       A(KK)=AK
357       B(KK)=BK
358       K1=KK
359       K2=KK+KSPNN
360       J=1
361   660 K1=K1+KSPAN
362       K2=K2-KSPAN
363       JJ=J
364       AK=AA
365       BK=BB
366       AJ=0.0
367       BJ=0.0
368       K=1
369   670 K=K+1
370       AK=AT(K)*CK(JJ)+AK
371       BK=BT(K)*CK(JJ)+BK
372       K=K+1
373       AJ=AT(K)*SK(JJ)+AJ
374       BJ=BT(K)*SK(JJ)+BJ
375       JJ=JJ+J
376       IF(JJ .GT. JF) JJ=JJ-JF
377       IF(K .LT. JF) GO TO 670
378       K=JF-J
379       A(K1)=AK-BJ
380       B(K1)=BK+AJ
381       A(K2)=AK+BJ
382       B(K2)=BK-AJ
383       J=J+1
384       IF(J .LT. K) GO TO 660
385       KK=KK+KSPNN
386       IF(KK .LE. NN) GO TO 640
387       KK=KK-NN
388       IF(KK .LE. KSPAN) GO TO 640
389 C     MULTIPLY BY ROTATION FACTOR (EXCEPT FOR FACTORS OF 2 AND 4)
390   700 IF(I .EQ. M) GO TO 800
391       KK=JC+1
392   710 C2=1.0-CD
393       S1=SD
394   720 C1=C2
395       S2=S1
396       KK=KK+KSPAN
397   730 AK=A(KK)
398       A(KK)=C2*AK-S2*B(KK)
399       B(KK)=S2*AK+C2*B(KK)
400       KK=KK+KSPNN
401       IF(KK .LE. NT) GO TO 730
402       AK=S1*S2
403       S2=S1*C2+C1*S2
404       C2=C1*C2-AK
405       KK=KK-NT+KSPAN
406       IF(KK .LE. KSPNN) GO TO 730
407       C2=C1-(CD*C1+SD*S1)
408       S1=S1+(SD*C1-CD*S1)
409 C     THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
410 C     ERROR. IF ROUNDED ARITHMETIC IS USED, THEY MAY
411 C     BE DELETED.
412 C     C1=0.5/(C2**2+S1**2)+0.5
413 C     S1=C1*S1
414 C     C2=C1*C2
415       KK=KK-KSPNN+JC
416       IF(KK .LE. KSPAN) GO TO 720
417       KK=KK-KSPAN+JC+INC
418       IF(KK .LE. JC+JC) GO TO 710
419       GO TO 100
420 C     PERMUTE THE RESULTS TO NORMAL ORDER---DONE IN TWO STAGES
421 C     PERMUTATION FOR SQUARE FACTORS OF N
422   800 NP(1)=KS
423       IF(KT .EQ. 0) GO TO 890
424       K=KT+KT+1
425       IF(M .LT. K) K=K-1
426       J=1
427       NP(K+1)=JC
428   810 NP(J+1)=NP(J)/NFAC(J)
429       NP(K)=NP(K+1)*NFAC(J)
430       J=J+1
431       K=K-1
432       IF(J .LT. K) GO TO 810
433       K3=NP(K+1)
434       KSPAN=NP(2)
435       KK=JC+1
436       K2=KSPAN+1
437       J=1
438       IF(N .NE. NTOT) GO TO 850
439 C     PERMUTATION FOR SINGLE-VARIATE TRANSFORM (OPTIONAL CODE)
440   820 AK=A(KK)
441       A(KK)=A(K2)
442       A(K2)=AK
443       BK=B(KK)
444       B(KK)=B(K2)
445       B(K2)=BK
446       KK=KK+INC
447       K2=KSPAN+K2
448       IF(K2 .LT. KS) GO TO 820
449   830 K2=K2-NP(J)
450       J=J+1
451       K2=NP(J+1)+K2
452       IF(K2 .GT. NP(J)) GO TO 830
453       J=1
454   840 IF(KK .LT. K2) GO TO 820
455       KK=KK+INC
456       K2=KSPAN+K2
457       IF(K2 .LT. KS) GO TO 840
458       IF(KK .LT. KS) GO TO 830
459       JC=K3
460       GO TO 890
461 C     PERMUTATION FOR MULTIVARIATE TRANSFORM
462   850 K=KK+JC
463   860 AK=A(KK)
464       A(KK)=A(K2)
465       A(K2)=AK
466       BK=B(KK)
467       B(KK)=B(K2)
468       B(K2)=BK
469       KK=KK+INC
470       K2=K2+INC
471       IF(KK .LT. K) GO TO 860
472       KK=KK+KS-JC
473       K2=K2+KS-JC
474       IF(KK .LT. NT) GO TO 850
475       K2=K2-NT+KSPAN
476       KK=KK-NT+JC
477       IF(K2 .LT. KS) GO TO 850
478   870 K2=K2-NP(J)
479       J=J+1
480       K2=NP(J+1)+K2
481       IF(K2 .GT. NP(J)) GO TO 870
482       J=1
483   880 IF(KK .LT. K2) GO TO 850
484       KK=KK+JC
485       K2=KSPAN+K2
486       IF(K2 .LT. KS) GO TO 880
487       IF(KK .LT. KS) GO TO 870
488       JC=K3
489   890 IF(2*KT+1 .GE. M) RETURN
490       KSPNN=NP(KT+1)
491 C     PERMUTATION FOR SQUARE-FREE FACTORS OF N
492       J=M-KT
493       NFAC(J+1)=1
494   900 NFAC(J)=NFAC(J)*NFAC(J+1)
495       J=J-1
496       IF(J .NE. KT) GO TO 900
497       KT=KT+1
498       NN=NFAC(KT)-1
499       IF(NN .GT. MAXP) GO TO 998
500       JJ=0
501       J=0
502       GO TO 906
503   902 JJ=JJ-K2
504       K2=KK
505       K=K+1
506       KK=NFAC(K)
507   904 JJ=KK+JJ
508       IF(JJ .GE. K2) GO TO 902
509       NP(J)=JJ
510   906 K2=NFAC(KT)
511       K=KT+1
512       KK=NFAC(K)
513       J=J+1
514       IF(J .LE. NN) GO TO 904
515 C     DETERMINE THE PERMUTATION CYCLES OF LENGTH GREATER THAN 1
516       J=0
517       GO TO 914
518   910 K=KK
519       KK=NP(K)
520       NP(K)=-KK
521       IF(KK .NE. J) GO TO 910
522       K3=KK
523   914 J=J+1
524       KK=NP(J)
525       IF(KK .LT. 0) GO TO 914
526       IF(KK .NE. J) GO TO 910
527       NP(J)=-J
528       IF(J .NE. NN) GO TO 914
529       MAXF=INC*MAXF
530 C     REORDER A AND B, FOLLOWING THE PERMUTATION CYCLES
531       GO TO 950
532   924 J=J-1
533       IF(NP(J) .LT. 0) GO TO 924
534       JJ=JC
535   926 KSPAN=JJ
536       IF(JJ .GT. MAXF) KSPAN=MAXF
537       JJ=JJ-KSPAN
538       K=NP(J)
539       KK=JC*K+II+JJ
540       K1=KK+KSPAN
541       K2=0
542   928 K2=K2+1
543       AT(K2)=A(K1)
544       BT(K2)=B(K1)
545       K1=K1-INC
546       IF(K1 .NE. KK) GO TO 928
547   932 K1=KK+KSPAN
548       K2=K1-JC*(K+NP(K))
549       K=-NP(K)
550   936 A(K1)=A(K2)
551       B(K1)=B(K2)
552       K1=K1-INC
553       K2=K2-INC
554       IF(K1 .NE. KK) GO TO 936
555       KK=K2
556       IF(K .NE. J) GO TO 932
557       K1=KK+KSPAN
558       K2=0
559   940 K2=K2+1
560       A(K1)=AT(K2)
561       B(K1)=BT(K2)
562       K1=K1-INC
563       IF(K1 .NE. KK) GO TO 940
564       IF(JJ .NE. 0) GO TO 926
565       IF(J .NE. 1) GO TO 924
566   950 J=K3+1
567       NT=NT-KSPNN
568       II=NT-INC+1
569       IF(NT .GE. 0) GO TO 924
570       RETURN
571 C     ERROR FINISH, INSUFFICIENT ARRAY STORAGE
572   998 ISN=0
573       WRITE(6,999)
574       STOP
575   999 FORMAT('0ARRAY BOUNDS EXCEEDED WITHIN SUBROUTINE CFT')
576       END