C C fftlib.f: fast-fourier transform library routines C Copyright (C) Lynn Ten Eyck C C This library is free software: you can redistribute it and/or C modify it under the terms of the GNU Lesser General Public License C version 3, modified in accordance with the provisions of the C license to address the requirements of UK law. C C You should have received a copy of the modified GNU Lesser General C Public License along with this library. If not, copies may be C downloaded from http://www.ccp4.ac.uk/ccp4license.php C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU Lesser General Public License for more details. C C-- FFT81 F77TRNFM.FOR 13/09/85 JWC C C C**** FOLLOWING ARE ROUTINES USED BY TEN EYCK'S FFT PROGRAMS*** C C SUBROUTINE CMPLFT(X,Y,N,D) C =============================== C C C Complex finite discrete fourier transform C transforms one dimension of multi-dimensional data C modified by L. F. TEN EYCK from a one-dimensional version written C by G. T. SANDE, 1969. C C This program calculates the transform C (X(T) + I*Y(T))*(COS(2*PI*T/N) - I*SIN(2*PI*T/N)) C C INDEXING -- the arrangement of the multi-dimensional data is C specified by the integer array D, the values of which are used as C control parameters in do loops. when it is desired to cover all C elements of the data for which the subscript being transformed has C the value I0, the following is used. C C I1 = (I0 - 1)*D(2) + 1 C DO 100 I2 = I1, D(1), D(3) C I3 = I2 + D(4) - 1 C DO 100 I = I2, I3, D(5) C . C . C 100 CONTINUE C C with this indexing it is possible to use a number of arrangements C of the data, including normal fortran complex numbers (d(5) = 2) C or separate storage of real and imaginary parts. C C C---- PMAX is the largest prime factor that will be tolerated by this C program. C C---- TWOGRP is the largest power of two that is treated as a special C case. C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL X(*),Y(*) INTEGER D(5) C .. C .. Local Scalars .. INTEGER PMAX,PSYM,TWOGRP LOGICAL ERROR CHARACTER EMESS*80 C .. C .. Local Arrays .. INTEGER FACTOR(15),SYM(15),UNSYM(15) C .. C .. External Subroutines .. EXTERNAL DIPRP,MDFTKD,SRFP C .. PMAX = 19 TWOGRP = 8 C IF (N.GT.1) THEN CALL SRFP(N,PMAX,TWOGRP,FACTOR,SYM,PSYM,UNSYM,ERROR) IF (ERROR) THEN C WRITE (EMESS,FMT=9000) N call ccperr(1, EMESS) ELSE C CALL MDFTKD(N,FACTOR,D,X,Y) CALL DIPRP(N,SYM,PSYM,UNSYM,D,X,Y) END IF END IF C C---- Format statements C 9000 FORMAT ('FFTLIB: invalid number of points for CMPL FT. N =',I10) END C C SUBROUTINE SRFP(PTS,PMAX,TWOGRP,FACTOR,SYM,PSYM,UNSYM,ERROR) C ================================================================== C C C---- Symmetrized reordering factoring programme C C C .. Scalar Arguments .. INTEGER PMAX,PSYM,PTS,TWOGRP LOGICAL ERROR C .. C .. Array Arguments .. INTEGER FACTOR(15),SYM(15),UNSYM(15) C .. C .. Local Scalars .. INTEGER F,J,JJ,N,NEST,P,PTWO,Q,R C .. C .. Local Arrays .. INTEGER PP(14),QQ(7) C .. NEST = 14 C N = PTS PSYM = 1 F = 2 P = 0 Q = 0 10 CONTINUE IF (N.LE.1) THEN GO TO 60 ELSE DO 20 J = F,PMAX IF (N.EQ. (N/J)*J) GO TO 30 20 CONTINUE GO TO 40 30 IF (2*P+Q.GE.NEST) THEN GO TO 50 ELSE F = J N = N/F IF (N.EQ. (N/F)*F) THEN N = N/F P = P + 1 PP(P) = F PSYM = PSYM*F ELSE Q = Q + 1 QQ(Q) = F END IF GO TO 10 END IF END IF C 40 CONTINUE WRITE (6,FMT=9000) PMAX,PTS ERROR = .TRUE. GO TO 100 C 50 CONTINUE WRITE (6,FMT=9010) NEST,PTS ERROR = .TRUE. GO TO 100 C 60 CONTINUE R = 1 IF (Q.EQ.0) R = 0 IF (P.GE.1) THEN DO 70 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) 70 CONTINUE END IF IF (Q.GE.1) THEN DO 80 J = 1,Q JJ = P + J UNSYM(J) = QQ(J) FACTOR(JJ) = QQ(J) 80 CONTINUE SYM(P+1) = PTS/PSYM**2 END IF JJ = 2*P + Q FACTOR(JJ+1) = 0 PTWO = 1 J = 0 90 CONTINUE J = J + 1 IF (FACTOR(J).NE.0) THEN IF (FACTOR(J).EQ.2) THEN PTWO = PTWO*2 FACTOR(J) = 1 IF (PTWO.LT.TWOGRP) THEN IF (FACTOR(J+1).EQ.2) GO TO 90 END IF FACTOR(J) = PTWO PTWO = 1 END IF GO TO 90 END IF IF (P.EQ.0) R = 0 JJ = 2*P + R SYM(JJ+1) = 0 IF (Q.LE.1) Q = 0 UNSYM(Q+1) = 0 ERROR = .FALSE. C 100 CONTINUE C C---- Format statements C 9000 FORMAT (' FFTLIB: Largest factor exceeds ',I3,'. N = ',I6,'.') 9010 FORMAT (' FFTLIB: Factor count exceeds ',I3,'. N = ',I6,'.') END C C SUBROUTINE MDFTKD(N,FACTOR,DIM,X,Y) C ========================================== C C C---- Multi-dimensional complex fourier transform kernel driver C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL X(*),Y(*) INTEGER DIM(5),FACTOR(15) C .. C .. Local Scalars .. INTEGER F,M,P,R,S C .. C .. External Subroutines .. EXTERNAL R2CFTK,R3CFTK,R4CFTK,R5CFTK,R8CFTK,RPCFTK C .. S = DIM(2) F = 0 M = N 10 CONTINUE F = F + 1 P = FACTOR(F) IF (P.EQ.0) THEN RETURN ELSE M = M/P R = M*S IF (P.LE.8) THEN GO TO (10,20,30,40,50,80,70,60) P GO TO 80 C 20 CONTINUE CALL R2CFTK(N,M,X(1),Y(1),X(R+1),Y(R+1),DIM) GO TO 10 C 30 CONTINUE CALL R3CFTK(N,M,X(1),Y(1),X(R+1),Y(R+1),X(2*R+1),Y(2*R+1), + DIM) GO TO 10 C 40 CONTINUE CALL R4CFTK(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 10 C 50 CONTINUE CALL R5CFTK(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 10 C 60 CONTINUE CALL R8CFTK(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 10 END IF C 70 CONTINUE CALL RPCFTK(N,M,P,R,X,Y,DIM) GO TO 10 END IF C 80 CONTINUE WRITE (6,FMT=9000) C C---- Format statements C 9000 FORMAT ('0TRANSFER ERROR DETECTED IN MDFTKD',//) END C C SUBROUTINE R2CFTK(N,M,X0,Y0,X1,Y1,DIM) C ============================================== C C C---- Radix 2 multi-dimensional complex fourier transform kernel C C .. Scalar Arguments .. INTEGER M,N C .. C .. Array Arguments .. REAL X0(*),X1(*),Y0(*),Y1(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL ANGLE,C,FM2,IS,IU,RS,RU,S,TWOPI INTEGER J,K,K0,K1,K2,KK,L,L1,M2,MM2,MOVER2,NS,NT,SEP, + SIZE LOGICAL FOLD,ZERO C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. C .. Data statements .. DATA TWOPI/6.283185/ C .. C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M2 = M*2 FM2 = REAL(M2) MOVER2 = M/2 + 1 MM2 = SEP*M2 C DO 50 J = 1,MOVER2 FOLD = J .GT. 1 .AND. 2*J .LT. M + 2 K0 = (J-1)*SEP + 1 ZERO = j. eq. 1 IF (.NOT.ZERO) THEN ANGLE = TWOPI*REAL(j-1)/FM2 C = COS(ANGLE) S = SIN(ANGLE) END IF 10 CONTINUE C DO 40 KK = K0,NS,MM2 DO 30 L = KK,NT,L1 K1 = L + SIZE DO 20 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) THEN X1(K) = RU Y1(K) = IU ELSE X1(K) = RU*C + IU*S Y1(K) = IU*C - RU*S END IF 20 CONTINUE 30 CONTINUE 40 CONTINUE IF (FOLD) THEN FOLD = .FALSE. K0 = (M+1-J)*SEP + 1 C = -C GO TO 10 END IF 50 CONTINUE C END C C SUBROUTINE R3CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,DIM) C ====================================================== C C C---- Radix 3 multi-dimensional complex fourier transform kernel C C .. Scalar Arguments .. INTEGER M,N C .. C .. Array Arguments .. REAL X0(*),X1(*),X2(*),Y0(*),Y1(*),Y2(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL A,ANGLE,B,C1,C2,FM3,I0,I1,I2,IA,IB,IS,R0, + R1,R2,RA,RB,RS,S1,S2,T,TWOPI INTEGER J,K,K0,K1,K2,KK,L,L1,M3,MM3,MOVER2,NS,NT,SEP, + SIZE LOGICAL FOLD,ZERO C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. C .. Data statements .. DATA TWOPI/6.283185/,A/-0.5/,B/0.86602540/ C .. C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M3 = M*3 FM3 = REAL(M3) MM3 = SEP*M3 MOVER2 = M/2 + 1 C DO 50 J = 1,MOVER2 FOLD = J .GT. 1 .AND. 2*J .LT. M + 2 K0 = (J-1)*SEP + 1 ZERO = j .eq. 1 IF (.NOT.ZERO) THEN ANGLE = TWOPI*REAL(j-1)/FM3 C1 = COS(ANGLE) S1 = SIN(ANGLE) C2 = C1*C1 - S1*S1 S2 = S1*C1 + C1*S1 END IF 10 CONTINUE C DO 40 KK = K0,NS,MM3 DO 30 L = KK,NT,L1 K1 = L + SIZE DO 20 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 = RS*A + R0 IA = IS*A + I0 RB = (X1(K)-X2(K))*B IB = (Y1(K)-Y2(K))*B IF (ZERO) THEN X1(K) = RA + IB Y1(K) = IA - RB X2(K) = RA - IB Y2(K) = IA + RB ELSE 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 END IF 20 CONTINUE 30 CONTINUE 40 CONTINUE IF (FOLD) THEN 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 GO TO 10 END IF 50 CONTINUE C END C C SUBROUTINE R4CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,DIM) C ============================================================== C C C---- Radix 4 multi-dimensional complex fourier transform kernel C C .. Scalar Arguments .. INTEGER M,N C .. C .. Array Arguments .. REAL X0(*),X1(*),X2(*),X3(*),Y0(*),Y1(*), + Y2(*),Y3(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL ANGLE,C1,C2,C3,FM4,I1,I2,I3,IS0,IS1,IU0, + IU1,R1,R2,R3,RS0,RS1,RU0,RU1,S1,S2,S3,T,TWOPI INTEGER J,K,K0,K1,K2,KK,L,L1,M4,MM4,MOVER2,NS,NT,SEP, + SIZE LOGICAL FOLD,ZERO C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. C .. Data statements .. DATA TWOPI/6.283185/ C .. C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M4 = M*4 FM4 = REAL(M4) MM4 = SEP*M4 MOVER2 = M/2 + 1 C DO 50 J = 1,MOVER2 FOLD = J .GT. 1 .AND. 2*J .LT. M + 2 K0 = (J-1)*SEP + 1 ZERO = j .eq. 1 IF (.NOT.ZERO) THEN ANGLE = TWOPI*REAL(j-1)/FM4 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 END IF 10 CONTINUE C DO 40 KK = K0,NS,MM4 DO 30 L = KK,NT,L1 K1 = L + SIZE DO 20 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) THEN X2(K) = RU0 + IU1 Y2(K) = IU0 - RU1 X1(K) = RS0 - RS1 Y1(K) = IS0 - IS1 X3(K) = RU0 - IU1 Y3(K) = IU0 + RU1 ELSE 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 END IF 20 CONTINUE 30 CONTINUE 40 CONTINUE IF (FOLD) THEN FOLD = .FALSE. K0 = (M+1-J)*SEP + 1 T = C1 C1 = S1 S1 = T C2 = -C2 T = C3 C3 = -S3 S3 = -T GO TO 10 END IF 50 CONTINUE C END C C SUBROUTINE R5CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4,DIM) C ================================================================= C C C---- Radix 5 multi-dimensional complex fourier transform kernel C C .. Scalar Arguments .. INTEGER M,N C .. C .. Array Arguments .. REAL X0(*),X1(*),X2(*),X3(*),X4(*),Y0(*), + Y1(*),Y2(*),Y3(*),Y4(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL A1,A2,ANGLE,B1,B2,C1,C2,C3,C4,FM5,I0,I1,I2, + I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2,R0,R1,R2, + R3,R4,RA1,RA2,RB1,RB2,RS1,RS2,RU1,RU2,S1,S2,S3, + S4,T,TWOPI INTEGER J,K,K0,K1,K2,KK,L,L1,M5,MM5,MOVER2,NS,NT,SEP, + SIZE LOGICAL FOLD,ZERO C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. C .. Data statements .. DATA TWOPI/6.283185/,A1/0.30901699/,B1/0.95105652/, + A2/-0.80901699/,B2/0.58778525/ C .. C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M5 = M*5 FM5 = REAL(M5) MM5 = SEP*M5 MOVER2 = M/2 + 1 C DO 50 J = 1,MOVER2 FOLD = J .GT. 1 .AND. 2*J .LT. M + 2 K0 = (J-1)*SEP + 1 ZERO = j .eq. 1 IF (.NOT.ZERO) THEN ANGLE = TWOPI*REAL(j-1)/FM5 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 END IF 10 CONTINUE C DO 40 KK = K0,NS,MM5 DO 30 L = KK,NT,L1 K1 = L + SIZE DO 20 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 = RS1*A1 + R0 + RS2*A2 IA1 = IS1*A1 + I0 + IS2*A2 RA2 = RS1*A2 + R0 + RS2*A1 IA2 = IS1*A2 + I0 + IS2*A1 RB1 = RU1*B1 + RU2*B2 IB1 = IU1*B1 + IU2*B2 RB2 = RU1*B2 - RU2*B1 IB2 = IU1*B2 - IU2*B1 IF (ZERO) THEN 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 ELSE 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 END IF 20 CONTINUE 30 CONTINUE 40 CONTINUE IF (FOLD) THEN 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 GO TO 10 END IF 50 CONTINUE C END C C SUBROUTINE R8CFTK(N,M,X0,Y0,X1,Y1,X2,Y2,X3,Y3,X4,Y4,X5,Y5,X6,Y6, + X7,Y7,DIM) C =============================================================== C C C---- Radix 8 multi-dimensional complex fourier transform kernel C C .. Scalar Arguments .. INTEGER M,N C .. C .. Array Arguments .. REAL X0(*),X1(*),X2(*),X3(*),X4(*),X5(*), + X6(*),X7(*),Y0(*),Y1(*),Y2(*),Y3(*), + Y4(*),Y5(*),Y6(*),Y7(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL ANGLE,C1,C2,C3,C4,C5,C6,C7,E,FM8,I1,I2,I3, + I4,I5,I6,I7,IS0,IS1,IS2,IS3,ISS0,ISS1,ISU0,ISU1, + IU0,IU1,IU2,IU3,IUS0,IUS1,IUU0,IUU1,R1,R2,R3,R4, + R5,R6,R7,RS0,RS1,RS2,RS3,RSS0,RSS1,RSU0,RSU1, + RU0,RU1,RU2,RU3,RUS0,RUS1,RUU0,RUU1,S1,S2,S3,S4, + S5,S6,S7,T,TWOPI INTEGER J,K,K0,K1,K2,KK,L,L1,M8,MM8,MOVER2,NS,NT,SEP, + SIZE LOGICAL FOLD,ZERO C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. C .. Data statements .. DATA TWOPI/6.283185/,E/0.70710678/ C .. C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP M8 = M*8 FM8 = REAL(M8) MM8 = SEP*M8 MOVER2 = M/2 + 1 C DO 50 J = 1,MOVER2 FOLD = J .GT. 1 .AND. 2*J .LT. M + 2 K0 = (J-1)*SEP + 1 ZERO = j .eq. 1 IF (.NOT.ZERO) THEN ANGLE = TWOPI*REAL(j-1)/FM8 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 END IF 10 CONTINUE C DO 40 KK = K0,NS,MM8 DO 30 L = KK,NT,L1 K1 = L + SIZE DO 20 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) THEN 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 ELSE 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 END IF 20 CONTINUE 30 CONTINUE 40 CONTINUE IF (FOLD) THEN 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 GO TO 10 END IF 50 CONTINUE C END C C SUBROUTINE RPCFTK(N,M,P,R,X,Y,DIM) C ========================================== C C C---- Radix prime multi-dimensional complex fourier transform kernel C CMDW: Note, routine works beyond the nominal bounds of X and Y. C We think this is deliberate, so don't be tempted to "fix" it. C This routine is therefore incompatible with "bounds-checking" C options of compilers. C C .. Scalar Arguments .. INTEGER M,N,P,R C .. C .. Array Arguments .. REAL X(R,P),Y(R,P) INTEGER DIM(5) C .. C .. Local Scalars .. REAL ANGLE,FMP,FP,FU,IS,IU,RS,RU,T,TWOPI,XT,YT INTEGER J,JJ,K,K0,K1,K2,KK,L,L1,MMP,MOVER2,MP,NS,NT,PM, + PP,SEP,SIZE,U,V LOGICAL FOLD,ZERO C .. C .. Local Arrays .. REAL A(18),AA(9,9),B(18),BB(9,9),C(18),IA(9),IB(9), + RA(9),RB(9),S(18) C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. C .. Data statements .. DATA TWOPI/6.283185/ C .. C C NT = DIM(1) SEP = DIM(2) L1 = DIM(3) SIZE = DIM(4) - 1 K2 = DIM(5) NS = N*SEP MOVER2 = M/2 + 1 MP = M*P FMP = REAL(MP) MMP = SEP*MP PP = P/2 PM = P - 1 FP = REAL(P) FU = 0.0 DO 10 U = 1,PP FU = FU + 1.0 ANGLE = TWOPI*FU/FP JJ = P - U A(U) = COS(ANGLE) B(U) = SIN(ANGLE) A(JJ) = A(U) B(JJ) = -B(U) 10 CONTINUE DO 30 U = 1,PP DO 20 V = 1,PP JJ = U*V - U*V/P*P AA(V,U) = A(JJ) BB(V,U) = B(JJ) 20 CONTINUE 30 CONTINUE C DO 140 J = 1,MOVER2 FOLD = J .GT. 1 .AND. 2*J .LT. M + 2 K0 = (J-1)*SEP + 1 ZERO = j .eq. 1 IF (.NOT.ZERO) THEN ANGLE = TWOPI*REAL(j-1)/FMP C(1) = COS(ANGLE) S(1) = SIN(ANGLE) DO 40 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) 40 CONTINUE END IF 50 CONTINUE C DO 120 KK = K0,NS,MMP DO 110 L = KK,NT,L1 K1 = L + SIZE DO 100 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 60 U = 1,PP RA(U) = AA(U,1)*RS + XT IA(U) = AA(U,1)*IS + YT RB(U) = BB(U,1)*RU IB(U) = BB(U,1)*IU 60 CONTINUE XT = XT + RS YT = YT + IS DO 80 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 70 V = 1,PP RA(V) = AA(V,U)*RS + RA(V) IA(V) = AA(V,U)*IS + IA(V) RB(V) = BB(V,U)*RU + RB(V) IB(V) = BB(V,U)*IU + IB(V) 70 CONTINUE 80 CONTINUE X(K,1) = XT Y(K,1) = YT DO 90 U = 1,PP JJ = P - U IF (ZERO) THEN 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) ELSE XT = RA(U) + IB(U) YT = IA(U) - RB(U) X(K,U+1) = C(U)*XT + S(U)*YT Y(K,U+1) = C(U)*YT - S(U)*XT XT = RA(U) - IB(U) YT = IA(U) + RB(U) X(K,JJ+1) = C(JJ)*XT + S(JJ)*YT Y(K,JJ+1) = C(JJ)*YT - S(JJ)*XT END IF 90 CONTINUE 100 CONTINUE 110 CONTINUE 120 CONTINUE IF (FOLD) THEN FOLD = .FALSE. K0 = (M+1-J)*SEP + 1 DO 130 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 130 CONTINUE GO TO 50 END IF 140 CONTINUE C END C C SUBROUTINE HERMFT(X,Y,N,DIM) C ================================== C C C C---- Hermitian symmetric fourier transform C C Given the unique terms of a hermitian symmetric sequence of length C 2N this subroutine calculates the 2N real numbers which are its C fourier transform. The even numbered elements of the transform C (0, 2, 4, . . ., 2n-2) are returned in X and the odd numbered C elements (1, 3, 5, . . ., 2n-1) in Y. C C A finite hermitian sequence of length 2n contains n + 1 unique C real numbers and n - 1 unique imaginary numbers. For convenience C the real value for X(n) is stored at Y(0). C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL X(*),Y(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL A,ANGLE,B,C,CO,D,E,F,SI,TWON,TWOPI INTEGER D2,D3,D4,D5,I,I0,I1,I2,J,K,K1,NOVER2,NT C .. C .. External Subroutines .. EXTERNAL CMPLFT C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. TWOPI = 6.283185 TWON = REAL(2*N) C NT = DIM(1) D2 = DIM(2) D3 = DIM(3) D4 = DIM(4) - 1 D5 = DIM(5) C DO 20 I0 = 1,NT,D3 I1 = I0 + D4 DO 10 I = I0,I1,D5 A = X(I) B = Y(I) X(I) = A + B Y(I) = A - B 10 CONTINUE 20 CONTINUE C NOVER2 = N/2 + 1 IF (NOVER2.GE.2) THEN DO 50 I0 = 2,NOVER2 ANGLE = REAL(I0-1)*TWOPI/TWON CO = COS(ANGLE) SI = SIN(ANGLE) K = (N+2-2*I0)*D2 K1 = (I0-1)*D2 + 1 DO 40 I1 = K1,NT,D3 I2 = I1 + D4 DO 30 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 30 CONTINUE 40 CONTINUE 50 CONTINUE C CALL CMPLFT(X,Y,N,DIM) END IF C C END C C SUBROUTINE INV21(X,Y,N,D) C ========================= C C C C---- Inverts fourier transform along a screw C diad. the result is scaled by n. C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL X(*),Y(*) INTEGER D(5) C .. C .. Local Scalars .. REAL A,B,C,C1,PI,R,S,S1 INTEGER D1,D2,D3,D4,D5,I,J,J1,J2,J3,K,KK,L,LL,M,NOVER2 C .. C .. External Subroutines .. EXTERNAL CMPLFT C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. PI = 3.141593 C D1 = D(1) D2 = D(2) D3 = D(3) D4 = D(4) - 1 D5 = D(5) C NOVER2 = N/2 LL = N*D2 KK = NOVER2*D2 DO 20 J1 = 1,D1,D3 J2 = J1 + D4 DO 10 J = J1,J2,D5 L = LL + J K = KK + J X(L) = X(J) + X(K) X(K) = X(K) + Y(K) Y(L) = 0.0 Y(K) = 0.0 10 CONTINUE 20 CONTINUE C C1 = COS(PI/REAL(N)) S1 = SIN(PI/REAL(N)) C = 1.0 S = 0.0 DO 50 I = 2,NOVER2 KK = (N+2-2*I)*D2 LL = (N+1-I)*D2 R = C*C1 - S*S1 S = C*S1 + S*C1 C = R J1 = (I-1)*D2 + 1 DO 40 J2 = J1,D1,D3 J3 = J2 + D4 DO 30 J = J2,J3,D5 L = J + LL K = J + KK X(L) = X(L) + X(J) + X(K) X(J) = Y(J)*S + X(J) X(K) = Y(K)*S + X(K) Y(J) = Y(J)*C Y(K) = -Y(K)*C 30 CONTINUE 40 CONTINUE 50 CONTINUE C CALL CMPLFT(X,Y,N,D) C DO 80 I = 1,NOVER2 KK = (N+1-2*I)*D2 LL = I*D2 + KK J1 = (I-1)*D2 + 1 DO 70 J2 = J1,D1,D3 J3 = J2 + D4 DO 60 J = J2,J3,D5 K = J + KK L = J + LL A = X(J) - X(L) B = Y(J) + Y(L) X(J) = X(L) Y(J) = -Y(L) X(L) = X(K) + A Y(L) = Y(K) - B X(K) = A Y(K) = B 60 CONTINUE 70 CONTINUE 80 CONTINUE C M = N - 2 DO 130 I = 1,M K = I 90 CONTINUE J = K K = J/2 IF (2*K.NE.J) K = N - 1 - K IF (K-I) 90,130,100 100 KK = (K-I)*D2 J1 = I*D2 + 1 DO 120 J2 = J1,D1,D3 J3 = J2 + D4 DO 110 J = J2,J3,D5 K = J + KK A = X(K) B = Y(K) X(K) = X(J) Y(K) = Y(J) X(J) = A Y(J) = B 110 CONTINUE 120 CONTINUE 130 CONTINUE C C END C C SUBROUTINE REALFT(EVEN,ODD,N,DIM) C ====================================== C C C C REAL FOURIER TRANSFORM C C Given a real sequence of length 2n this subroutine calculates the C unique part of the fourier transform. The fourier transform has C n + 1 unique real parts and n - 1 unique imaginary parts. Since C the real part at x(n) is frequently of interest, this subroutine C stores it at x(n) rather than in y(0). Therefore x and y must be C of length n + 1 instead of n. Note that this storage arrangement C is different from that employed by the hermitian fourier transform C subroutine. C C For convenience the data is presented in two parts, the first C containing the even numbered real terms and the second containing C the odd numbered terms (numbering starting at 0). On return the C real part of the transform replaces the even terms and the C imaginary part of the transform replaces the odd terms. C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL EVEN(*),ODD(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL A,ANGLE,B,C,CO,D,E,F,SI,TWON,TWOPI INTEGER D2,D3,D4,D5,I,I0,I1,I2,J,K,L,NOVER2,NT C .. C .. External Subroutines .. EXTERNAL CMPLFT C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. TWOPI = 6.283185 TWON = REAL(2*N) C CALL CMPLFT(EVEN,ODD,N,DIM) C NT = DIM(1) D2 = DIM(2) D3 = DIM(3) D4 = DIM(4) - 1 D5 = DIM(5) NOVER2 = N/2 + 1 C IF (NOVER2.GE.2) THEN DO 30 I = 2,NOVER2 ANGLE = REAL(I-1)*TWOPI/TWON CO = COS(ANGLE) SI = SIN(ANGLE) I0 = (I-1)*D2 + 1 J = (N+2-2*I)*D2 DO 20 I1 = I0,NT,D3 I2 = I1 + D4 DO 10 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 10 CONTINUE 20 CONTINUE 30 CONTINUE END IF C IF (N.GE.1) THEN J = N*D2 DO 50 I1 = 1,NT,D3 I2 = I1 + D4 DO 40 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 40 CONTINUE 50 CONTINUE END IF C C END C C SUBROUTINE RSYMFT(X,N,DIM) C =============================== C C C C REAL SYMMETRIC MULTIDIMENSIONAL FOURIER TRANSFORM C C N must be a multiple of 4. The two unique elements are stored at C X(1) and X(n+1). C C Decimation in frequency applied to a real symmetric sequence of C length 2n gives a real symmetric sequence of length n, the C transform of which gives the even numbered fourier coefficients, C and a hermitian symmetric sequence of length n, the transform of C which gives the odd numbered fourier coefficients. The sum of C the two sequences is a hermitian symmetric sequence of length n, C which may be stored in n/2 complex locations. The transform of C this sequence is n real numbers representing the term by term sum C of the even and odd numbered fourier coefficients. This symmetric C sequence may be solved if any of the fourier coefficients are C known. For this purpose x0, which is simply the sum of the C original sequence, is computed and saved in x(n+1). C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL X(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL A,ANGLE,B,C,CO,D,SI,TWON,TWOPI INTEGER D1,D2,D3,D4,D5,I,I0,I1,I2,II,J,J0,J1,K,K0,K1, + K2,L,M,MJ,MK,ML,MM,NN,NOVER2,NOVER4, + TWOD2 CHARACTER EMESS*80 C .. C .. External Subroutines .. EXTERNAL HERMFT C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. IF (N.NE.1) THEN NOVER2 = N/2 NOVER4 = N/4 IF (4*NOVER4.NE.N) THEN C WRITE (EMESS,FMT=9000) N call ccperr(1, EMESS) ELSE D1 = DIM(1) D2 = DIM(2) D3 = DIM(3) D4 = DIM(4) - 1 D5 = DIM(5) TWOPI = 6.283185 TWON = REAL(2*N) TWOD2 = 2*D2 C K0 = N*D2 + 1 DO 20 K1 = K0,D1,D3 K2 = K1 + D4 DO 10 K = K1,K2,D5 X(K) = X(K)/2.0 10 CONTINUE 20 CONTINUE C DO 50 I = 2,NOVER2 ANGLE = REAL(I-1)*TWOPI/TWON CO = COS(ANGLE) SI = SIN(ANGLE) K0 = (I-1)*D2 + 1 J0 = (N+2-2*I)*D2 J1 = (N+1-I)*D2 DO 40 K1 = K0,D1,D3 K2 = K1 + D4 DO 30 K = K1,K2,D5 L = K + J0 NN = K + J1 A = X(L) + X(K) B = X(L) - X(K) X(K) = A - B*CO X(L) = B*SI X(NN) = X(NN) + A 30 CONTINUE 40 CONTINUE 50 CONTINUE C IF (NOVER4.NE.1) THEN J0 = NOVER4 - 1 DO 80 I = 1,J0 K0 = (NOVER2+I)*D2 + 1 J1 = (NOVER2-2*I)*D2 DO 70 K1 = K0,D1,D3 K2 = K1 + D4 DO 60 K = K1,K2,D5 L = K + J1 A = X(K) X(K) = X(L) X(L) = A 60 CONTINUE 70 CONTINUE 80 CONTINUE END IF C J0 = NOVER2*D2 J1 = N*D2 DO 100 K1 = 1,D1,D3 K2 = K1 + D4 DO 90 K = K1,K2,D5 I = K + J0 L = K + J1 X(I) = X(I)*2.0 X(L) = X(K) + X(I) + X(L)*2.0 X(K) = X(K)*2.0 90 CONTINUE 100 CONTINUE C K = NOVER2*D2 + 1 CALL HERMFT(X(1),X(K),NOVER2,DIM) C C---- Solve the equations for all of the sequences C I0 = 1 - D2 MK = NOVER2*D2 MJ = MK + D2 ML = N*D2 + D2 MM = ML DO 130 II = 1,NOVER4 I0 = I0 + D2 MJ = MJ - TWOD2 ML = ML - TWOD2 MM = MM - D2 DO 120 I1 = I0,D1,D3 I2 = I1 + D4 DO 110 I = I1,I2,D5 J = I + MJ K = I + MK L = I + ML M = I + MM A = X(I) - X(M) B = X(L) - A C = X(K) - B D = X(J) - C X(I) = X(M) X(J) = A X(K) = B X(L) = C X(M) = D 110 CONTINUE 120 CONTINUE 130 CONTINUE C C---- The results are now in a scrambled digit reversed order, i.e. C x(1), x(5), x(9), ..., x(10), x(6), x(2), ..., x(3), x(7), x(11), C ..., x(12), x(8), x(4). the following section of program follows C the permutation cycles and does the necessary interchanges. C IF (NOVER4.NE.1) THEN NN = N - 2 DO 170 I = 1,NN K = I 140 CONTINUE C K0 = K/4 L = K - K0*4 IF (L.NE. (L/2)*2) K0 = NOVER4 - 1 - K0 K = L*NOVER4 + K0 IF (K.LT.I) GO TO 140 IF (K.NE.I) THEN C K0 = I*D2 + 1 J0 = (K-I)*D2 DO 160 K1 = K0,D1,D3 K2 = K1 + D4 DO 150 K = K1,K2,D5 L = K + J0 A = X(K) X(K) = X(L) X(L) = A 150 CONTINUE 160 CONTINUE END IF 170 CONTINUE END IF END IF END IF C C---- Format statements C 9000 FORMAT ('FFTLIB: N not a multiple of 4 in R SYM FT. N =',I10,//) END C C SUBROUTINE SDIAD(X,Y,N,DIM) C =============================== C C C C This subroutine computes half the fourier synthesis along a screw C diad lying along a crystallographic axis given half the fourier C coefficients. That is, it assumes that f(t) = conjg(f(-t)) for t C even and f(t) = -conjg(f(-t)) for t odd. n is the length of the C desired half of the transform. The location x(n+1) is required as C a scratch location and therefore a value is also returned in C x(n+1) and y(n+1). The value of the second half of the transform C may be generated from the first half by the formula x(n+t) = x(t), C y(n+t) = -y(t). In other words, the last half of the transform is C the complex conjugate of the first half. C C The transform is calculated by forming the sum of the even terms C and the odd terms in place, using the symmetry relations to C obtain the values for negative subscripts. The transform of the C resulting sequence may be separated by using the fact that the C transform of the even terms is real, while the prodct of the C transform of the odd terms and (cos(pi*t/n) - i*sin(pi*t/n)) is C imaginary. The scratch location is required because the formula C for separating the two transforms breaks down when t = n/2. C C C Corrections from A.D.MCLACHLAN 1980, put here sep 1985 C errors in original algorithm for the scratch location which C assumed f(n)=0 C C C C .. Scalar Arguments .. INTEGER N C .. C .. Array Arguments .. REAL X(*),Y(*) INTEGER DIM(5) C .. C .. Local Scalars .. REAL A,ANGLE,C,S,TWON,TWOPI INTEGER D1,D2,D3,D4,D5,I,J,K,K0,K1,K2,L,M,MN,NN,NOVER2 LOGICAL FOLD CHARACTER EMESS*80 C .. C .. External Subroutines .. EXTERNAL CMPLFT C .. C .. Intrinsic Functions .. INTRINSIC COS,REAL,SIN C .. NOVER2 = N/2 IF (2*NOVER2.NE.N) THEN C WRITE (EMESS,FMT=9000) N call ccperr(1, EMESS) ELSE TWON = REAL(2*N) TWOPI = 6.2831852 D1 = DIM(1) D2 = DIM(2) D3 = DIM(3) D4 = DIM(4) - 1 D5 = DIM(5) C S = -1.0 IF (NOVER2.EQ. (2* (NOVER2/2))) S = -S K0 = (N-1)*D2 + 1 DO 20 K1 = K0,D1,D3 K2 = K1 + D4 DO 10 K = K1,K2,D5 L = K + D2 Y(L) = X(K)*S 10 CONTINUE 20 CONTINUE S = 1.0 NN = N - 2 DO 50 I = 1,NN,2 S = -S MN = (N+1-I)*D2 K0 = (I-1)*D2 + 1 DO 40 K1 = K0,D1,D3 K2 = K1 + D4 DO 30 K = K1,K2,D5 J = K + D2 L = 2*D2 + K M = K + MN Y(M) = X(J)*S + Y(M) X(K) = X(K) + X(J) X(J) = X(L) - X(J) Y(K) = Y(K) + Y(J) Y(J) = Y(J) - Y(L) 30 CONTINUE 40 CONTINUE 50 CONTINUE K0 = (N-2)*D2 + 1 DO 70 K1 = K0,D1,D3 K2 = K1 + D4 DO 60 K = K1,K2,D5 L = K + D2 X(K) = X(K) + X(L) Y(K) = Y(K) + Y(L) J = L + D2 X(L) = X(J) - X(L) 60 CONTINUE 70 CONTINUE C C---- Reorder scrambled fourier coefficients C DO 110 I = 1,NN K = I 80 CONTINUE K = 2*K IF (K.GT.N-1) K = 2*N - 1 - K IF (K.LT.I) GO TO 80 IF (K.NE.I) THEN J = (K-I)*D2 K0 = I*D2 + 1 DO 100 K1 = K0,D1,D3 K2 = K1 + D4 DO 90 K = K1,K2,D5 L = K + J A = X(K) X(K) = X(L) X(L) = A A = Y(K) Y(K) = Y(L) Y(L) = A 90 CONTINUE 100 CONTINUE END IF 110 CONTINUE C CALL CMPLFT(X,Y,N,DIM) C M = NOVER2 - 1 DO 150 I = 1,M ANGLE = REAL(I)*TWOPI/TWON C = COS(ANGLE) S = SIN(ANGLE) K0 = I*D2 + 1 FOLD = .TRUE. 120 CONTINUE C DO 140 K1 = K0,D1,D3 K2 = K1 + D4 DO 130 K = K1,K2,D5 A = Y(K)/C X(K) = X(K) + S*A Y(K) = A 130 CONTINUE 140 CONTINUE IF (FOLD) THEN C C = -C K0 = (N-I)*D2 + 1 FOLD = .FALSE. GO TO 120 END IF 150 CONTINUE C M = NOVER2*D2 K0 = M + 1 DO 170 K1 = K0,D1,D3 K2 = K1 + D4 DO 160 K = K1,K2,D5 J = K - M L = K + M A = Y(L)*2.0 X(K) = X(K) + A Y(K) = A X(L) = X(J) Y(L) = -Y(J) 160 CONTINUE 170 CONTINUE C END IF C C---- Format statements C 9000 FORMAT ('FFT error: SDIAD: N odd. N =',I10) END C C SUBROUTINE DIPRP(PTS,SYM,PSYM,UNSYM,DIM,X,Y) C ================================================= C C C---- Double in place reordering programme C C C .. Scalar Arguments .. INTEGER PSYM,PTS C .. C .. Array Arguments .. REAL X(*),Y(*) INTEGER DIM(5),SYM(15),UNSYM(15) C .. C .. Local Scalars .. REAL T INTEGER A,AL,B,BL,BS,C,CL,CS,D,DELTA,DK,DL,DS,E,EL,ES,F, + FL,FS,G,GL,GS,H,HL,HS,I,IL,IS,J,JJ,JL,JS,K,KK,KL, + KS,L,LK,LL,LS,M,ML,MODS,MS,MULT,N,NEST,NL,NS,NT, + P,P0,P1,P2,P3,P4,P5,PUNSYM,SEP,SIZE,TEST LOGICAL ONEMOD C .. C .. Local Arrays .. INTEGER MODULO(14),S(14),U(14) C .. C .. Equivalences .. 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)) C .. NEST = 14 C NT = DIM(1) SEP = DIM(2) P2 = DIM(3) SIZE = DIM(4) - 1 P4 = DIM(5) IF (SYM(1).NE.0) THEN DO 10 J = 1,NEST U(J) = 1 S(J) = 1 10 CONTINUE N = PTS DO 20 J = 1,NEST IF (SYM(J).EQ.0) THEN GO TO 30 ELSE JJ = NEST + 1 - J U(JJ) = N S(JJ) = N/SYM(J) N = N/SYM(J) END IF 20 CONTINUE C 30 JJ = 0 DO 190 A = 1,AL DO 180 B = A,BL,BS DO 170 C = B,CL,CS DO 160 D = C,DL,DS DO 150 E = D,EL,ES DO 140 F = E,FL,FS DO 130 G = F,GL,GS DO 120 H = G,HL,HS DO 110 I = H,IL,IS DO 100 J = I,JL,JS DO 90 K = J,KL,KS DO 80 L = K,LL,LS DO 70 M = L,ML,MS DO 60 N = M,NL,NS JJ = JJ + 1 IF (JJ.LT.N) THEN DELTA = (N-JJ)*SEP P1 = (JJ-1)*SEP + 1 DO 50 P0 = P1,NT,P2 P3 = P0 + SIZE DO 40 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 40 CONTINUE 50 CONTINUE END IF 60 CONTINUE 70 CONTINUE 80 CONTINUE 90 CONTINUE 100 CONTINUE 110 CONTINUE 120 CONTINUE 130 CONTINUE 140 CONTINUE 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE END IF C IF (UNSYM(1).NE.0) THEN PUNSYM = PTS/PSYM**2 MULT = PUNSYM/UNSYM(1) TEST = (UNSYM(1)*UNSYM(2)-1)*MULT*PSYM LK = MULT DK = MULT DO 200 K = 2,NEST IF (UNSYM(K).EQ.0) THEN GO TO 210 ELSE LK = UNSYM(K-1)*LK DK = DK/UNSYM(K) U(K) = (LK-DK)*PSYM MODS = K END IF 200 CONTINUE 210 ONEMOD = MODS .LT. 3 IF (.NOT.ONEMOD) THEN DO 220 J = 3,MODS JJ = MODS + 3 - J MODULO(JJ) = U(J) 220 CONTINUE END IF MODULO(2) = U(2) JL = (PUNSYM-3)*PSYM MS = PUNSYM*PSYM C DO 290 J = PSYM,JL,PSYM K = J 230 CONTINUE C K = K*MULT IF (.NOT.ONEMOD) THEN DO 240 I = 3,MODS K = K - (K/MODULO(I))*MODULO(I) 240 CONTINUE END IF IF (K.GE.TEST) THEN K = K - (K/MODULO(2))*MODULO(2) + MODULO(2) ELSE K = K - (K/MODULO(2))*MODULO(2) END IF IF (K.LT.J) GO TO 230 C IF (K.NE.J) THEN DELTA = (K-J)*SEP DO 280 L = 1,PSYM DO 270 M = L,PTS,MS P1 = (M+J-1)*SEP + 1 DO 260 P0 = P1,NT,P2 P3 = P0 + SIZE DO 250 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 250 CONTINUE 260 CONTINUE 270 CONTINUE 280 CONTINUE END IF 290 CONTINUE END IF C END C C Obscure routines only used by SFALL C =============================== SUBROUTINE PRMVCI(PERM,JV,N,N1) C =============================== C C---- Permute vector JV(N,3) by permutation matrix PERM C N1 is first dimension of JV C C .. Scalar Arguments .. INTEGER N,N1 C .. C .. Array Arguments .. REAL PERM(4,4) INTEGER JV(N1,3) C .. C .. Local Scalars .. INTEGER I C .. C .. Local Arrays .. REAL BV(3) C .. C .. Intrinsic Functions .. INTRINSIC NINT C .. C C---- Permute C DO 10 I = 1,3 BV(I) = PERM(I,1)*FLOAT(JV(N,1)) + PERM(I,2)*FLOAT(JV(N,2)) + + PERM(I,3)*FLOAT(JV(N,3)) 10 CONTINUE C C---- Copy back C DO 20 I = 1,3 JV(N,I) = NINT(BV(I)) 20 CONTINUE C END C C =============================== SUBROUTINE PRMVCR(PERM,AV,N,N1) C =============================== C C---- Permute vector AV(N,3) by permutation vector KP C N1 is first dimension of AV C C .. Scalar Arguments .. INTEGER N,N1 C .. C .. Array Arguments .. REAL AV(N1,3),PERM(4,4) C .. C .. Local Scalars .. INTEGER I C .. C .. Local Arrays .. REAL BV(3) C .. C C---- Permute C DO 10 I = 1,3 BV(I) = PERM(I,1)*AV(N,1) + PERM(I,2)*AV(N,2) + + PERM(I,3)*AV(N,3) 10 CONTINUE C C---- Copy back C DO 20 I = 1,3 AV(N,I) = BV(I) 20 CONTINUE C END C