program example C To use this distribution, clip off the main program and the 'BUILD.COM' C file. The remainder will be the FT code. You can fragment it further C if you want to do so. BTW: The BUILD.COM is for a VAX. I know. Yuck. C But this compiles just fine on Ultrix, and most everything else, so don't C get your jeans in an uproar. real a(10200) complex c(1) equivalence (a,c) C Forward 2d transform N = 100 M = 100 call real2f(a,N,M, +1) C This has created a 2-d transform, 102 x 100 reals (51 x 100 complex) C for an N x M array the transform is stored as a (N/2+1) x M complex C In X orders go from 0 to N/2 (0 to 50 in our example here: 50 is the C Nyquist order. In Y the orders go from 0 to M/2, and then down from C -M/2 to -1. You can reorder the array to get the zero order in the C middle, but I do not use it like that, and you are probably better C off leaving well enough alone adn indexing the array with something C like 'val at order (i,j) = c(i+1+(N+2)*mod(j+M,M))' C ** NOTE ** N returns as N/2+1 i.e. the X dimension of the complex array C generated C Reverse Transform call real2f(a,N,M, -1) call exit end C --------------------------- Clip Here ----------------------------- $for/ext CMPLFT.FOR $for/ext DIPRP.FOR $for/ext HERMFT.FOR $for/ext MDFTKD.FOR $for/ext MDP_ERROR.FOR ! <-- This is a dummy routine $for/ext R2CFTK.FOR $for/ext R3CFTK.FOR $for/ext R4CFTK.FOR $for/ext R5CFTK.FOR $for/ext R8CFTK.FOR $for/ext REAL2F.FOR $for/ext REAL3F.FOR $for/ext REALFT.FOR $for/ext RPCFTK.FOR $for/ext SRFP.FOR $library/create ten_eyke_fft.olb *.obj $Delete/noconfirm *.obj;* $write sys$output "FFT Library complete" $exit $! ------------------------- Clip here _________________________________ SUBROUTINE CMPLFT (X, Y, N, D) INTEGER N INTEGER D(5) REAL X(1), Y(1) C LOGICAL ERROR INTEGER P MAX, P SYM, TWO GRP INTEGER FACTOR(15), SYM(15), UN SYM(15) C P MAX = 41 TWO GRP = 8 C IF (N .LE. 1) GO TO 100 CALL S R FP (N, P MAX, TWO GRP, FACTOR, SYM, P SYM, UN SYM, ERROR) IF (ERROR) GO TO 200 C CALL MD FT KD (N, FACTOR, D, X, Y) CALL D IP RP (N, SYM, P SYM, UN SYM, D, X, Y) C 100 CONTINUE RETURN C 200 CONTINUE call mdp_error(13,0) return C c 300 FORMAT ('0INVALID NUMBER OF POINTS FOR CMPL FT. N =', I10//) C END SUBROUTINE DIPRP (PTS,SYM,P SYM,UN SYM,DIM,X,Y) C DOUBLE IN PLACE REORDERING PROGRAMME C INTEGER PTS,P SYM REAL X (4000), Y (4000) INTEGER SYM(15), UN SYM(15), DIM(5) INTEGER*4 K C REAL T LOGICAL ONE MOD INTEGER MODULO (14) INTEGER DK,JJ,KK,LK,MODS,MULT,P UN SYM,TEST INTEGER NT, SEP, P, P0, P1, P2, P3, P4, P5, SIZE C INTEGER S (14), U (14) INTEGER A,B,C,D,E,F,G,H,I,J,L,M,N INTEGER BS,CS,DS,ES,FS,GS,HS,IS,JS,KS,LS,MS,NS INTEGER AL,BL,CL,DL,EL,FL,GL,HL,IL,JL,KL,LL,ML,NL EQUIVALENCE (AL,U(1)),(BS,S(2)),(BL,U(2)) EQUIVALENCE (CS,S(3)),(CL,U(3)),(DS,S(4)),(DL,U(4)) EQUIVALENCE (ES,S(5)),(EL,U(5)),(FS,S(6)),(FL,U(6)) EQUIVALENCE (GS,S(7)),(GL,U(7)),(HS,S(8)),(HL,U(8)) EQUIVALENCE (IS,S(9)),(IL,U(9)),(JS,S(10)),(JL,U(10)) EQUIVALENCE (KS,S(11)),(KL,U(11)),(LS,S(12)),(LL,U(12)) EQUIVALENCE (MS,S(13)),(ML,U(13)),(NS,S(14)),(NL,U(14)) PARAMETER NEST = 14 C NT = DIM(1) SEP = DIM(2) P2 = DIM(3) SIZE = DIM(4) - 1 P4 = DIM(5) IF (SYM(1).EQ.0) GO TO 500 DO 100 J=1,NEST U(J)=1 S(J)=1 100 CONTINUE N=PTS DO 200 J=1,NEST IF (SYM(J).EQ.0) GO TO 300 JJ=NEST+1-J U(JJ)=N S(JJ)=N/SYM(J) N=N/SYM(J) 200 CONTINUE 300 CONTINUE C JJ=0 DO 400 A=1,AL DO 400 B=A,BL,BS DO 400 C=B,CL,CS DO 400 D=C,DL,DS DO 400 E=D,EL,ES DO 400 F=E,FL,FS DO 400 G=F,GL,GS DO 400 H=G,HL,HS DO 400 I=H,IL,IS DO 400 J=I,JL,JS DO 400 KK=J,KL,KS DO 400 L=KK,LL,LS DO 400 M=L,ML,MS DO 400 N=M,NL,NS JJ=JJ+1 IF (JJ.GE.N) GO TO 400 DELTA = (N-JJ)*SEP P1 = (JJ-1)*SEP + 1 DO 350 P0 = P1, NT, P2 P3 = P0 + SIZE DO 350 P = P0, P3, P4 P5 = P + DELTA T = X(P) X(P) = X(P5) X(P5) = T T = Y(P) Y(P) = Y(P5) Y(P5) = T 350 CONTINUE 400 CONTINUE C 500 CONTINUE IF (UN SYM(1).EQ.0) GO TO 1900 P UN SYM=PTS/P SYM**2 MULT=P UN SYM/UN SYM(1) TEST=(UN SYM(1)*UN SYM(2)-1)*MULT*P SYM LK=MULT DK=MULT DO 600 K=2,NEST IF (UN SYM(K).EQ.0) GO TO 700 LK=LK*UN SYM(K-1) DK=DK/UN SYM(K) U(K)=(LK-DK)*P SYM MODS=K 600 CONTINUE 700 CONTINUE ONE MOD=MODS.LT.3 IF (ONE MOD) GO TO 900 DO 800 J=3,MODS JJ=MODS+3-J MODULO(JJ)=U(J) 800 CONTINUE 900 CONTINUE MODULO(2)=U(2) JL=(P UN SYM-3)*P SYM MS=P UN SYM*P SYM C DO 1800 J=P SYM,JL,P SYM K=J C 1000 K=K*MULT IF (ONE MOD) GO TO 1200 DO 1100 I=3,MODS 1100 K=K-(K/MODULO(I))*MODULO(I) C 1200 IF (K.GE.TEST) GO TO 1300 K=K-(K/MODULO(2))*MODULO(2) GO TO 1400 1300 K=K-(K/MODULO(2))*MODULO(2)+MODULO(2) 1400 IF (K.LT.J) GO TO 1000 C IF (K.EQ.J) GO TO 1700 DELTA = (K-J)*SEP DO 1600 L=1,P SYM DO 1500 M=L,PTS,MS P1 = (M+J-1)*SEP + 1 DO 1500 P0 = P1, NT, P2 P3 = P0 + SIZE DO 1500 JJ = P0, P3, P4 KK = JJ + DELTA T=X(JJ) X(JJ)=X(KK) X(KK)=T T=Y(JJ) Y(JJ)=Y(KK) Y(KK)=T 1500 CONTINUE 1600 CONTINUE 1700 CONTINUE 1800 CONTINUE C 1900 CONTINUE RETURN END SUBROUTINE HERMFT (X, Y, N, DIM) INTEGER N INTEGER DIM(5) REAL X(1), Y(1) C REAL A, B, C, D, E, F, ANGLE, CO, SI, TWO N, TWO PI INTEGER NT, D2, D3, D4, D5 INTEGER I, J, K, N OVER 2, I0, I1, I2 INTEGER K1 C TWO PI = 6.283185 TWO N = FLOAT(2*N) C NT = DIM(1) D2 = DIM(2) D3 = DIM(3) D4 = DIM(4) - 1 D5 = DIM(5) C DO 100 I0 = 1, NT, D3 I1 = I0 + D4 DO 100 I = I0, I1, D5 A = X(I) B = Y(I) X(I) = A + B Y(I) = A - B 100 CONTINUE C N OVER 2 = N/2 + 1 IF (N OVER 2 .LT. 2) GO TO 500 DO 400 I0 = 2, N OVER 2 ANGLE = TWO PI*FLOAT(I0-1)/TWO N CO = COS(ANGLE) SI = SIN(ANGLE) K = (N + 2 - 2*I0)*D2 K1 = (I0 - 1)*D2 + 1 DO 300 I1 = K1, NT, D3 I2 = I1 + D4 DO 200 I = I1, I2, D5 J = I + K A = X(I) + X(J) B = X(I) - X(J) C = Y(I) + Y(J) D = Y(I) - Y(J) E = B*CO + C*SI F = B*SI - C*CO X(I) = A + F X(J) = A - F Y(I) = E + D Y(J) = E - D 200 CONTINUE 300 CONTINUE 400 CONTINUE C CALL CMPL FT (X, Y, N, DIM) C 500 CONTINUE RETURN C END SUBROUTINE MDFTKD (N, FACTOR, DIM, X, Y) C MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL DRIVER C INTEGER N INTEGER FACTOR (10), DIM(5) REAL X(10), Y(10) C INTEGER F, M, P, R, S C S = DIM(2) F = 0 M = N 100 CONTINUE F = F + 1 P = FACTOR(F) IF (P.EQ.0) RETURN M = M/P R = M*S IF (P.GT.8) GO TO 700 GO TO (100, 200, 300, 400, 500, 800, 700, 600), P GO TO 800 C 200 CONTINUE CALL R2 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), DIM) GO TO 100 C 300 CONTINUE CALL R3 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1) ., DIM) GO TO 100 C 400 CONTINUE CALL R4 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1) ., X(3*R+1), Y(3*R+1), DIM) GO TO 100 C 500 CONTINUE CALL R5 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1) ., X(3*R+1), Y(3*R+1), X(4*R+1), Y(4*R+1), DIM) GO TO 100 C 600 CONTINUE CALL R8 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1) ., X(3*R+1), Y(3*R+1), X(4*R+1), Y(4*R+1), X(5*R+1), Y(5*R+1), .X(6*R+1), Y(6*R+1), X(7*R+1), Y(7*R+1), DIM) GO TO 100 C 700 CONTINUE CALL RP CFTK (N, M, P, R, X, Y, DIM) GO TO 100 C 800 CONTINUE CALL MDP_ERROR(13,0) RETURN C 900 FORMAT ('0TRANSFER ERROR DETECTED IN MDFTKD'//) END subroutine mdp_error entry mdp_error_t return end C++******************************************************************* C C $$ R2CFTK: C C CALL R2CFTK(N,M,X0,Y0,X1,Y1,DIM) C N C M C X0 C Y0 C X1 C Y1 C DIM C C--******************************************************************* C C SUBROUTINE R2 CFTK (N, M, X0, Y0, X1, Y1, DIM) C RADIX 2 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C INTEGER N, M, DIM(5) REAL X0 (10), Y0 (10), X1 (10), Y1 (10) C LOGICAL FOLD,ZERO INTEGER J,K,K0,M2,M OVER 2 INTEGER K1, K2, KK, L, L1, MM2, NT, SIZE, SEP INTEGER NS REAL ANGLE,C,IS,IU,RS,RU,S,TWOPI REAL FJM1, FM2 DATA TWO PI/6.283185/ C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M2=M*2 FM2 = FLOAT(M2) M OVER 2=M/2+1 MM2 = SEP*M2 C FJM1 = -1.0 DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*SEP + 1 FJM1 = FJM1 + 1.0 ANGLE = TWO PI*FJM1/FM2 ZERO=ANGLE.EQ.0.0 IF (ZERO) GO TO 200 C=COS(ANGLE) S=SIN(ANGLE) GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*SEP + 1 C=-C 200 CONTINUE C DO 500 KK = K0, NS, MM2 DO 440 L = KK, NT, L1 K1 = L + SIZE DO 420 K = L, K1, K2 RS=X0(K)+X1(K) IS=Y0(K)+Y1(K) RU=X0(K)-X1(K) IU=Y0(K)-Y1(K) X0(K)=RS Y0(K)=IS IF (ZERO) GO TO 300 X1(K)=RU*C+IU*S Y1(K)=IU*C-RU*S GO TO 400 300 CONTINUE X1(K)=RU Y1(K)=IU 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C RETURN END C++******************************************************************* C C $$ R3CFTK: C C CALL R3CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,DIM) C N C M C X0 C Y0 C X1 C Y1 C X2 C Y2 C DIM C C--******************************************************************* C C SUBROUTINE R3 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, DIM) C RADIX 3 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C INTEGER N, M, DIM(5) REAL X0 (10), Y0 (10), X1 (10), Y1 (10), X2 (10), Y2 (10) C LOGICAL FOLD,ZERO INTEGER J,K,K0,M3,M OVER 2 INTEGER K1, K2, KK, L, L1, MM3, NT, SIZE, SEP INTEGER NS REAL ANGLE,A,B,C1,C2,S1,S2,T,TWOPI REAL I0,I1,I2,IA,IB,IS,R0,R1,R2,RA,RB,RS REAL FJM1, FM3 DATA TWO PI/6.283185/, A/-0.5/, B/0.86602540/ C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M3=M*3 FM3 = FLOAT(M3) MM3 = SEP*M3 M OVER 2=M/2+1 C FJM1 = -1.0 DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*SEP + 1 FJM1 = FJM1 + 1.0 ANGLE = TWO PI*FJM1/FM3 ZERO=ANGLE.EQ.0.0 IF (ZERO) GO TO 200 C1=COS(ANGLE) S1=SIN(ANGLE) C2=C1*C1-S1*S1 S2=S1*C1+C1*S1 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*SEP + 1 T=C1*A+S1*B S1=C1*B-S1*A C1=T T=C2*A-S2*B S2=-C2*B-S2*A C2=T 200 CONTINUE C DO 500 KK = K0, NS, MM3 DO 440 L = KK, NT, L1 K1 = L + SIZE DO 420 K = L, K1, K2 R0=X0(K) I0=Y0(K) RS=X1(K)+X2(K) IS=Y1(K)+Y2(K) X0(K)=R0+RS Y0(K)=I0+IS RA=R0+RS*A IA=I0+IS*A RB=(X1(K)-X2(K))*B IB=(Y1(K)-Y2(K))*B IF (ZERO) GO TO 300 R1=RA+IB I1=IA-RB R2=RA-IB I2=IA+RB X1(K)=R1*C1+I1*S1 Y1(K)=I1*C1-R1*S1 X2(K)=R2*C2+I2*S2 Y2(K)=I2*C2-R2*S2 GO TO 400 300 CONTINUE X1(K)=RA+IB Y1(K)=IA-RB X2(K)=RA-IB Y2(K)=IA+RB 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C RETURN END C++******************************************************************* C C $$ R4CFTK: C C CALL R4CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,DIM) C N C M C X0 C Y0 C X1 C Y1 C X2 C Y2 C X3 C Y3 C DIM C C--******************************************************************* C C SUBROUTINE R4 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, DIM) C RADIX 4 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C INTEGER N, M, DIM(5) REAL X0 (10), Y0 (10), X1 (10), Y1 (10) REAL X2 (10), Y2 (10), X3 (10), Y3 (10) C LOGICAL FOLD,ZERO INTEGER J,K,K0,M4,M OVER 2 INTEGER K1, K2, KK, L, L1, MM4, NT, SIZE, SEP INTEGER NS REAL ANGLE,C1,C2,C3,S1,S2,S3,T,TWOPI REAL I1,I2,I3,IS0,IS1,IU0,IU1,R1,R2,R3,RS0,RS1,RU0,RU1 REAL FJM1, FM4 DATA TWO PI/6.283185/ C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M4=M*4 FM4 = FLOAT(M4) MM4 = SEP*M4 M OVER 2=M/2+1 C FJM1 = -1.0 DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*SEP + 1 FJM1 = FJM1 + 1.0 ANGLE = TWO PI*FJM1/FM4 ZERO=ANGLE.EQ.0.0 IF (ZERO) GO TO 200 C1=COS(ANGLE) S1=SIN(ANGLE) C2=C1*C1-S1*S1 S2=S1*C1+C1*S1 C3=C2*C1-S2*S1 S3=S2*C1+C2*S1 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*SEP + 1 T=C1 C1=S1 S1=T C2=-C2 T=C3 C3=-S3 S3=-T 200 CONTINUE C DO 500 KK = K0, NS, MM4 DO 440 L = KK, NT, L1 K1 = L + SIZE DO 420 K = L, K1, K2 RS0=X0(K)+X2(K) IS0=Y0(K)+Y2(K) RU0=X0(K)-X2(K) IU0=Y0(K)-Y2(K) RS1=X1(K)+X3(K) IS1=Y1(K)+Y3(K) RU1=X1(K)-X3(K) IU1=Y1(K)-Y3(K) X0(K)=RS0+RS1 Y0(K)=IS0+IS1 IF (ZERO) GO TO 300 R1=RU0+IU1 I1=IU0-RU1 R2=RS0-RS1 I2=IS0-IS1 R3=RU0-IU1 I3=IU0+RU1 X2(K)=R1*C1+I1*S1 Y2(K)=I1*C1-R1*S1 X1(K)=R2*C2+I2*S2 Y1(K)=I2*C2-R2*S2 X3(K)=R3*C3+I3*S3 Y3(K)=I3*C3-R3*S3 GO TO 400 300 CONTINUE X2(K)=RU0+IU1 Y2(K)=IU0-RU1 X1(K)=RS0-RS1 Y1(K)=IS0-IS1 X3(K)=RU0-IU1 Y3(K)=IU0+RU1 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C RETURN END C++******************************************************************* C C $$ R5CFTK: C C CALL R5CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4,DIM) C N C M C X0 C Y0 C X1 C Y1 C X2 C Y2 C X3 C Y3 C X4 C Y4 C DIM C C--******************************************************************* C C SUBROUTINE R5 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, X4, Y4, . DIM) C RADIX 5 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C INTEGER N, M, DIM(5) REAL X0 (10), Y0 (10), X1 (10), Y1 (10), X2 (10), Y2 (10) REAL X3 (10), Y3 (10), X4 (10), Y4 (10) C LOGICAL FOLD,ZERO INTEGER J,K,K0,M5,M OVER 2 INTEGER K1, K2, KK, L, L1, MM5, NT, SIZE, SEP INTEGER NS REAL ANGLE,A1,A2,B1,B2,C1,C2,C3,C4,S1,S2,S3,S4,T,TWOPI REAL R0,R1,R2,R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2 REAL I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2 REAL FJM1, FM5 DATA TWO PI/6.283185/, A1/0.30901699/, B1/0.95105652/, . A2/-0.80901699/, B2/0.58778525/ C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M5=M*5 FM5 = FLOAT(M5) MM5 = SEP*M5 M OVER 2=M/2+1 C FJM1 = -1.0 DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*SEP + 1 FJM1 = FJM1 + 1.0 ANGLE = TWO PI*FJM1/FM5 ZERO=ANGLE.EQ.0.0 IF (ZERO) GO TO 200 C1=COS(ANGLE) S1=SIN(ANGLE) C2=C1*C1-S1*S1 S2=S1*C1+C1*S1 C3=C2*C1-S2*S1 S3=S2*C1+C2*S1 C4=C2*C2-S2*S2 S4=S2*C2+C2*S2 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*SEP + 1 T=C1*A1+S1*B1 S1=C1*B1-S1*A1 C1=T T=C2*A2+S2*B2 S2=C2*B2-S2*A2 C2=T T=C3*A2-S3*B2 S3=-C3*B2-S3*A2 C3=T T=C4*A1-S4*B1 S4=-C4*B1-S4*A1 C4=T 200 CONTINUE C DO 500 KK = K0, NS, MM5 DO 440 L = KK, NT, L1 K1 = L + SIZE DO 420 K = L, K1, K2 R0=X0(K) I0=Y0(K) RS1=X1(K)+X4(K) IS1=Y1(K)+Y4(K) RU1=X1(K)-X4(K) IU1=Y1(K)-Y4(K) RS2=X2(K)+X3(K) IS2=Y2(K)+Y3(K) RU2=X2(K)-X3(K) IU2=Y2(K)-Y3(K) X0(K)=R0+RS1+RS2 Y0(K)=I0+IS1+IS2 RA1=R0+RS1*A1+RS2*A2 IA1=I0+IS1*A1+IS2*A2 RA2=R0+RS1*A2+RS2*A1 IA2=I0+IS1*A2+IS2*A1 RB1=RU1*B1+RU2*B2 IB1=IU1*B1+IU2*B2 RB2=RU1*B2-RU2*B1 IB2=IU1*B2-IU2*B1 IF (ZERO) GO TO 300 R1=RA1+IB1 I1=IA1-RB1 R2=RA2+IB2 I2=IA2-RB2 R3=RA2-IB2 I3=IA2+RB2 R4=RA1-IB1 I4=IA1+RB1 X1(K)=R1*C1+I1*S1 Y1(K)=I1*C1-R1*S1 X2(K)=R2*C2+I2*S2 Y2(K)=I2*C2-R2*S2 X3(K)=R3*C3+I3*S3 Y3(K)=I3*C3-R3*S3 X4(K)=R4*C4+I4*S4 Y4(K)=I4*C4-R4*S4 GO TO 400 300 CONTINUE X1(K)=RA1+IB1 Y1(K)=IA1-RB1 X2(K)=RA2+IB2 Y2(K)=IA2-RB2 X3(K)=RA2-IB2 Y3(K)=IA2+RB2 X4(K)=RA1-IB1 Y4(K)=IA1+RB1 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C RETURN END C++******************************************************************* C C $$ R8CFTK: C C CALL R8CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4,X5,Y5, C X6,Y6,X7,Y7,DIM) C C N C M C X0 C Y0 C X1 C Y1 C X2 C Y2 C X3 C Y3 C X4 C Y4 C X5 C Y5 C X6 C Y6 C X7 C Y7 C DIM C C--******************************************************************* C C SUBROUTINE R8 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, X4, Y4, . X5, Y5, X6, Y6, X7, Y7, DIM) C RADIX 8 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C INTEGER N, M, DIM(5) REAL X0 (10), Y0 (10), X1 (10), Y1 (10), X2 (10), Y2 (10) REAL X3 (10), Y3 (10), X4 (10), Y4 (10), X5 (10), Y5 (10) REAL X6 (10), Y6 (10), X7 (10), Y7 (10) C LOGICAL FOLD,ZERO INTEGER J,K,K0,M8,M OVER 2 INTEGER K1, K2, KK, L, L1, MM8, NT, SIZE, SEP INTEGER NS REAL ANGLE,C1,C2,C3,C4,C5,C6,C7,E,S1,S2,S3,S4,S5,S6,S7,T,TWOPI REAL R1,R2,R3,R4,R5,R6,R7,RS0,RS1,RS2,RS3,RU0,RU1,RU2,RU3 REAL I1,I2,I3,I4,I5,I6,I7,IS0,IS1,IS2,IS3,IU0,IU1,IU2,IU3 REAL RSS0,RSS1,RSU0,RSU1,RUS0,RUS1,RUU0,RUU1 REAL ISS0,ISS1,ISU0,ISU1,IUS0,IUS1,IUU0,IUU1 REAL FJM1, FM8 DATA TWO PI/6.283185/, E/0.70710678/ C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M8=M*8 FM8 = FLOAT(M8) MM8 = SEP*M8 M OVER 2=M/2+1 C FJM1 = -1.0 DO 600 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*SEP + 1 FJM1 = FJM1 + 1.0 ANGLE = TWO PI*FJM1/FM8 ZERO=ANGLE.EQ.0.0 IF (ZERO) GO TO 200 C1=COS(ANGLE) S1=SIN(ANGLE) C2=C1*C1-S1*S1 S2=S1*C1+C1*S1 C3=C2*C1-S2*S1 S3=S2*C1+C2*S1 C4=C2*C2-S2*S2 S4=S2*C2+C2*S2 C5=C4*C1-S4*S1 S5=S4*C1+C4*S1 C6=C4*C2-S4*S2 S6=S4*C2+C4*S2 C7=C4*C3-S4*S3 S7=S4*C3+C4*S3 GO TO 200 100 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*SEP + 1 T=(C1+S1)*E S1=(C1-S1)*E C1=T T=S2 S2=C2 C2=T T=(-C3+S3)*E S3=(C3+S3)*E C3=T C4=-C4 T=-(C5+S5)*E S5=(-C5+S5)*E C5=T T=-S6 S6=-C6 C6=T T=(C7-S7)*E S7=-(C7+S7)*E C7=T 200 CONTINUE C DO 500 KK = K0, NS, MM8 DO 440 L = KK, NT, L1 K1 = L + SIZE DO 420 K = L, K1, K2 RS0=X0(K)+X4(K) IS0=Y0(K)+Y4(K) RU0=X0(K)-X4(K) IU0=Y0(K)-Y4(K) RS1=X1(K)+X5(K) IS1=Y1(K)+Y5(K) RU1=X1(K)-X5(K) IU1=Y1(K)-Y5(K) RS2=X2(K)+X6(K) IS2=Y2(K)+Y6(K) RU2=X2(K)-X6(K) IU2=Y2(K)-Y6(K) RS3=X3(K)+X7(K) IS3=Y3(K)+Y7(K) RU3=X3(K)-X7(K) IU3=Y3(K)-Y7(K) RSS0=RS0+RS2 ISS0=IS0+IS2 RSU0=RS0-RS2 ISU0=IS0-IS2 RSS1=RS1+RS3 ISS1=IS1+IS3 RSU1=RS1-RS3 ISU1=IS1-IS3 RUS0=RU0-IU2 IUS0=IU0+RU2 RUU0=RU0+IU2 IUU0=IU0-RU2 RUS1=RU1-IU3 IUS1=IU1+RU3 RUU1=RU1+IU3 IUU1=IU1-RU3 T=(RUS1+IUS1)*E IUS1=(IUS1-RUS1)*E RUS1=T T=(RUU1+IUU1)*E IUU1=(IUU1-RUU1)*E RUU1=T X0(K)=RSS0+RSS1 Y0(K)=ISS0+ISS1 IF (ZERO) GO TO 300 R1=RUU0+RUU1 I1=IUU0+IUU1 R2=RSU0+ISU1 I2=ISU0-RSU1 R3=RUS0+IUS1 I3=IUS0-RUS1 R4=RSS0-RSS1 I4=ISS0-ISS1 R5=RUU0-RUU1 I5=IUU0-IUU1 R6=RSU0-ISU1 I6=ISU0+RSU1 R7=RUS0-IUS1 I7=IUS0+RUS1 X4(K)=R1*C1+I1*S1 Y4(K)=I1*C1-R1*S1 X2(K)=R2*C2+I2*S2 Y2(K)=I2*C2-R2*S2 X6(K)=R3*C3+I3*S3 Y6(K)=I3*C3-R3*S3 X1(K)=R4*C4+I4*S4 Y1(K)=I4*C4-R4*S4 X5(K)=R5*C5+I5*S5 Y5(K)=I5*C5-R5*S5 X3(K)=R6*C6+I6*S6 Y3(K)=I6*C6-R6*S6 X7(K)=R7*C7+I7*S7 Y7(K)=I7*C7-R7*S7 GO TO 400 300 CONTINUE X4(K)=RUU0+RUU1 Y4(K)=IUU0+IUU1 X2(K)=RSU0+ISU1 Y2(K)=ISU0-RSU1 X6(K)=RUS0+IUS1 Y6(K)=IUS0-RUS1 X1(K)=RSS0-RSS1 Y1(K)=ISS0-ISS1 X5(K)=RUU0-RUU1 Y5(K)=IUU0-IUU1 X3(K)=RSU0-ISU1 Y3(K)=ISU0+RSU1 X7(K)=RUS0-IUS1 Y7(K)=IUS0+RUS1 400 CONTINUE 420 CONTINUE 440 CONTINUE 500 CONTINUE IF (FOLD) GO TO 100 600 CONTINUE C RETURN END SUBROUTINE REAL2F(A,NX,NY,JD) C ************************************** Copyright 1983 P.R. Smith. INTEGER ID(5) REAL A(1) C C TWO DIMENSIONAL TRANSFORM OF REAL OBJECT C IF(.NOT.(MOD(NX,2).EQ.0.OR.JD.LT.0))THEN CALL MDP_ERROR(13,0) RETURN END IF C C SWITCH TO FORWARD OR REVERSE TRANSFORM C IF(JD)1,2,2 C C REAL TRANSFORM C 2 NXX=NX+2 NS=NXX*NY C C PACK OUT ROWS OF 'A' WITH ZEROS C DO I=NY,1,-1 L=I LF=NXX*(L-1) LN=NX*(L-1) DO J=NX,1,-1 A(LF+J)=A(LN+J) END DO END DO C C NOW PERFORM THE REAL F.T. AS STEP ONE NT=NX/2 ID(1)=NS ID(2)=2 ID(3)=NS ID(4)=NS ID(5)=NXX C CALL REALFT(A(1),A(2),NT,ID) C C NOW THE COMPLEX F.T. C ID(1)=NS ID(2)=NXX ID(3)=NS ID(4)=NXX ID(5)=2 C CALL CMPLFT(A(1),A(2),NY,ID) C NX=NT+1 RETURN C C HERMITIAN SYMMETRIC TRANSFORM C 1 CONTINUE NXX=2*NX NT=NXX-2 NS=NY*NXX C C COMPLEX TRANSFORM C ID(1)=NS ID(2)=NXX ID(3)=NS ID(4)=NXX ID(5)=2 C C CONJUGATE ARRAY AND TRANSFORM C DO10I=2,NS,2 10 A(I)=-A(I) C CALL CMPLFT(A(1),A(2),NY,ID) C C COPY LAST ELEMENT OVER AND COMPRESS C DO12I=1,NY LY=NXX*(I-1) LM=NT*(I-1) A(2+LY)=A(NXX-1+LY) DO13J=1,NT 13 A(LM+J)=A(LY+J) 12 CONTINUE C C NOW THE HERMITIAN F.T. C NS=NY*NT ID(1)=NS ID(2)=2 ID(3)=NS ID(4)=NS ID(5)=NT C CALL HERMFT(A(1),A(2),NT/2,ID) C NX=NT RETURN END SUBROUTINE REAL3F(A,NX,NY,NZ,JD) C ************************************** Copyright 1990 P.R. Smith. INTEGER ID(5) REAL A(*) C C THREE DIMENSIONAL TRANSFORM OF REAL OBJECT C IF(.NOT.(MOD(NX,2).EQ.0.OR.JD.LT.0))THEN CALL MDP_ERROR(13,0) RETURN END IF C C SWITCH TO FORWARD OR REVERSE TRANSFORM C IF(JD)1,2,2 C C REAL TRANSFORM C 2 NXX=NX+2 NS=NXX*NY*NZ NXY=NXX*NY NYNZ=NY*NZ C C PACK OUT ROWS OF 'A' WITH ZEROS C DO I=NYNZ,1,-1 L=I LF=NXX*(L-1) LN=NX*(L-1) DO J=NX,1,-1 A(LF+J)=A(LN+J) END DO END DO C C NOW PERFORM THE REAL F.T. AS STEP ONE C NT=NX/2 ID(1)=NS ID(2)=2 ID(3)=NS ID(4)=NS ID(5)=NXX C CALL REALFT(A(1),A(2),NT,ID) C C NOW THE COMPLEX F.T. IN "Y" C ID(1)=NXY ID(2)=NXX ID(3)=NXY ID(4)=NXX ID(5)=2 C DO I=1,NZ NZ1=NXY*(I-1)+1 CALL CMPLFT(A(NZ1),A(NZ1+1),NY,ID) END DO C C NOW THE COMPLEX F.T. IN "Z" C ID(1)=NS ID(2)=NXY ID(3)=NS ID(4)=NXY ID(5)=2 C CALL CMPLFT(A(1),A(2),NZ,ID) C NX=NT+1 RETURN C C HERMITIAN SYMMETRIC TRANSFORM C 1 CONTINUE NXX=2*NX NT=NXX-2 NS=NY*NXX*NZ NXY=NXX*NY C C COMPLEX TRANSFORM "Z" C ID(1)=NS ID(2)=NXY ID(3)=NS ID(4)=NXY ID(5)=2 C C CONJUGATE ARRAY AND TRANSFORM C DO I=2,NS,2 A(I)=-A(I) END DO C CALL CMPLFT(A(1),A(2),NZ,ID) C C COMPLEX TRANSFORM "Y" C ID(1)=NXY ID(2)=NXX ID(3)=NXY ID(4)=NXX ID(5)=2 C DO I=1,NZ NZ1=NXY*(I-1)+1 CALL CMPLFT(A(NZ1),A(NZ1+1),NY,ID) END DO C C COPY LAST ELEMENT OVER AND COMPRESS C DO I=1,NY*NZ LY=NXX*(I-1) LM=NT*(I-1) A(2+LY)=A(NXX-1+LY) DO J=1,NT A(LM+J)=A(LY+J) END DO END DO C C NOW THE HERMITIAN F.T. C NS=NY*NT*NZ ID(1)=NS ID(2)=2 ID(3)=NS ID(4)=NS ID(5)=NT C CALL HERMFT(A(1),A(2),NT/2,ID) C NX=NT RETURN END SUBROUTINE REALFT (EVEN, ODD, N, DIM) INTEGER N INTEGER DIM(5) REAL EVEN(1), ODD(1) C REAL A, B, C, D, E, F, ANGLE, CO, SI, TWO PI, TWO N INTEGER NT, D2, D3, D4, D5 INTEGER I, J, K, L, N OVER 2, I0, I1, I2 C TWO PI = 6.283185 TWO N = FLOAT(2*N) C CALL CMPL FT (EVEN, ODD, N, DIM) C NT = DIM(1) D2 = DIM(2) D3 = DIM(3) D4 = DIM(4) - 1 D5 = DIM(5) N OVER 2 = N/2 + 1 C IF (N OVER 2 .LT. 2) GO TO 400 DO 300 I = 2, N OVER 2 ANGLE = TWO PI*FLOAT(I-1)/TWO N CO = COS(ANGLE) SI = SIN(ANGLE) I0 = (I - 1)*D2 + 1 J = (N + 2 - 2*I)*D2 DO 200 I1 = I0, NT, D3 I2 = I1 + D4 DO 100 K = I1, I2, D5 L = K + J A = (EVEN(L) + EVEN(K))/2.0 C = (EVEN(L) - EVEN(K))/2.0 B = (ODD(L) + ODD(K))/2.0 D = (ODD(L) - ODD(K))/2.0 E = C*SI + B*CO F = C*CO - B*SI EVEN(K) = A + E EVEN(L) = A - E ODD(K) = F - D ODD(L) = F + D 100 CONTINUE 200 CONTINUE 300 CONTINUE C 400 CONTINUE IF (N .LT. 1) GO TO 600 J = N*D2 DO 500 I1 = 1, NT, D3 I2 = I1 + D4 DO 500 K = I1, I2, D5 L = K + J EVEN(L) = EVEN(K) - ODD(K) ODD(L) = 0.0 EVEN(K) = EVEN(K) + ODD(K) ODD(K) = 0.0 500 CONTINUE C 600 CONTINUE RETURN END SUBROUTINE RPCFTK (N, M, P, R, X, Y, DIM) C RADIX PRIME MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL C INTEGER N, M, P, R, DIM(5) REAL X(R,P), Y(R,P) C LOGICAL FOLD,ZERO REAL ANGLE,IS,IU,RS,RU,T,TWOPI,XT,YT REAL FU, FP, FJM1, FMP INTEGER J,JJ,K0,K,M OVER 2,MP,PM,PP,U,V INTEGER K1, K2, KK, L, L1, MMP, NT, SIZE, SEP INTEGER NS C REAL AA (20,20), BB (20,20) REAL A (40), B (40), C (40), S (40) REAL IA (20), IB (20), RA (20), RB (20) DATA TWO PI/6.283185/ C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M OVER 2=M/2+1 MP=M*P FMP = FLOAT(MP) MMP = SEP*MP PP=P/2 PM=P-1 FP = FLOAT(P) FU = 0.0 DO 100 U=1,PP FU = FU + 1.0 ANGLE = TWO PI*FU/FP JJ=P-U A(U)=COS(ANGLE) B(U)=SIN(ANGLE) A(JJ)=A(U) B(JJ)=-B(U) 100 CONTINUE DO 300 U=1,PP DO 200 V=1,PP JJ=U*V-U*V/P*P AA(V,U)=A(JJ) BB(V,U)=B(JJ) 200 CONTINUE 300 CONTINUE C FJM1 = -1.0 DO 1500 J=1,M OVER 2 FOLD=J.GT.1 .AND. 2*J.LT.M+2 K0 = (J-1)*SEP + 1 FJM1 = FJM1 + 1.0 ANGLE = TWO PI*FJM1/FMP ZERO=ANGLE.EQ.0.0 IF (ZERO) GO TO 700 C(1)=COS(ANGLE) S(1)=SIN(ANGLE) DO 400 U=2,PM C(U)=C(U-1)*C(1)-S(U-1)*S(1) S(U)=S(U-1)*C(1)+C(U-1)*S(1) 400 CONTINUE GO TO 700 500 CONTINUE FOLD=.FALSE. K0 = (M+1-J)*SEP + 1 DO 600 U=1,PM T=C(U)*A(U)+S(U)*B(U) S(U)=-S(U)*A(U)+C(U)*B(U) C(U)=T 600 CONTINUE 700 CONTINUE C DO 1400 KK = K0, NS, MMP DO 1340 L = KK, NT, L1 K1 = L + SIZE DO 1320 K = L, K1, K2 XT=X(K,1) YT=Y(K,1) RS=X(K,2)+X(K,P) IS=Y(K,2)+Y(K,P) RU=X(K,2)-X(K,P) IU=Y(K,2)-Y(K,P) DO 800 U=1,PP RA(U)=XT+RS*AA(U,1) IA(U)=YT+IS*AA(U,1) RB(U)=RU*BB(U,1) IB(U)=IU*BB(U,1) 800 CONTINUE XT=XT+RS YT=YT+IS DO 1000 U=2,PP JJ=P-U RS=X(K,U+1)+X(K,JJ+1) IS=Y(K,U+1)+Y(K,JJ+1) RU=X(K,U+1)-X(K,JJ+1) IU=Y(K,U+1)-Y(K,JJ+1) XT=XT+RS YT=YT+IS DO 900 V=1,PP RA(V)=RA(V)+RS*AA(V,U) IA(V)=IA(V)+IS*AA(V,U) RB(V)=RB(V)+RU*BB(V,U) IB(V)=IB(V)+IU*BB(V,U) 900 CONTINUE 1000 CONTINUE X(K,1)=XT Y(K,1)=YT DO 1300 U=1,PP JJ=P-U IF (ZERO) GO TO 1100 XT=RA(U)+IB(U) YT=IA(U)-RB(U) X(K,U+1)=XT*C(U)+YT*S(U) Y(K,U+1)=YT*C(U)-XT*S(U) XT=RA(U)-IB(U) YT=IA(U)+RB(U) X(K,JJ+1)=XT*C(JJ)+YT*S(JJ) Y(K,JJ+1)=YT*C(JJ)-XT*S(JJ) GO TO 1200 1100 CONTINUE X(K,U+1)=RA(U)+IB(U) Y(K,U+1)=IA(U)-RB(U) X(K,JJ+1)=RA(U)-IB(U) Y(K,JJ+1)=IA(U)+RB(U) 1200 CONTINUE 1300 CONTINUE 1320 CONTINUE 1340 CONTINUE 1400 CONTINUE IF (FOLD) GO TO 500 1500 CONTINUE C RETURN END SUBROUTINE SRFP (PTS,PMAX,TWO GRP,FACTOR,SYM,P SYM,UN SYM,ERROR) C SYMMETRIZED REORDERING FACTORING PROGRAMME C LOGICAL ERROR INTEGER PTS,PMAX,TWO GRP,P SYM INTEGER FACTOR (10), SYM (10), UN SYM (10) C INTEGER PP(14), QQ (7) INTEGER F,J,JJ,N,NEST,P,P TWO,Q,R C NEST=14 C N=PTS P SYM=1 F=2 P=0 Q=0 100 CONTINUE IF (N.LE.1) GO TO 500 DO 200 J=F,PMAX IF (N.EQ.(N/J)*J) GO TO 300 200 CONTINUE GO TO 1400 300 CONTINUE IF (2*P+Q.GE.NEST) GO TO 1600 F=J N=N/F IF (N.EQ.(N/F)*F) GO TO 400 Q=Q+1 QQ(Q)=F GO TO 100 400 CONTINUE N=N/F P=P+1 PP(P)=F P SYM=P SYM*F GO TO 100 C 500 CONTINUE R=1 IF (Q.EQ.0) R=0 IF (P.LT.1) GO TO 700 DO 600 J=1,P JJ=P+1-J SYM(J)=PP(JJ) FACTOR(J)=PP(JJ) JJ=P+Q+J FACTOR(JJ)=PP(J) JJ=P+R+J SYM(JJ)=PP(J) 600 CONTINUE 700 CONTINUE IF (Q.LT.1) GO TO 900 DO 800 J=1,Q JJ=P+J UN SYM(J)=QQ(J) FACTOR(JJ)=QQ(J) 800 CONTINUE SYM(P+1)=PTS/P SYM**2 900 CONTINUE JJ=2*P+Q FACTOR(JJ+1)=0 P TWO=1 J=0 1000 CONTINUE J=J+1 IF (FACTOR(J).EQ.0) GO TO 1200 IF (FACTOR(J).NE.2) GO TO 1000 P TWO=P TWO*2 FACTOR(J)=1 IF (P TWO.GE.TWO GRP) GO TO 1100 IF (FACTOR(J+1).EQ.2) GO TO 1000 1100 CONTINUE FACTOR(J)=P TWO P TWO=1 GO TO 1000 1200 CONTINUE IF (P.EQ.0) R=0 JJ=2*P+R SYM(JJ+1)=0 IF (Q.LE.1) Q=0 UN SYM(Q+1)=0 ERROR=.FALSE. C 1300 CONTINUE RETURN C 1400 CONTINUE call mdp_error_t(56,'Largest factor exceeds maximum') ERROR=.TRUE. GO TO 1300 C 1500 FORMAT (24H0LARGEST FACTOR EXCEEDS ,I3,7H. N = ,I6,1H.) C 1600 CONTINUE call mdp_error_t(56,'Factor count exceeded') ERROR=.TRUE. GO TO 1300 C 1700 FORMAT (22H0FACTOR COUNT EXCEEDS ,I3,7H. N = ,I6,1H.) END