SUBROUTINE DCFFTB(N,C,WSAVE) C***BEGIN PROLOGUE DCFFTB C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 860115 (YYMMDD) C***CATEGORY NO. J1A2 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE Unnormalized inverse of DCFFTF. C***DESCRIPTION C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1988 C C Subroutine DCFFTB computes the backward complex discrete Fourier C transform (the Fourier synthesis). Equivalently, DCFFTB computes C a complex periodic sequence from its Fourier coefficients. C The transform is defined below at output parameter C. C C A call of DCFFTF followed by a call of DCFFTB will multiply the C sequence by N. C C The array WSAVE which is used by subroutine DCFFTB must be C initialized by calling subroutine DCFFTI(N,WSAVE). C C Input Parameters C C C N the length of the complex sequence C. The method is C more efficient when N is the product of small primes. C C C a complex array of length N which contains the sequence C C WSAVE a d.p. work array which must be dimensioned at least 4*N+15 C in the program that calls DCFFTB. The WSAVE array must be C initialized by calling subroutine DCFFTI(N,WSAVE), and a C different WSAVE array must be used for each different C value of N. This initialization does not have to be C repeated so long as N remains unchanged. Thus subsequent C transforms can be obtained faster than the first. C The same WSAVE array can be used by DCFFTF and DCFFTB. C C Output Parameters C C C For J=1,...,N C C C(J)=the sum from K=1,...,N of C C C(K)*EXP(I*J*K*2*PI/N) C C where I=SQRT(-1) C C WSAVE contains initialization calculations which must not be C destroyed between calls of subroutine DCFFTF or DCFFTB C C * References * C * * C * 1. P.N. Swarztrauber, Vectorizing the FFTs, in Parallel * C * Computations (G. Rodrigue, ed.), Academic Press, 1982, * C * pp. 51-83. * C * 2. B.L. Buzbee, The SLATEC Common Math Library, in Sources * C * and Development of Mathematical Software (W. Cowell, ed.), * C * Prentice-Hall, 1984, pp. 302-318. * C * * C ********************************************************************* C C***REFERENCES (NONE) C***ROUTINES CALLED DCFTB1 C***END PROLOGUE DCFFTB IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION C(*) ,WSAVE(*) C***FIRST EXECUTABLE STATEMENT DCFFTB IF (N .EQ. 1) RETURN IW1 = N+N+1 IW2 = IW1+N+N CALL DCFTB1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2)) RETURN END SUBROUTINE DCFTB1(N,C,CH,WA,IFAC) C***BEGIN PROLOGUE DCFTB1 C***REFER TO DCFFTB C***ROUTINES CALLED DPASSB,DPASB2,DPASB3,DPASB4,DPASB5 C***END PROLOGUE DCFTB1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) C***FIRST EXECUTABLE STATEMENT DCFTB1 NF = IFAC(2) NA = 0 L1 = 1 IW = 1 DO 116 K1=1,NF IP = IFAC(K1+2) L2 = IP*L1 IDO = N/L2 IDOT = IDO+IDO IDL1 = IDOT*L1 IF (IP .NE. 4) GO TO 103 IX2 = IW+IDOT IX3 = IX2+IDOT IF (NA .NE. 0) GO TO 101 CALL DPASB4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) GO TO 102 101 CALL DPASB4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) 102 NA = 1-NA GO TO 115 103 IF (IP .NE. 2) GO TO 106 IF (NA .NE. 0) GO TO 104 CALL DPASB2 (IDOT,L1,C,CH,WA(IW)) GO TO 105 104 CALL DPASB2 (IDOT,L1,CH,C,WA(IW)) 105 NA = 1-NA GO TO 115 106 IF (IP .NE. 3) GO TO 109 IX2 = IW+IDOT IF (NA .NE. 0) GO TO 107 CALL DPASB3 (IDOT,L1,C,CH,WA(IW),WA(IX2)) GO TO 108 107 CALL DPASB3 (IDOT,L1,CH,C,WA(IW),WA(IX2)) 108 NA = 1-NA GO TO 115 109 IF (IP .NE. 5) GO TO 112 IX2 = IW+IDOT IX3 = IX2+IDOT IX4 = IX3+IDOT IF (NA .NE. 0) GO TO 110 CALL DPASB5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 111 110 CALL DPASB5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 111 NA = 1-NA GO TO 115 112 IF (NA .NE. 0) GO TO 113 CALL DPASSB (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) GO TO 114 113 CALL DPASSB (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) 114 IF (NAC .NE. 0) NA = 1-NA 115 L1 = L2 IW = IW+(IP-1)*IDOT 116 CONTINUE IF (NA .EQ. 0) RETURN N2 = N+N DO 117 I=1,N2 C(I) = CH(I) 117 CONTINUE RETURN END SUBROUTINE DPASSB(NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) C***BEGIN PROLOGUE DPASSB C***REFER TO DCFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASSB IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,WA(*) ,C2(IDL1,IP), 2 CH2(IDL1,IP) C***FIRST EXECUTABLE STATEMENT DPASSB IDOT = IDO/2 IPP2 = IP+2 IPPH = (IP+1)/2 IDP = IP*IDO C IF (IDO .LT. L1) GO TO 106 DO 103 J=2,IPPH JC = IPP2-J DO 102 K=1,L1 CDIR$ IVDEP DO 101 I=1,IDO CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) 101 CONTINUE 102 CONTINUE 103 CONTINUE DO 105 K=1,L1 CDIR$ IVDEP DO 104 I=1,IDO CH(I,K,1) = CC(I,1,K) 104 CONTINUE 105 CONTINUE GO TO 112 106 DO 109 J=2,IPPH JC = IPP2-J DO 108 I=1,IDO CDIR$ IVDEP DO 107 K=1,L1 CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) 107 CONTINUE 108 CONTINUE 109 CONTINUE DO 111 I=1,IDO CDIR$ IVDEP DO 110 K=1,L1 CH(I,K,1) = CC(I,1,K) 110 CONTINUE 111 CONTINUE 112 IDL = 2-IDO INC = 0 DO 116 L=2,IPPH LC = IPP2-L IDL = IDL+IDO CDIR$ IVDEP DO 113 IK=1,IDL1 C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2) C2(IK,LC) = WA(IDL)*CH2(IK,IP) 113 CONTINUE IDLJ = IDL INC = INC+IDO DO 115 J=3,IPPH JC = IPP2-J IDLJ = IDLJ+INC IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP WAR = WA(IDLJ-1) WAI = WA(IDLJ) CDIR$ IVDEP DO 114 IK=1,IDL1 C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J) C2(IK,LC) = C2(IK,LC)+WAI*CH2(IK,JC) 114 CONTINUE 115 CONTINUE 116 CONTINUE DO 118 J=2,IPPH CDIR$ IVDEP DO 117 IK=1,IDL1 CH2(IK,1) = CH2(IK,1)+CH2(IK,J) 117 CONTINUE 118 CONTINUE DO 120 J=2,IPPH JC = IPP2-J CDIR$ IVDEP DO 119 IK=2,IDL1,2 CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC) CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC) CH2(IK,J) = C2(IK,J)+C2(IK-1,JC) CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC) 119 CONTINUE 120 CONTINUE NAC = 1 IF (IDO .EQ. 2) RETURN NAC = 0 DO 121 IK=1,IDL1 C2(IK,1) = CH2(IK,1) 121 CONTINUE DO 123 J=2,IP CDIR$ IVDEP DO 122 K=1,L1 C1(1,K,J) = CH(1,K,J) C1(2,K,J) = CH(2,K,J) 122 CONTINUE 123 CONTINUE IF (IDOT .GT. L1) GO TO 127 IDIJ = 0 DO 126 J=2,IP IDIJ = IDIJ+2 DO 125 I=4,IDO,2 IDIJ = IDIJ+2 CDIR$ IVDEP DO 124 K=1,L1 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) 124 CONTINUE 125 CONTINUE 126 CONTINUE RETURN 127 IDJ = 2-IDO DO 130 J=2,IP IDJ = IDJ+IDO DO 129 K=1,L1 IDIJ = IDJ CDIR$ IVDEP DO 128 I=4,IDO,2 IDIJ = IDIJ+2 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) 128 CONTINUE 129 CONTINUE 130 CONTINUE RETURN END SUBROUTINE DPASB2(IDO,L1,CC,CH,WA1) C***BEGIN PROLOGUE DPASB2 C***REFER TO DCFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASB2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 1 WA1(*) C***FIRST EXECUTABLE STATEMENT DPASB2 IF (IDO .GT. 2) GO TO 102 DO 101 K=1,L1 CH(1,K,1) = CC(1,1,K)+CC(1,2,K) CH(1,K,2) = CC(1,1,K)-CC(1,2,K) CH(2,K,1) = CC(2,1,K)+CC(2,2,K) CH(2,K,2) = CC(2,1,K)-CC(2,2,K) 101 CONTINUE RETURN 102 IF(IDO/2.LT.L1) GO TO 105 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=2,IDO,2 CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) TR2 = CC(I-1,1,K)-CC(I-1,2,K) CH(I,K,1) = CC(I,1,K)+CC(I,2,K) TI2 = CC(I,1,K)-CC(I,2,K) CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2 CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2 103 CONTINUE 104 CONTINUE RETURN 105 DO 107 I=2,IDO,2 CDIR$ IVDEP DO 106 K=1,L1 CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) TR2 = CC(I-1,1,K)-CC(I-1,2,K) CH(I,K,1) = CC(I,1,K)+CC(I,2,K) TI2 = CC(I,1,K)-CC(I,2,K) CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2 CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2 106 CONTINUE 107 CONTINUE RETURN END SUBROUTINE DPASB3(IDO,L1,CC,CH,WA1,WA2) C***BEGIN PROLOGUE DPASB3 C***REFER TO DCFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASB3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 1 WA1(*) ,WA2(*) C***FIRST EXECUTABLE STATEMENT DPASB3 TAUR = -.5D0 TAUI = .5D0*SQRT(3.0D0) IF (IDO .NE. 2) GO TO 102 DO 101 K=1,L1 TR2 = CC(1,2,K)+CC(1,3,K) CR2 = CC(1,1,K)+TAUR*TR2 CH(1,K,1) = CC(1,1,K)+TR2 TI2 = CC(2,2,K)+CC(2,3,K) CI2 = CC(2,1,K)+TAUR*TI2 CH(2,K,1) = CC(2,1,K)+TI2 CR3 = TAUI*(CC(1,2,K)-CC(1,3,K)) CI3 = TAUI*(CC(2,2,K)-CC(2,3,K)) CH(1,K,2) = CR2-CI3 CH(1,K,3) = CR2+CI3 CH(2,K,2) = CI2+CR3 CH(2,K,3) = CI2-CR3 101 CONTINUE RETURN 102 IF(IDO/2.LT.L1) GO TO 105 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=2,IDO,2 TR2 = CC(I-1,2,K)+CC(I-1,3,K) CR2 = CC(I-1,1,K)+TAUR*TR2 CH(I-1,K,1) = CC(I-1,1,K)+TR2 TI2 = CC(I,2,K)+CC(I,3,K) CI2 = CC(I,1,K)+TAUR*TI2 CH(I,K,1) = CC(I,1,K)+TI2 CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) DR2 = CR2-CI3 DR3 = CR2+CI3 DI2 = CI2+CR3 DI3 = CI2-CR3 CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 103 CONTINUE 104 CONTINUE RETURN 105 DO 107 I=2,IDO,2 CDIR$ IVDEP DO 106 K=1,L1 TR2 = CC(I-1,2,K)+CC(I-1,3,K) CR2 = CC(I-1,1,K)+TAUR*TR2 CH(I-1,K,1) = CC(I-1,1,K)+TR2 TI2 = CC(I,2,K)+CC(I,3,K) CI2 = CC(I,1,K)+TAUR*TI2 CH(I,K,1) = CC(I,1,K)+TI2 CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) DR2 = CR2-CI3 DR3 = CR2+CI3 DI2 = CI2+CR3 DI3 = CI2-CR3 CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 106 CONTINUE 107 CONTINUE RETURN END SUBROUTINE DPASB4(IDO,L1,CC,CH,WA1,WA2,WA3) C***BEGIN PROLOGUE DPASB4 C***REFER TO DCFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASB4 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 1 WA1(*) ,WA2(*) ,WA3(*) C***FIRST EXECUTABLE STATEMENT DPASB4 IF (IDO .NE. 2) GO TO 102 DO 101 K=1,L1 TI1 = CC(2,1,K)-CC(2,3,K) TI2 = CC(2,1,K)+CC(2,3,K) TR4 = CC(2,4,K)-CC(2,2,K) TI3 = CC(2,2,K)+CC(2,4,K) TR1 = CC(1,1,K)-CC(1,3,K) TR2 = CC(1,1,K)+CC(1,3,K) TI4 = CC(1,2,K)-CC(1,4,K) TR3 = CC(1,2,K)+CC(1,4,K) CH(1,K,1) = TR2+TR3 CH(1,K,3) = TR2-TR3 CH(2,K,1) = TI2+TI3 CH(2,K,3) = TI2-TI3 CH(1,K,2) = TR1+TR4 CH(1,K,4) = TR1-TR4 CH(2,K,2) = TI1+TI4 CH(2,K,4) = TI1-TI4 101 CONTINUE RETURN 102 IF(IDO/2.LT.L1) GO TO 105 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=2,IDO,2 TI1 = CC(I,1,K)-CC(I,3,K) TI2 = CC(I,1,K)+CC(I,3,K) TI3 = CC(I,2,K)+CC(I,4,K) TR4 = CC(I,4,K)-CC(I,2,K) TR1 = CC(I-1,1,K)-CC(I-1,3,K) TR2 = CC(I-1,1,K)+CC(I-1,3,K) TI4 = CC(I-1,2,K)-CC(I-1,4,K) TR3 = CC(I-1,2,K)+CC(I-1,4,K) CH(I-1,K,1) = TR2+TR3 CR3 = TR2-TR3 CH(I,K,1) = TI2+TI3 CI3 = TI2-TI3 CR2 = TR1+TR4 CR4 = TR1-TR4 CI2 = TI1+TI4 CI4 = TI1-TI4 CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2 CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2 CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3 CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3 CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4 CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4 103 CONTINUE 104 CONTINUE RETURN 105 DO 107 I=2,IDO,2 CDIR$ IVDEP DO 106 K=1,L1 TI1 = CC(I,1,K)-CC(I,3,K) TI2 = CC(I,1,K)+CC(I,3,K) TI3 = CC(I,2,K)+CC(I,4,K) TR4 = CC(I,4,K)-CC(I,2,K) TR1 = CC(I-1,1,K)-CC(I-1,3,K) TR2 = CC(I-1,1,K)+CC(I-1,3,K) TI4 = CC(I-1,2,K)-CC(I-1,4,K) TR3 = CC(I-1,2,K)+CC(I-1,4,K) CH(I-1,K,1) = TR2+TR3 CR3 = TR2-TR3 CH(I,K,1) = TI2+TI3 CI3 = TI2-TI3 CR2 = TR1+TR4 CR4 = TR1-TR4 CI2 = TI1+TI4 CI4 = TI1-TI4 CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2 CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2 CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3 CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3 CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4 CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4 106 CONTINUE 107 CONTINUE RETURN END SUBROUTINE DPASB5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4) C***BEGIN PROLOGUE DPASB5 C***REFER TO DCFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASB5 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*) C***FIRST EXECUTABLE STATEMENT DPASB5 PI = 4.0D0*ATAN(1.0D0) TR11 = SIN(.1D0*PI) TI11 = SIN(.4D0*PI) TR12 = -SIN(.3D0*PI) TI12 = SIN(.2D0*PI) IF (IDO .NE. 2) GO TO 102 DO 101 K=1,L1 TI5 = CC(2,2,K)-CC(2,5,K) TI2 = CC(2,2,K)+CC(2,5,K) TI4 = CC(2,3,K)-CC(2,4,K) TI3 = CC(2,3,K)+CC(2,4,K) TR5 = CC(1,2,K)-CC(1,5,K) TR2 = CC(1,2,K)+CC(1,5,K) TR4 = CC(1,3,K)-CC(1,4,K) TR3 = CC(1,3,K)+CC(1,4,K) CH(1,K,1) = CC(1,1,K)+TR2+TR3 CH(2,K,1) = CC(2,1,K)+TI2+TI3 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3 CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3 CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3 CR5 = TI11*TR5+TI12*TR4 CI5 = TI11*TI5+TI12*TI4 CR4 = TI12*TR5-TI11*TR4 CI4 = TI12*TI5-TI11*TI4 CH(1,K,2) = CR2-CI5 CH(1,K,5) = CR2+CI5 CH(2,K,2) = CI2+CR5 CH(2,K,3) = CI3+CR4 CH(1,K,3) = CR3-CI4 CH(1,K,4) = CR3+CI4 CH(2,K,4) = CI3-CR4 CH(2,K,5) = CI2-CR5 101 CONTINUE RETURN 102 IF(IDO/2.LT.L1) GO TO 105 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=2,IDO,2 TI5 = CC(I,2,K)-CC(I,5,K) TI2 = CC(I,2,K)+CC(I,5,K) TI4 = CC(I,3,K)-CC(I,4,K) TI3 = CC(I,3,K)+CC(I,4,K) TR5 = CC(I-1,2,K)-CC(I-1,5,K) TR2 = CC(I-1,2,K)+CC(I-1,5,K) TR4 = CC(I-1,3,K)-CC(I-1,4,K) TR3 = CC(I-1,3,K)+CC(I-1,4,K) CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 CH(I,K,1) = CC(I,1,K)+TI2+TI3 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 CR5 = TI11*TR5+TI12*TR4 CI5 = TI11*TI5+TI12*TI4 CR4 = TI12*TR5-TI11*TR4 CI4 = TI12*TI5-TI11*TI4 DR3 = CR3-CI4 DR4 = CR3+CI4 DI3 = CI3+CR4 DI4 = CI3-CR4 DR5 = CR2+CI5 DR2 = CR2-CI5 DI5 = CI2-CR5 DI2 = CI2+CR5 CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4 CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4 CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5 CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5 103 CONTINUE 104 CONTINUE RETURN 105 DO 107 I=2,IDO,2 CDIR$ IVDEP DO 106 K=1,L1 TI5 = CC(I,2,K)-CC(I,5,K) TI2 = CC(I,2,K)+CC(I,5,K) TI4 = CC(I,3,K)-CC(I,4,K) TI3 = CC(I,3,K)+CC(I,4,K) TR5 = CC(I-1,2,K)-CC(I-1,5,K) TR2 = CC(I-1,2,K)+CC(I-1,5,K) TR4 = CC(I-1,3,K)-CC(I-1,4,K) TR3 = CC(I-1,3,K)+CC(I-1,4,K) CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 CH(I,K,1) = CC(I,1,K)+TI2+TI3 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 CR5 = TI11*TR5+TI12*TR4 CI5 = TI11*TI5+TI12*TI4 CR4 = TI12*TR5-TI11*TR4 CI4 = TI12*TI5-TI11*TI4 DR3 = CR3-CI4 DR4 = CR3+CI4 DI3 = CI3+CR4 DI4 = CI3-CR4 DR5 = CR2+CI5 DR2 = CR2-CI5 DI5 = CI2-CR5 DI2 = CI2+CR5 CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2 CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2 CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3 CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3 CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4 CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4 CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5 CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5 106 CONTINUE 107 CONTINUE RETURN END SUBROUTINE DCFFTF(N,C,WSAVE) C***BEGIN PROLOGUE DCFFTF C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 860115 (YYMMDD) C***CATEGORY NO. J1A2 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE Forward transform of a complex, periodic sequence. C***DESCRIPTION C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1988 C C Subroutine DCFFTF computes the forward complex discrete Fourier C transform (the Fourier analysis). Equivalently, DCFFTF computes C the Fourier coefficients of a complex periodic sequence. C The transform is defined below at output parameter C. C C The transform is not normalized. To obtain a normalized transform C the output must be divided by N. Otherwise a call of DCFFTF C followed by a call of DCFFTB will multiply the sequence by N. C C The array WSAVE which is used by subroutine DCFFTF must be C initialized by calling subroutine DCFFTI(N,WSAVE). C C Input Parameters C C C N the length of the complex sequence C. The method is C more efficient when N is the product of small primes. C C C a complex array of length N which contains the sequence C C WSAVE a d.p. work array which must be dimensioned at least 4*N+15 C in the program that calls DCFFTF. The WSAVE array must be C initialized by calling subroutine DCFFTI(N,WSAVE), and a C different WSAVE array must be used for each different C value of N. This initialization does not have to be C repeated so long as N remains unchanged. Thus subsequent C transforms can be obtained faster than the first. C The same WSAVE array can be used by DCFFTF and DCFFTB. C C Output Parameters C C C for J=1,...,N C C C(J)=the sum from K=1,...,N of C C C(K)*EXP(-I*J*K*2*PI/N) C C where I=SQRT(-1) C C WSAVE contains initialization calculations which must not be C destroyed between calls of subroutine DCFFTF or DCFFTB C C * References * C * * C * 1. P.N. Swarztrauber, Vectorizing the FFTs, in Parallel * C * Computations (G. Rodrigue, ed.), Academic Press, 1982, * C * pp. 51-83. * C * 2. B.L. Buzbee, The SLATEC Common Math Library, in Sources * C * and Development of Mathematical Software (W. Cowell, ed.), * C * Prentice-Hall, 1984, pp. 302-318. * C * * C ********************************************************************* C C***REFERENCES (NONE) C***ROUTINES CALLED DCFTF1 C***END PROLOGUE DCFFTF IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION C(*) ,WSAVE(*) C***FIRST EXECUTABLE STATEMENT DCFFTF IF (N .EQ. 1) RETURN IW1 = N+N+1 IW2 = IW1+N+N CALL DCFTF1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2)) RETURN END SUBROUTINE DCFTF1(N,C,CH,WA,IFAC) C***BEGIN PROLOGUE DCFTF1 C***REFER TO DCFFTF C***ROUTINES CALLED DPASSF,DPASF2,DPASF3,DPASF4,DPASF5 C***END PROLOGUE DCFTF1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) C***FIRST EXECUTABLE STATEMENT DCFTF1 NF = IFAC(2) NA = 0 L1 = 1 IW = 1 DO 116 K1=1,NF IP = IFAC(K1+2) L2 = IP*L1 IDO = N/L2 IDOT = IDO+IDO IDL1 = IDOT*L1 IF (IP .NE. 4) GO TO 103 IX2 = IW+IDOT IX3 = IX2+IDOT IF (NA .NE. 0) GO TO 101 CALL DPASF4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) GO TO 102 101 CALL DPASF4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) 102 NA = 1-NA GO TO 115 103 IF (IP .NE. 2) GO TO 106 IF (NA .NE. 0) GO TO 104 CALL DPASF2 (IDOT,L1,C,CH,WA(IW)) GO TO 105 104 CALL DPASF2 (IDOT,L1,CH,C,WA(IW)) 105 NA = 1-NA GO TO 115 106 IF (IP .NE. 3) GO TO 109 IX2 = IW+IDOT IF (NA .NE. 0) GO TO 107 CALL DPASF3 (IDOT,L1,C,CH,WA(IW),WA(IX2)) GO TO 108 107 CALL DPASF3 (IDOT,L1,CH,C,WA(IW),WA(IX2)) 108 NA = 1-NA GO TO 115 109 IF (IP .NE. 5) GO TO 112 IX2 = IW+IDOT IX3 = IX2+IDOT IX4 = IX3+IDOT IF (NA .NE. 0) GO TO 110 CALL DPASF5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 111 110 CALL DPASF5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 111 NA = 1-NA GO TO 115 112 IF (NA .NE. 0) GO TO 113 CALL DPASSF (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) GO TO 114 113 CALL DPASSF (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) 114 IF (NAC .NE. 0) NA = 1-NA 115 L1 = L2 IW = IW+(IP-1)*IDOT 116 CONTINUE IF (NA .EQ. 0) RETURN N2 = N+N DO 117 I=1,N2 C(I) = CH(I) 117 CONTINUE RETURN END SUBROUTINE DPASSF(NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) C***BEGIN PROLOGUE DPASSF C***REFER TO DCFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASSF IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,WA(*) ,C2(IDL1,IP), 2 CH2(IDL1,IP) C***FIRST EXECUTABLE STATEMENT DPASSF IDOT = IDO/2 IPP2 = IP+2 IPPH = (IP+1)/2 IDP = IP*IDO C IF (IDO .LT. L1) GO TO 106 DO 103 J=2,IPPH JC = IPP2-J DO 102 K=1,L1 CDIR$ IVDEP DO 101 I=1,IDO CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) 101 CONTINUE 102 CONTINUE 103 CONTINUE DO 105 K=1,L1 CDIR$ IVDEP DO 104 I=1,IDO CH(I,K,1) = CC(I,1,K) 104 CONTINUE 105 CONTINUE GO TO 112 106 DO 109 J=2,IPPH JC = IPP2-J DO 108 I=1,IDO CDIR$ IVDEP DO 107 K=1,L1 CH(I,K,J) = CC(I,J,K)+CC(I,JC,K) CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K) 107 CONTINUE 108 CONTINUE 109 CONTINUE DO 111 I=1,IDO CDIR$ IVDEP DO 110 K=1,L1 CH(I,K,1) = CC(I,1,K) 110 CONTINUE 111 CONTINUE 112 IDL = 2-IDO INC = 0 DO 116 L=2,IPPH LC = IPP2-L IDL = IDL+IDO CDIR$ IVDEP DO 113 IK=1,IDL1 C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2) C2(IK,LC) = -WA(IDL)*CH2(IK,IP) 113 CONTINUE IDLJ = IDL INC = INC+IDO DO 115 J=3,IPPH JC = IPP2-J IDLJ = IDLJ+INC IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP WAR = WA(IDLJ-1) WAI = WA(IDLJ) CDIR$ IVDEP DO 114 IK=1,IDL1 C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J) C2(IK,LC) = C2(IK,LC)-WAI*CH2(IK,JC) 114 CONTINUE 115 CONTINUE 116 CONTINUE DO 118 J=2,IPPH CDIR$ IVDEP DO 117 IK=1,IDL1 CH2(IK,1) = CH2(IK,1)+CH2(IK,J) 117 CONTINUE 118 CONTINUE DO 120 J=2,IPPH JC = IPP2-J CDIR$ IVDEP DO 119 IK=2,IDL1,2 CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC) CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC) CH2(IK,J) = C2(IK,J)+C2(IK-1,JC) CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC) 119 CONTINUE 120 CONTINUE NAC = 1 IF (IDO .EQ. 2) RETURN NAC = 0 CDIR$ IVDEP DO 121 IK=1,IDL1 C2(IK,1) = CH2(IK,1) 121 CONTINUE DO 123 J=2,IP CDIR$ IVDEP DO 122 K=1,L1 C1(1,K,J) = CH(1,K,J) C1(2,K,J) = CH(2,K,J) 122 CONTINUE 123 CONTINUE IF (IDOT .GT. L1) GO TO 127 IDIJ = 0 DO 126 J=2,IP IDIJ = IDIJ+2 DO 125 I=4,IDO,2 IDIJ = IDIJ+2 CDIR$ IVDEP DO 124 K=1,L1 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J) C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J) 124 CONTINUE 125 CONTINUE 126 CONTINUE RETURN 127 IDJ = 2-IDO DO 130 J=2,IP IDJ = IDJ+IDO DO 129 K=1,L1 IDIJ = IDJ CDIR$ IVDEP DO 128 I=4,IDO,2 IDIJ = IDIJ+2 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J) C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J) 128 CONTINUE 129 CONTINUE 130 CONTINUE RETURN END SUBROUTINE DPASF2(IDO,L1,CC,CH,WA1) C***BEGIN PROLOGUE DPASF2 C***REFER TO DCFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASF2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 1 WA1(*) C***FIRST EXECUTABLE STATEMENT DPASF2 IF (IDO .GT. 2) GO TO 102 DO 101 K=1,L1 CH(1,K,1) = CC(1,1,K)+CC(1,2,K) CH(1,K,2) = CC(1,1,K)-CC(1,2,K) CH(2,K,1) = CC(2,1,K)+CC(2,2,K) CH(2,K,2) = CC(2,1,K)-CC(2,2,K) 101 CONTINUE RETURN 102 IF(IDO/2.LT.L1) GO TO 105 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=2,IDO,2 CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) TR2 = CC(I-1,1,K)-CC(I-1,2,K) CH(I,K,1) = CC(I,1,K)+CC(I,2,K) TI2 = CC(I,1,K)-CC(I,2,K) CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2 CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2 103 CONTINUE 104 CONTINUE RETURN 105 DO 107 I=2,IDO,2 CDIR$ IVDEP DO 106 K=1,L1 CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K) TR2 = CC(I-1,1,K)-CC(I-1,2,K) CH(I,K,1) = CC(I,1,K)+CC(I,2,K) TI2 = CC(I,1,K)-CC(I,2,K) CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2 CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2 106 CONTINUE 107 CONTINUE RETURN END SUBROUTINE DPASF3(IDO,L1,CC,CH,WA1,WA2) C***BEGIN PROLOGUE DPASF3 C***REFER TO DCFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASF3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 1 WA1(*) ,WA2(*) C***FIRST EXECUTABLE STATEMENT DPASF3 TAUR = -.5D0 TAUI = -.5D0*SQRT(3.0D0) IF (IDO .NE. 2) GO TO 102 DO 101 K=1,L1 TR2 = CC(1,2,K)+CC(1,3,K) CR2 = CC(1,1,K)+TAUR*TR2 CH(1,K,1) = CC(1,1,K)+TR2 TI2 = CC(2,2,K)+CC(2,3,K) CI2 = CC(2,1,K)+TAUR*TI2 CH(2,K,1) = CC(2,1,K)+TI2 CR3 = TAUI*(CC(1,2,K)-CC(1,3,K)) CI3 = TAUI*(CC(2,2,K)-CC(2,3,K)) CH(1,K,2) = CR2-CI3 CH(1,K,3) = CR2+CI3 CH(2,K,2) = CI2+CR3 CH(2,K,3) = CI2-CR3 101 CONTINUE RETURN 102 IF(IDO/2.LT.L1) GO TO 105 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=2,IDO,2 TR2 = CC(I-1,2,K)+CC(I-1,3,K) CR2 = CC(I-1,1,K)+TAUR*TR2 CH(I-1,K,1) = CC(I-1,1,K)+TR2 TI2 = CC(I,2,K)+CC(I,3,K) CI2 = CC(I,1,K)+TAUR*TI2 CH(I,K,1) = CC(I,1,K)+TI2 CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) DR2 = CR2-CI3 DR3 = CR2+CI3 DI2 = CI2+CR3 DI3 = CI2-CR3 CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 103 CONTINUE 104 CONTINUE RETURN 105 DO 107 I=2,IDO,2 CDIR$ IVDEP DO 106 K=1,L1 TR2 = CC(I-1,2,K)+CC(I-1,3,K) CR2 = CC(I-1,1,K)+TAUR*TR2 CH(I-1,K,1) = CC(I-1,1,K)+TR2 TI2 = CC(I,2,K)+CC(I,3,K) CI2 = CC(I,1,K)+TAUR*TI2 CH(I,K,1) = CC(I,1,K)+TI2 CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K)) CI3 = TAUI*(CC(I,2,K)-CC(I,3,K)) DR2 = CR2-CI3 DR3 = CR2+CI3 DI2 = CI2+CR3 DI3 = CI2-CR3 CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 106 CONTINUE 107 CONTINUE RETURN END SUBROUTINE DPASF4(IDO,L1,CC,CH,WA1,WA2,WA3) C***BEGIN PROLOGUE DPASF4 C***REFER TO DCFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASF4 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 1 WA1(*) ,WA2(*) ,WA3(*) C***FIRST EXECUTABLE STATEMENT DPASF4 IF (IDO .NE. 2) GO TO 102 DO 101 K=1,L1 TI1 = CC(2,1,K)-CC(2,3,K) TI2 = CC(2,1,K)+CC(2,3,K) TR4 = CC(2,2,K)-CC(2,4,K) TI3 = CC(2,2,K)+CC(2,4,K) TR1 = CC(1,1,K)-CC(1,3,K) TR2 = CC(1,1,K)+CC(1,3,K) TI4 = CC(1,4,K)-CC(1,2,K) TR3 = CC(1,2,K)+CC(1,4,K) CH(1,K,1) = TR2+TR3 CH(1,K,3) = TR2-TR3 CH(2,K,1) = TI2+TI3 CH(2,K,3) = TI2-TI3 CH(1,K,2) = TR1+TR4 CH(1,K,4) = TR1-TR4 CH(2,K,2) = TI1+TI4 CH(2,K,4) = TI1-TI4 101 CONTINUE RETURN 102 IF(IDO/2.LT.L1) GO TO 105 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=2,IDO,2 TI1 = CC(I,1,K)-CC(I,3,K) TI2 = CC(I,1,K)+CC(I,3,K) TI3 = CC(I,2,K)+CC(I,4,K) TR4 = CC(I,2,K)-CC(I,4,K) TR1 = CC(I-1,1,K)-CC(I-1,3,K) TR2 = CC(I-1,1,K)+CC(I-1,3,K) TI4 = CC(I-1,4,K)-CC(I-1,2,K) TR3 = CC(I-1,2,K)+CC(I-1,4,K) CH(I-1,K,1) = TR2+TR3 CR3 = TR2-TR3 CH(I,K,1) = TI2+TI3 CI3 = TI2-TI3 CR2 = TR1+TR4 CR4 = TR1-TR4 CI2 = TI1+TI4 CI4 = TI1-TI4 CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2 CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2 CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3 CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3 CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4 CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4 103 CONTINUE 104 CONTINUE RETURN 105 DO 107 I=2,IDO,2 CDIR$ IVDEP DO 106 K=1,L1 TI1 = CC(I,1,K)-CC(I,3,K) TI2 = CC(I,1,K)+CC(I,3,K) TI3 = CC(I,2,K)+CC(I,4,K) TR4 = CC(I,2,K)-CC(I,4,K) TR1 = CC(I-1,1,K)-CC(I-1,3,K) TR2 = CC(I-1,1,K)+CC(I-1,3,K) TI4 = CC(I-1,4,K)-CC(I-1,2,K) TR3 = CC(I-1,2,K)+CC(I-1,4,K) CH(I-1,K,1) = TR2+TR3 CR3 = TR2-TR3 CH(I,K,1) = TI2+TI3 CI3 = TI2-TI3 CR2 = TR1+TR4 CR4 = TR1-TR4 CI2 = TI1+TI4 CI4 = TI1-TI4 CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2 CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2 CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3 CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3 CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4 CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4 106 CONTINUE 107 CONTINUE RETURN END SUBROUTINE DPASF5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4) C***BEGIN PROLOGUE DPASF5 C***REFER TO DCFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DPASF5 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*) C***FIRST EXECUTABLE STATEMENT DPASF5 PI = 4.0D0*ATAN(1.0D0) TR11 = SIN(.1D0*PI) TI11 = -SIN(.4D0*PI) TR12 = -SIN(.3D0*PI) TI12 = -SIN(.2D0*PI) IF (IDO .NE. 2) GO TO 102 DO 101 K=1,L1 TI5 = CC(2,2,K)-CC(2,5,K) TI2 = CC(2,2,K)+CC(2,5,K) TI4 = CC(2,3,K)-CC(2,4,K) TI3 = CC(2,3,K)+CC(2,4,K) TR5 = CC(1,2,K)-CC(1,5,K) TR2 = CC(1,2,K)+CC(1,5,K) TR4 = CC(1,3,K)-CC(1,4,K) TR3 = CC(1,3,K)+CC(1,4,K) CH(1,K,1) = CC(1,1,K)+TR2+TR3 CH(2,K,1) = CC(2,1,K)+TI2+TI3 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3 CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3 CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3 CR5 = TI11*TR5+TI12*TR4 CI5 = TI11*TI5+TI12*TI4 CR4 = TI12*TR5-TI11*TR4 CI4 = TI12*TI5-TI11*TI4 CH(1,K,2) = CR2-CI5 CH(1,K,5) = CR2+CI5 CH(2,K,2) = CI2+CR5 CH(2,K,3) = CI3+CR4 CH(1,K,3) = CR3-CI4 CH(1,K,4) = CR3+CI4 CH(2,K,4) = CI3-CR4 CH(2,K,5) = CI2-CR5 101 CONTINUE RETURN 102 IF(IDO/2.LT.L1) GO TO 105 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=2,IDO,2 TI5 = CC(I,2,K)-CC(I,5,K) TI2 = CC(I,2,K)+CC(I,5,K) TI4 = CC(I,3,K)-CC(I,4,K) TI3 = CC(I,3,K)+CC(I,4,K) TR5 = CC(I-1,2,K)-CC(I-1,5,K) TR2 = CC(I-1,2,K)+CC(I-1,5,K) TR4 = CC(I-1,3,K)-CC(I-1,4,K) TR3 = CC(I-1,3,K)+CC(I-1,4,K) CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 CH(I,K,1) = CC(I,1,K)+TI2+TI3 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 CR5 = TI11*TR5+TI12*TR4 CI5 = TI11*TI5+TI12*TI4 CR4 = TI12*TR5-TI11*TR4 CI4 = TI12*TI5-TI11*TI4 DR3 = CR3-CI4 DR4 = CR3+CI4 DI3 = CI3+CR4 DI4 = CI3-CR4 DR5 = CR2+CI5 DR2 = CR2-CI5 DI5 = CI2-CR5 DI2 = CI2+CR5 CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4 CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4 CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5 CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5 103 CONTINUE 104 CONTINUE RETURN 105 DO 107 I=2,IDO,2 CDIR$ IVDEP DO 106 K=1,L1 TI5 = CC(I,2,K)-CC(I,5,K) TI2 = CC(I,2,K)+CC(I,5,K) TI4 = CC(I,3,K)-CC(I,4,K) TI3 = CC(I,3,K)+CC(I,4,K) TR5 = CC(I-1,2,K)-CC(I-1,5,K) TR2 = CC(I-1,2,K)+CC(I-1,5,K) TR4 = CC(I-1,3,K)-CC(I-1,4,K) TR3 = CC(I-1,3,K)+CC(I-1,4,K) CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 CH(I,K,1) = CC(I,1,K)+TI2+TI3 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 CR5 = TI11*TR5+TI12*TR4 CI5 = TI11*TI5+TI12*TI4 CR4 = TI12*TR5-TI11*TR4 CI4 = TI12*TI5-TI11*TI4 DR3 = CR3-CI4 DR4 = CR3+CI4 DI3 = CI3+CR4 DI4 = CI3-CR4 DR5 = CR2+CI5 DR2 = CR2-CI5 DI5 = CI2-CR5 DI2 = CI2+CR5 CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2 CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2 CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3 CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3 CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4 CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4 CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5 CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5 106 CONTINUE 107 CONTINUE RETURN END SUBROUTINE DCFFTI(N,WSAVE) C***BEGIN PROLOGUE DCFFTI C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 860115 (YYMMDD) C***CATEGORY NO. J1A2 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE Initialize for DCFFTF and DCFFTB. C***DESCRIPTION C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1988 C C Subroutine DCFFTI initializes the array WSAVE which is used in C both DCFFTF and DCFFTB. The prime factorization of N together with C a tabulation of the trigonometric functions are computed and C stored in WSAVE. C C Input Parameter C C N the length of the sequence to be transformed C C Output Parameter C C WSAVE a work array which must be dimensioned at least 4*N+15. C The same work array can be used for both DCFFTF and DCFFTB C as long as N remains unchanged. Different WSAVE arrays C are required for different values of N. The contents of C WSAVE must not be changed between calls of DCFFTF or DCFFTB. C***REFERENCES (NONE) C***ROUTINES CALLED DCFTI1 C***END PROLOGUE DCFFTI IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WSAVE(*) C***FIRST EXECUTABLE STATEMENT DCFFTI IF (N .EQ. 1) RETURN IW1 = N+N+1 IW2 = IW1+N+N CALL DCFTI1 (N,WSAVE(IW1),WSAVE(IW2)) RETURN END SUBROUTINE DCFTI1(N,WA,IFAC) C***BEGIN PROLOGUE DCFTI1 C***REFER TO DCFFTI C***ROUTINES CALLED (NONE) C***END PROLOGUE DCFTI1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/ C***FIRST EXECUTABLE STATEMENT DCFTI1 NL = N NF = 0 J = 0 101 J = J+1 IF (J-4) 102,102,103 102 NTRY = NTRYH(J) GO TO 104 103 NTRY = NTRY+2 104 NQ = NL/NTRY NR = NL-NTRY*NQ IF (NR) 101,105,101 105 NF = NF+1 IFAC(NF+2) = NTRY NL = NQ IF (NTRY .NE. 2) GO TO 107 IF (NF .EQ. 1) GO TO 107 DO 106 I=2,NF IB = NF-I+2 IFAC(IB+2) = IFAC(IB+1) 106 CONTINUE IFAC(3) = 2 107 IF (NL .NE. 1) GO TO 104 IFAC(1) = N IFAC(2) = NF TPI = 8.0D0*ATAN(1.0D0) ARGH = TPI/DBLE(N) I = 2 L1 = 1 DO 110 K1=1,NF IP = IFAC(K1+2) LD = 0 L2 = L1*IP IDO = N/L2 IDOT = IDO+IDO+2 IPM = IP-1 DO 109 J=1,IPM I1 = I WA(I-1) = 1.0D0 WA(I) = 0.0D0 LD = LD+L1 FI = 0.0D0 ARGLD = DBLE(LD)*ARGH DO 108 II=4,IDOT,2 I = I+2 FI = FI+1.0D0 ARG = FI*ARGLD WA(I-1) = COS(ARG) WA(I) = SIN(ARG) 108 CONTINUE IF (IP .LE. 5) GO TO 109 WA(I1-1) = WA(I-1) WA(I1) = WA(I) 109 CONTINUE L1 = L2 110 CONTINUE RETURN END SUBROUTINE DCFT2D(N,F,LDF,W,FORWD) C***BEGIN PROLOGUE DCFT2D C***DATE WRITTEN 870811 (YYMMDD) C***REVISION DATE 870811 (YYMMDD) C***CATEGORY NO. J1B C***KEYWORDS TWO DIMENSIONAL FOURIER TRANSFORM, FFT C***AUTHOR KAHANER, DAVID K., (NBS) C***PURPOSE Two dimensional complex fast Fourier transform. C***DESCRIPTION C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1988 C Two dimensional fast Fourier transform, forward or backward C of complex N*N matrix F. C C Input: C N: (INTEGER) Number of rows and columns in the matrix F to be C transformed. You must set N > 0, NOT CHECKED. C F: (COMPLEX*16) Array of N*N complex values to be transformed. C This array is overwritten on output. C LDF: (INTEGER) Leading (first) dimension of the complex array F C in the subroutine that calls DCFT2D. For example, C if you declare F by either C COMPLEX*16 F(0:15,0:20) OR F(16,21) then C set LDF=16. You must have LDF >= N, NOT CHECKED. C W: (COMPLEX*16) Array for internal use as work storage. C Must be dimensioned by calling program to be C at least 6N+15 D. P. words or 3N+8 COMPLEX*16 words. C FORWD: (LOGICAL) Direction of transform. Set to .TRUE. for C forward transform, set to .FALSE. for backward C transform. C C C Output: C F: (COMPLEX*16) Forward or reverse transformed input matrix. C Output is unscaled, that is, a call to DCFT2D C with FORWD=.TRUE. followed by a call to DCFT2D C with FORWD=.FALSE. returns original data C multiplied by N*N. C C C Remark: C For some applications it is desirable to have the transform scaled so C the center of the N by N frequency square corresponds to zero C frequency. The user can do this replacing the original input data C F(I,J) by F(I,J)*(-1.)**(I+J), I,J =0,...,N-1.0D0 C C C***REFERENCES (NONE) C***ROUTINES CALLED DCFFTI, DCFFTF, DCFFTB C***END PROLOGUE DCFT2D COMPLEX*16 F(0:LDF-1,0:*), W(0:*) LOGICAL FORWD C C NOTE: DOUBLE PRECISION COMPLEX IS NOT STANDARD FORTRAN. C SOME COMPILER VENDORS MAY FAIL TO IMPLEMENT IT C OR MAY IMPLEMENT IN A DIFFERENT MANNER. IF YOUR C COMPILER DOES NOT ACCEPT THE COMPLEX*16 STATEMENT C WE RECOMMEND THAT YOU CHECK WITH YOUR VENDOR AND C MODIFY THE ABOVE ACCORDINGLY. C C***FIRST EXECUTABLE STATEMENT DCFT2D C Find transform of each row, then of each column C First row transforms CALL DCFFTI(N,W(N)) DO 10 I=0,N-1 C Place row in beginning of array W DO 5 J=0,N-1 W(J) = F(I,J) 5 CONTINUE C Compute fft of row IF(FORWD) THEN CALL DCFFTF(N,W,W(N)) ELSE CALL DCFFTB(N,W,W(N)) ENDIF C Copy back to row of f DO 6 J=0,N-1 F(I,J)=W(J) 6 CONTINUE 10 CONTINUE C Column transforms DO 20 J=0,N-1 C Pass column of F to DCFFTF or DCFFTB by passing F(0,J). C Similarly DCFFTF or DCFFTB places results back there. IF(FORWD) THEN CALL DCFFTF(N, F(0,J), W(N)) ELSE CALL DCFFTB(N, F(0,J), W(N)) ENDIF 20 CONTINUE RETURN END SUBROUTINE DEZFTB(N,R,AZERO,A,B,WSAVE) C***BEGIN PROLOGUE DEZFTB C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 860115 (YYMMDD) C***CATEGORY NO. J1A1 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE A Simplified double precision, periodic, C backward transform C***DESCRIPTION C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1988 C C Subroutine DEZFTB computes a d. p. periodic sequence from its C Fourier coefficients (Fourier synthesis). The transform is C defined below at Output Parameter R. DEZFTB is a simplified C but slower version of DRFFTB. C C Input Parameters C C N the length of the output array R. The method is most C efficient when N is the product of small primes. C C AZERO the constant Fourier coefficient C C A,B arrays which contain the remaining Fourier coefficients. C These arrays are not destroyed. C C The length of these arrays depends on whether N is even or C odd. C C If N is even, N/2 locations are required. C If N is odd, (N-1)/2 locations are required C C WSAVE a work array which must be dimensioned at least 3*N+15 C in the program that calls DEZFTB. The WSAVE array must be C initialized by calling subroutine DEZFTI(N,WSAVE), and a C different WSAVE array must be used for each different C value of N. This initialization does not have to be C repeated so long as N remains unchanged. Thus subsequent C transforms can be obtained faster than the first. C The same WSAVE array can be used by EZFFTF and DEZFTB. C C C Output Parameters C C R if N is even, define KMAX=N/2 C if N is odd, define KMAX=(N-1)/2 C C Then for I=1,...,N C C R(I)=AZERO plus the sum from K=1 to K=KMAX of C C A(K)*COS(K*(I-1)*2*PI/N)+B(K)*SIN(K*(I-1)*2*PI/N) C C ********************* Complex Notation ************************** C C For J=1,...,N C C R(J) equals the sum from K=-KMAX to K=KMAX of C C C(K)*EXP(I*K*(J-1)*2*PI/N) C C where C C C(K) = .5D0*CMPLX(A(K),-B(K)) for K=1,...,KMAX C C C(-K) = CONJG(C(K)) C C C(0) = AZERO C C and I=SQRT(-1) C C *************** Amplitude - Phase Notation *********************** C C For I=1,...,N C C R(I) equals AZERO plus the sum from K=1 to K=KMAX of C C ALPHA(K)*COS(K*(I-1)*2*PI/N+BETA(K)) C C where C C ALPHA(K) = SQRT(A(K)*A(K)+B(K)*B(K)) C C COS(BETA(K))=A(K)/ALPHA(K) C C SIN(BETA(K))=-B(K)/ALPHA(K) C C * * C * References * C * * C * 1. P.N. Swarztrauber, Vectorizing the FFTs, in Parallel * C * Computations (G. Rodrigue, ed.), Academic Press, 1982, * C * pp. 51-83. * C * 2. B.L. Buzbee, The SLATEC Common Math Library, in Sources * C * and Development of Mathematical Software (W. Cowell, ed.), * C * Prentice-Hall, 1984, pp. 302-318. * C * * C ********************************************************************* C C***REFERENCES (NONE) C***ROUTINES CALLED DRFFTB C***END PROLOGUE DEZFTB IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION R(*) ,A(*) ,B(*) ,WSAVE(*) C***FIRST EXECUTABLE STATEMENT DEZFTB IF (N-2) 101,102,103 101 R(1) = AZERO RETURN 102 R(1) = AZERO+A(1) R(2) = AZERO-A(1) RETURN 103 NS2 = (N-1)/2 DO 104 I=1,NS2 R(2*I) = .5D0*A(I) R(2*I+1) = -.5D0*B(I) 104 CONTINUE R(1) = AZERO IF (MOD(N,2) .EQ. 0) R(N) = A(NS2+1) CALL DRFFTB (N,R,WSAVE(N+1)) RETURN END SUBROUTINE DRFFTB(N,R,WSAVE) C***BEGIN PROLOGUE DRFFTB C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 890413 (YYMMDD) C***CATEGORY NO. J1A1 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE Backward transform of a dp coefficient array. C***DESCRIPTION C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1989 C C Subroutine DRFFTB computes the dp perodic sequence from its C Fourier coefficients (Fourier synthesis). The transform is defined C below at output parameter R. C C Input Parameters C C N the length of the array R to be transformed. The method C is most efficient when N is a product of small primes. C N may change so long as different work arrays are provided. C C R a dp array of length N which contains the sequence C to be transformed C C WSAVE a dp work array which must be dimensioned at least 2*N+15 C in the program that calls DRFFTB. The WSAVE array must be C initialized by calling subroutine DRFFTI(N,WSAVE), and a C different WSAVE array must be used for each different C value of N. This initialization does not have to be C repeated so long as N remains unchanged. Thus subsequent C transforms can be obtained faster than the first. C The same WSAVE array can be used by DRFFTF and DRFFTB. C C C Output Parameters C C R For N even and For I = 1,...,N C C R(I) = R(1)+(-1)**(I-1)*R(N) C C plus the sum from K=2 to K=N/2 of C C 2.*R(2*K-2)*COS((K-1)*(I-1)*2*PI/N) C C -2.*R(2*K-1)*SIN((K-1)*(I-1)*2*PI/N) C C For N odd and For I = 1,...,N C C R(I) = R(1) plus the sum from K=2 to K=(N+1)/2 of C C 2.*R(2*K-2)*COS((K-1)*(I-1)*2*PI/N) C C -2.*R(2*K-1)*SIN((K-1)*(I-1)*2*PI/N) C C ***** Note: C This transform is unnormalized since a call of DRFFTF C followed by a call of DRFFTB will multiply the input C sequence by N. C C WSAVE contains results which must not be destroyed between C calls of DRFFTB or DRFFTF. C C * * C * 1. P.N. Swarztrauber, Vectorizing the FFTs, in Parallel * C * Computations (G. Rodrigue, ed.), Academic Press, 1982, * C * pp. 51-83. * C * 2. B.L. Buzbee, The SLATEC Common Math Library, in Sources * C * and Development of Mathematical Software (W. Cowell, ed.), * C * Prentice-Hall, 1984, pp. 302-318. * C * * C ********************************************************************* C C***REFERENCES (NONE) C***ROUTINES CALLED DRFTB1 C***END PROLOGUE DRFFTB IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION R(*) ,WSAVE(*) C***FIRST EXECUTABLE STATEMENT DRFFTB IF (N .EQ. 1) RETURN CALL DRFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1)) RETURN END SUBROUTINE DRFTB1(N,C,CH,WA,IFAC) C***BEGIN PROLOGUE DRFTB1 C***REFER TO DRFFTB C***ROUTINES CALLED DRADB2,DRADB3,DRADB4,DRADB5,DRADBG C***END PROLOGUE DRFTB1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) C***FIRST EXECUTABLE STATEMENT DRFTB1 NF = IFAC(2) NA = 0 L1 = 1 IW = 1 DO 116 K1=1,NF IP = IFAC(K1+2) L2 = IP*L1 IDO = N/L2 IDL1 = IDO*L1 IF (IP .NE. 4) GO TO 103 IX2 = IW+IDO IX3 = IX2+IDO IF (NA .NE. 0) GO TO 101 CALL DRADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) GO TO 102 101 CALL DRADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) 102 NA = 1-NA GO TO 115 103 IF (IP .NE. 2) GO TO 106 IF (NA .NE. 0) GO TO 104 CALL DRADB2 (IDO,L1,C,CH,WA(IW)) GO TO 105 104 CALL DRADB2 (IDO,L1,CH,C,WA(IW)) 105 NA = 1-NA GO TO 115 106 IF (IP .NE. 3) GO TO 109 IX2 = IW+IDO IF (NA .NE. 0) GO TO 107 CALL DRADB3 (IDO,L1,C,CH,WA(IW),WA(IX2)) GO TO 108 107 CALL DRADB3 (IDO,L1,CH,C,WA(IW),WA(IX2)) 108 NA = 1-NA GO TO 115 109 IF (IP .NE. 5) GO TO 112 IX2 = IW+IDO IX3 = IX2+IDO IX4 = IX3+IDO IF (NA .NE. 0) GO TO 110 CALL DRADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 111 110 CALL DRADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) 111 NA = 1-NA GO TO 115 112 IF (NA .NE. 0) GO TO 113 CALL DRADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) GO TO 114 113 CALL DRADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) 114 IF (IDO .EQ. 1) NA = 1-NA 115 L1 = L2 IW = IW+(IP-1)*IDO 116 CONTINUE IF (NA .EQ. 0) RETURN DO 117 I=1,N C(I) = CH(I) 117 CONTINUE RETURN END SUBROUTINE DRADB2(IDO,L1,CC,CH,WA1) C***BEGIN PROLOGUE DRADB2 C***REFER TO DRFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADB2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) , 1 WA1(*) C***FIRST EXECUTABLE STATEMENT DRADB2 DO 101 K=1,L1 CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K) CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 IF((IDO-1)/2.LT.L1) GO TO 108 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=3,IDO,2 IC = IDP2-I CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K) TR2 = CC(I-1,1,K)-CC(IC-1,2,K) CH(I,K,1) = CC(I,1,K)-CC(IC,2,K) TI2 = CC(I,1,K)+CC(IC,2,K) CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2 CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2 103 CONTINUE 104 CONTINUE GO TO 111 108 DO 110 I=3,IDO,2 IC = IDP2-I CDIR$ IVDEP DO 109 K=1,L1 CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K) TR2 = CC(I-1,1,K)-CC(IC-1,2,K) CH(I,K,1) = CC(I,1,K)-CC(IC,2,K) TI2 = CC(I,1,K)+CC(IC,2,K) CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2 CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2 109 CONTINUE 110 CONTINUE 111 IF (MOD(IDO,2) .EQ. 1) RETURN 105 DO 106 K=1,L1 CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K) CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K)) 106 CONTINUE 107 RETURN END SUBROUTINE DRADB3(IDO,L1,CC,CH,WA1,WA2) C***BEGIN PROLOGUE DRADB3 C***REFER TO DRFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADB3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) , 1 WA1(*) ,WA2(*) C***FIRST EXECUTABLE STATEMENT DRADB3 TAUR = -.5D0 TAUI = .5D0*SQRT(3.0D0) DO 101 K=1,L1 TR2 = CC(IDO,2,K)+CC(IDO,2,K) CR2 = CC(1,1,K)+TAUR*TR2 CH(1,K,1) = CC(1,1,K)+TR2 CI3 = TAUI*(CC(1,3,K)+CC(1,3,K)) CH(1,K,2) = CR2-CI3 CH(1,K,3) = CR2+CI3 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 IF((IDO-1)/2.LT.L1) GO TO 104 DO 103 K=1,L1 CDIR$ IVDEP DO 102 I=3,IDO,2 IC = IDP2-I TR2 = CC(I-1,3,K)+CC(IC-1,2,K) CR2 = CC(I-1,1,K)+TAUR*TR2 CH(I-1,K,1) = CC(I-1,1,K)+TR2 TI2 = CC(I,3,K)-CC(IC,2,K) CI2 = CC(I,1,K)+TAUR*TI2 CH(I,K,1) = CC(I,1,K)+TI2 CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K)) CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K)) DR2 = CR2-CI3 DR3 = CR2+CI3 DI2 = CI2+CR3 DI3 = CI2-CR3 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3 102 CONTINUE 103 CONTINUE RETURN 104 DO 106 I=3,IDO,2 IC = IDP2-I CDIR$ IVDEP DO 105 K=1,L1 TR2 = CC(I-1,3,K)+CC(IC-1,2,K) CR2 = CC(I-1,1,K)+TAUR*TR2 CH(I-1,K,1) = CC(I-1,1,K)+TR2 TI2 = CC(I,3,K)-CC(IC,2,K) CI2 = CC(I,1,K)+TAUR*TI2 CH(I,K,1) = CC(I,1,K)+TI2 CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K)) CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K)) DR2 = CR2-CI3 DR3 = CR2+CI3 DI2 = CI2+CR3 DI3 = CI2-CR3 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3 105 CONTINUE 106 CONTINUE RETURN END SUBROUTINE DRADB4(IDO,L1,CC,CH,WA1,WA2,WA3) C***BEGIN PROLOGUE DRADB4 C***REFER TO DRFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADB4 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) , 1 WA1(*) ,WA2(*) ,WA3(*) C***FIRST EXECUTABLE STATEMENT DRADB4 SQRT2 = SQRT(2.0D0) DO 101 K=1,L1 TR1 = CC(1,1,K)-CC(IDO,4,K) TR2 = CC(1,1,K)+CC(IDO,4,K) TR3 = CC(IDO,2,K)+CC(IDO,2,K) TR4 = CC(1,3,K)+CC(1,3,K) CH(1,K,1) = TR2+TR3 CH(1,K,2) = TR1-TR4 CH(1,K,3) = TR2-TR3 CH(1,K,4) = TR1+TR4 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 IF((IDO-1)/2.LT.L1) GO TO 108 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=3,IDO,2 IC = IDP2-I TI1 = CC(I,1,K)+CC(IC,4,K) TI2 = CC(I,1,K)-CC(IC,4,K) TI3 = CC(I,3,K)-CC(IC,2,K) TR4 = CC(I,3,K)+CC(IC,2,K) TR1 = CC(I-1,1,K)-CC(IC-1,4,K) TR2 = CC(I-1,1,K)+CC(IC-1,4,K) TI4 = CC(I-1,3,K)-CC(IC-1,2,K) TR3 = CC(I-1,3,K)+CC(IC-1,2,K) CH(I-1,K,1) = TR2+TR3 CR3 = TR2-TR3 CH(I,K,1) = TI2+TI3 CI3 = TI2-TI3 CR2 = TR1-TR4 CR4 = TR1+TR4 CI2 = TI1+TI4 CI4 = TI1-TI4 CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2 CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2 CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3 CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3 CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4 CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4 103 CONTINUE 104 CONTINUE GO TO 111 108 DO 110 I=3,IDO,2 IC = IDP2-I CDIR$ IVDEP DO 109 K=1,L1 TI1 = CC(I,1,K)+CC(IC,4,K) TI2 = CC(I,1,K)-CC(IC,4,K) TI3 = CC(I,3,K)-CC(IC,2,K) TR4 = CC(I,3,K)+CC(IC,2,K) TR1 = CC(I-1,1,K)-CC(IC-1,4,K) TR2 = CC(I-1,1,K)+CC(IC-1,4,K) TI4 = CC(I-1,3,K)-CC(IC-1,2,K) TR3 = CC(I-1,3,K)+CC(IC-1,2,K) CH(I-1,K,1) = TR2+TR3 CR3 = TR2-TR3 CH(I,K,1) = TI2+TI3 CI3 = TI2-TI3 CR2 = TR1-TR4 CR4 = TR1+TR4 CI2 = TI1+TI4 CI4 = TI1-TI4 CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2 CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2 CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3 CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3 CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4 CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4 109 CONTINUE 110 CONTINUE 111 IF (MOD(IDO,2) .EQ. 1) RETURN 105 DO 106 K=1,L1 TI1 = CC(1,2,K)+CC(1,4,K) TI2 = CC(1,4,K)-CC(1,2,K) TR1 = CC(IDO,1,K)-CC(IDO,3,K) TR2 = CC(IDO,1,K)+CC(IDO,3,K) CH(IDO,K,1) = TR2+TR2 CH(IDO,K,2) = SQRT2*(TR1-TI1) CH(IDO,K,3) = TI2+TI2 CH(IDO,K,4) = -SQRT2*(TR1+TI1) 106 CONTINUE 107 RETURN END SUBROUTINE DRADB5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4) C***BEGIN PROLOGUE DRADB5 C***REFER TO DRFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADB5 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) , 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*) C***FIRST EXECUTABLE STATEMENT DRADB5 PI = 4.0D0*ATAN(1.0D0) TR11 = SIN(.1D0*PI) TI11 = SIN(.4D0*PI) TR12 = -SIN(.3D0*PI) TI12 = SIN(.2D0*PI) DO 101 K=1,L1 TI5 = CC(1,3,K)+CC(1,3,K) TI4 = CC(1,5,K)+CC(1,5,K) TR2 = CC(IDO,2,K)+CC(IDO,2,K) TR3 = CC(IDO,4,K)+CC(IDO,4,K) CH(1,K,1) = CC(1,1,K)+TR2+TR3 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3 CI5 = TI11*TI5+TI12*TI4 CI4 = TI12*TI5-TI11*TI4 CH(1,K,2) = CR2-CI5 CH(1,K,3) = CR3-CI4 CH(1,K,4) = CR3+CI4 CH(1,K,5) = CR2+CI5 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 IF((IDO-1)/2.LT.L1) GO TO 104 DO 103 K=1,L1 CDIR$ IVDEP DO 102 I=3,IDO,2 IC = IDP2-I TI5 = CC(I,3,K)+CC(IC,2,K) TI2 = CC(I,3,K)-CC(IC,2,K) TI4 = CC(I,5,K)+CC(IC,4,K) TI3 = CC(I,5,K)-CC(IC,4,K) TR5 = CC(I-1,3,K)-CC(IC-1,2,K) TR2 = CC(I-1,3,K)+CC(IC-1,2,K) TR4 = CC(I-1,5,K)-CC(IC-1,4,K) TR3 = CC(I-1,5,K)+CC(IC-1,4,K) CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 CH(I,K,1) = CC(I,1,K)+TI2+TI3 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 CR5 = TI11*TR5+TI12*TR4 CI5 = TI11*TI5+TI12*TI4 CR4 = TI12*TR5-TI11*TR4 CI4 = TI12*TI5-TI11*TI4 DR3 = CR3-CI4 DR4 = CR3+CI4 DI3 = CI3+CR4 DI4 = CI3-CR4 DR5 = CR2+CI5 DR2 = CR2-CI5 DI5 = CI2-CR5 DI2 = CI2+CR5 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3 CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4 CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4 CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5 CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5 102 CONTINUE 103 CONTINUE RETURN 104 DO 106 I=3,IDO,2 IC = IDP2-I CDIR$ IVDEP DO 105 K=1,L1 TI5 = CC(I,3,K)+CC(IC,2,K) TI2 = CC(I,3,K)-CC(IC,2,K) TI4 = CC(I,5,K)+CC(IC,4,K) TI3 = CC(I,5,K)-CC(IC,4,K) TR5 = CC(I-1,3,K)-CC(IC-1,2,K) TR2 = CC(I-1,3,K)+CC(IC-1,2,K) TR4 = CC(I-1,5,K)-CC(IC-1,4,K) TR3 = CC(I-1,5,K)+CC(IC-1,4,K) CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3 CH(I,K,1) = CC(I,1,K)+TI2+TI3 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3 CR5 = TI11*TR5+TI12*TR4 CI5 = TI11*TI5+TI12*TI4 CR4 = TI12*TR5-TI11*TR4 CI4 = TI12*TI5-TI11*TI4 DR3 = CR3-CI4 DR4 = CR3+CI4 DI3 = CI3+CR4 DI4 = CI3-CR4 DR5 = CR2+CI5 DR2 = CR2-CI5 DI5 = CI2-CR5 DI2 = CI2+CR5 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3 CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4 CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4 CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5 CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5 105 CONTINUE 106 CONTINUE RETURN END SUBROUTINE DRADBG(IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) C***BEGIN PROLOGUE DRADBG C***REFER TO DRFFTB C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADBG IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,C2(IDL1,IP), 2 CH2(IDL1,IP) ,WA(*) C***FIRST EXECUTABLE STATEMENT DRADBG TPI = 8.0D0*ATAN(1.0D0) ARG = TPI/DBLE(IP) DCP = COS(ARG) DSP = SIN(ARG) IDP2 = IDO+2 NBD = (IDO-1)/2 IPP2 = IP+2 IPPH = (IP+1)/2 IF (IDO .LT. L1) GO TO 103 DO 102 K=1,L1 DO 101 I=1,IDO CH(I,K,1) = CC(I,1,K) 101 CONTINUE 102 CONTINUE GO TO 106 103 DO 105 I=1,IDO DO 104 K=1,L1 CH(I,K,1) = CC(I,1,K) 104 CONTINUE 105 CONTINUE 106 DO 108 J=2,IPPH JC = IPP2-J J2 = J+J DO 107 K=1,L1 CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K) CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K) 107 CONTINUE 108 CONTINUE IF (IDO .EQ. 1) GO TO 116 IF (NBD .LT. L1) GO TO 112 DO 111 J=2,IPPH JC = IPP2-J DO 110 K=1,L1 CDIR$ IVDEP DO 109 I=3,IDO,2 IC = IDP2-I CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K) CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K) CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K) CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K) 109 CONTINUE 110 CONTINUE 111 CONTINUE GO TO 116 112 DO 115 J=2,IPPH JC = IPP2-J CDIR$ IVDEP DO 114 I=3,IDO,2 IC = IDP2-I DO 113 K=1,L1 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K) CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K) CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K) CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K) 113 CONTINUE 114 CONTINUE 115 CONTINUE 116 AR1 = 1.0D0 AI1 = 0.0D0 DO 120 L=2,IPPH LC = IPP2-L AR1H = DCP*AR1-DSP*AI1 AI1 = DCP*AI1+DSP*AR1 AR1 = AR1H DO 117 IK=1,IDL1 C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2) C2(IK,LC) = AI1*CH2(IK,IP) 117 CONTINUE DC2 = AR1 DS2 = AI1 AR2 = AR1 AI2 = AI1 DO 119 J=3,IPPH JC = IPP2-J AR2H = DC2*AR2-DS2*AI2 AI2 = DC2*AI2+DS2*AR2 AR2 = AR2H DO 118 IK=1,IDL1 C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J) C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC) 118 CONTINUE 119 CONTINUE 120 CONTINUE DO 122 J=2,IPPH DO 121 IK=1,IDL1 CH2(IK,1) = CH2(IK,1)+CH2(IK,J) 121 CONTINUE 122 CONTINUE DO 124 J=2,IPPH JC = IPP2-J DO 123 K=1,L1 CH(1,K,J) = C1(1,K,J)-C1(1,K,JC) CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC) 123 CONTINUE 124 CONTINUE IF (IDO .EQ. 1) GO TO 132 IF (NBD .LT. L1) GO TO 128 DO 127 J=2,IPPH JC = IPP2-J DO 126 K=1,L1 CDIR$ IVDEP DO 125 I=3,IDO,2 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC) CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC) CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC) CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC) 125 CONTINUE 126 CONTINUE 127 CONTINUE GO TO 132 128 DO 131 J=2,IPPH JC = IPP2-J DO 130 I=3,IDO,2 DO 129 K=1,L1 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC) CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC) CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC) CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC) 129 CONTINUE 130 CONTINUE 131 CONTINUE 132 CONTINUE IF (IDO .EQ. 1) RETURN DO 133 IK=1,IDL1 C2(IK,1) = CH2(IK,1) 133 CONTINUE DO 135 J=2,IP DO 134 K=1,L1 C1(1,K,J) = CH(1,K,J) 134 CONTINUE 135 CONTINUE IF (NBD .GT. L1) GO TO 139 IS = -IDO DO 138 J=2,IP IS = IS+IDO IDIJ = IS DO 137 I=3,IDO,2 IDIJ = IDIJ+2 DO 136 K=1,L1 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) 136 CONTINUE 137 CONTINUE 138 CONTINUE GO TO 143 139 IS = -IDO DO 142 J=2,IP IS = IS+IDO DO 141 K=1,L1 IDIJ = IS CDIR$ IVDEP DO 140 I=3,IDO,2 IDIJ = IDIJ+2 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J) C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J) 140 CONTINUE 141 CONTINUE 142 CONTINUE 143 RETURN END SUBROUTINE DEZFTF(N,R,AZERO,A,B,WSAVE) C***BEGIN PROLOGUE DEZFTF C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 860115 (YYMMDD) C***CATEGORY NO. J1A1 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE A simplified d. p., periodic, forward transform C***DESCRIPTION C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1988 C C Subroutine DEZFTF computes the Fourier coefficients of a d. p. C perodic sequence (Fourier analysis). The transform is defined C below at Output Parameters AZERO, A and B. DEZFTF is a simplified C but slower version of DRFFTF. C C Input Parameters C C N the length of the array R to be transformed. The method C is must efficient when N is the product of small primes. C C R a d. p. array of length N which contains the sequence C to be transformed. R is not destroyed. C C C WSAVE a work array which must be dimensioned at least 3*N+15 C in the program that calls DEZFTF. The WSAVE array must be C initialized by calling subroutine DEZFTI(N,WSAVE), and a C different WSAVE array must be used for each different C value of N. This initialization does not have to be C repeated so long as N remains unchanged. Thus subsequent C transforms can be obtained faster than the first. C The same WSAVE array can be used by DEZFTF and EZFFTB. C C Output Parameters C C AZERO the sum from I=1 to I=N of R(I)/N C C A,B for N even B(N/2)=0.0D0 and A(N/2) is the sum from I=1 to C I=N of (-1)**(I-1)*R(I)/N C C for N even define KMAX=N/2-1 C for N odd define KMAX=(N-1)/2 C C then for k=1,...,KMAX C C A(K) equals the sum from I=1 to I=N of C C 2.0D0/N*R(I)*COS(K*(I-1)*2*PI/N) C C B(K) equals the sum from I=1 to I=N of C C 2.0D0/N*R(I)*SIN(K*(I-1)*2*PI/N) C C * * C * References * C * * C * 1. P.N. Swarztrauber, Vectorizing the FFTs, in Parallel * C * Computations (G. Rodrigue, ed.), Academic Press, 1982, * C * pp. 51-83. * C * 2. B.L. Buzbee, The SLATEC Common Math Library, in Sources * C * and Development of Mathematical Software (W. Cowell, ed.), * C * Prentice-Hall, 1984, pp. 302-318. * C * * C C***REFERENCES (NONE) C***ROUTINES CALLED DRFFTF C***END PROLOGUE DEZFTF IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION R(*) ,A(*) ,B(*) ,WSAVE(*) C***FIRST EXECUTABLE STATEMENT DEZFTF IF (N-2) 101,102,103 101 AZERO = R(1) RETURN 102 AZERO = 0.5D0*(R(1)+R(2)) A(1) = 0.5D0*(R(1)-R(2)) RETURN 103 DO 104 I=1,N WSAVE(I) = R(I) 104 CONTINUE CALL DRFFTF (N,WSAVE,WSAVE(N+1)) CF = 2.0D0/DBLE(N) CFM = -CF AZERO = 0.5D0*CF*WSAVE(1) NS2 = (N+1)/2 NS2M = NS2-1 DO 105 I=1,NS2M A(I) = CF*WSAVE(2*I) B(I) = CFM*WSAVE(2*I+1) 105 CONTINUE IF (MOD(N,2) .EQ. 0) THEN A(NS2) = 0.5D0*CF*WSAVE(N) B(NS2) = 0.0D0 ENDIF RETURN END SUBROUTINE DRFFTF(N,R,WSAVE) C***BEGIN PROLOGUE DRFFTF C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 890413 (YYMMDD) C***CATEGORY NO. J1A1 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE Forward transform of a dp, periodic sequence. C***DESCRIPTION C From the book "Numerical Methods and Software" C by D. Kahaner, C. Moler, S. Nash C Prentice Hall 1988 C C Subroutine DRFFTF computes the Fourier coefficients of a dp C perodic sequence (Fourier analysis). The transform is defined C below at output parameter R. C C Input Parameters C C N the length of the array R to be transformed. The method C is most efficient when N is a product of small primes. C N may change so long as different work arrays are provided C C R a dp array of length N which contains the sequence C to be transformed C C WSAVE a dp work array which must be dimensioned at least 2*N+15 C in the program that calls DRFFTF. The WSAVE array must be C initialized by calling subroutine DRFFTI(N,WSAVE), and a C different WSAVE array must be used for each different C value of N. This initialization does not have to be C repeated so long as N remains unchanged. Thus subsequent C transforms can be obtained faster than the first. C the same WSAVE array can be used by DRFFTF and DRFFTB. C C C Output Parameters C C R R(1) = the sum from I=1 to I=N of R(I) C C If N is even set L = N/2; if N is odd set L = (N+1)/2 C C then for K = 2,...,L C C R(2*K-2) = the sum from I = 1 to I = N of C C R(I)*COS((K-1)*(I-1)*2*PI/N) C C R(2*K-1) = the sum from I = 1 to I = N of C C -R(I)*SIN((K-1)*(I-1)*2*PI/N) C C If N is even C C R(N) = the sum from I = 1 to I = N of C C (-1)**(I-1)*R(I) C C ***** Note: C This transform is unnormalized since a call of DRFFTF C followed by a call of DRFFTB will multiply the input C sequence by N. C C WSAVE contains results which must not be destroyed between C calls of DRFFTF or DRFFTB. C C * * C * 1. P.N. Swarztrauber, Vectorizing the FFTs, in Parallel * C * Computations (G. Rodrigue, ed.), Academic Press, 1982, * C * pp. 51-83. * C * 2. B.L. Buzbee, The SLATEC Common Math Library, in Sources * C * and Development of Mathematical Software (W. Cowell, ed.), * C * Prentice-Hall, 1984, pp. 302-318. * C * * C ********************************************************************* C C***REFERENCES (NONE) C***ROUTINES CALLED DRFTF1 C***END PROLOGUE DRFFTF IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION R(*) ,WSAVE(*) C***FIRST EXECUTABLE STATEMENT DRFFTF IF (N .EQ. 1) RETURN CALL DRFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1)) RETURN END SUBROUTINE DRFTF1(N,C,CH,WA,IFAC) C***BEGIN PROLOGUE DRFTF1 C***REFER TO DRFTF C***ROUTINES CALLED DRADF2,DRADF3,DRADF4,DRADF5,DRADFG C***END PROLOGUE DRFTF1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(*) ,C(*) ,WA(*) ,IFAC(*) C***FIRST EXECUTABLE STATEMENT DRFTF1 NF = IFAC(2) NA = 1 L2 = N IW = N DO 111 K1=1,NF KH = NF-K1 IP = IFAC(KH+3) L1 = L2/IP IDO = N/L2 IDL1 = IDO*L1 IW = IW-(IP-1)*IDO NA = 1-NA IF (IP .NE. 4) GO TO 102 IX2 = IW+IDO IX3 = IX2+IDO IF (NA .NE. 0) GO TO 101 CALL DRADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3)) GO TO 110 101 CALL DRADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3)) GO TO 110 102 IF (IP .NE. 2) GO TO 104 IF (NA .NE. 0) GO TO 103 CALL DRADF2 (IDO,L1,C,CH,WA(IW)) GO TO 110 103 CALL DRADF2 (IDO,L1,CH,C,WA(IW)) GO TO 110 104 IF (IP .NE. 3) GO TO 106 IX2 = IW+IDO IF (NA .NE. 0) GO TO 105 CALL DRADF3 (IDO,L1,C,CH,WA(IW),WA(IX2)) GO TO 110 105 CALL DRADF3 (IDO,L1,CH,C,WA(IW),WA(IX2)) GO TO 110 106 IF (IP .NE. 5) GO TO 108 IX2 = IW+IDO IX3 = IX2+IDO IX4 = IX3+IDO IF (NA .NE. 0) GO TO 107 CALL DRADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 110 107 CALL DRADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4)) GO TO 110 108 IF (IDO .EQ. 1) NA = 1-NA IF (NA .NE. 0) GO TO 109 CALL DRADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW)) NA = 1 GO TO 110 109 CALL DRADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW)) NA = 0 110 L2 = L1 111 CONTINUE IF (NA .EQ. 1) RETURN DO 112 I=1,N C(I) = CH(I) 112 CONTINUE RETURN END SUBROUTINE DRADF2(IDO,L1,CC,CH,WA1) C***BEGIN PROLOGUE DRADF2 C***REFER TO DRFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADF2 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) , 1 WA1(*) C***FIRST EXECUTABLE STATEMENT DRADF2 DO 101 K=1,L1 CH(1,1,K) = CC(1,K,1)+CC(1,K,2) CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 IF((IDO-1)/2.LT.L1) GO TO 108 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=3,IDO,2 IC = IDP2-I TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) CH(I,1,K) = CC(I,K,1)+TI2 CH(IC,2,K) = TI2-CC(I,K,1) CH(I-1,1,K) = CC(I-1,K,1)+TR2 CH(IC-1,2,K) = CC(I-1,K,1)-TR2 103 CONTINUE 104 CONTINUE GO TO 111 108 DO 110 I=3,IDO,2 IC = IDP2-I CDIR$ IVDEP DO 109 K=1,L1 TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) CH(I,1,K) = CC(I,K,1)+TI2 CH(IC,2,K) = TI2-CC(I,K,1) CH(I-1,1,K) = CC(I-1,K,1)+TR2 CH(IC-1,2,K) = CC(I-1,K,1)-TR2 109 CONTINUE 110 CONTINUE 111 IF (MOD(IDO,2) .EQ. 1) RETURN 105 DO 106 K=1,L1 CH(1,2,K) = -CC(IDO,K,2) CH(IDO,1,K) = CC(IDO,K,1) 106 CONTINUE 107 RETURN END SUBROUTINE DRADF3(IDO,L1,CC,CH,WA1,WA2) C***BEGIN PROLOGUE DRADF3 C***REFER TO DRFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADF3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) , 1 WA1(*) ,WA2(*) C***FIRST EXECUTABLE STATEMENT DRADF3 TAUR = -0.5D0 TAUI = 0.5D0*SQRT(3.0D0) DO 101 K=1,L1 CR2 = CC(1,K,2)+CC(1,K,3) CH(1,1,K) = CC(1,K,1)+CR2 CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2)) CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 IF((IDO-1)/2.LT.L1) GO TO 104 DO 103 K=1,L1 CDIR$ IVDEP DO 102 I=3,IDO,2 IC = IDP2-I DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) CR2 = DR2+DR3 CI2 = DI2+DI3 CH(I-1,1,K) = CC(I-1,K,1)+CR2 CH(I,1,K) = CC(I,K,1)+CI2 TR2 = CC(I-1,K,1)+TAUR*CR2 TI2 = CC(I,K,1)+TAUR*CI2 TR3 = TAUI*(DI2-DI3) TI3 = TAUI*(DR3-DR2) CH(I-1,3,K) = TR2+TR3 CH(IC-1,2,K) = TR2-TR3 CH(I,3,K) = TI2+TI3 CH(IC,2,K) = TI3-TI2 102 CONTINUE 103 CONTINUE RETURN 104 DO 106 I=3,IDO,2 IC = IDP2-I CDIR$ IVDEP DO 105 K=1,L1 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) CR2 = DR2+DR3 CI2 = DI2+DI3 CH(I-1,1,K) = CC(I-1,K,1)+CR2 CH(I,1,K) = CC(I,K,1)+CI2 TR2 = CC(I-1,K,1)+TAUR*CR2 TI2 = CC(I,K,1)+TAUR*CI2 TR3 = TAUI*(DI2-DI3) TI3 = TAUI*(DR3-DR2) CH(I-1,3,K) = TR2+TR3 CH(IC-1,2,K) = TR2-TR3 CH(I,3,K) = TI2+TI3 CH(IC,2,K) = TI3-TI2 105 CONTINUE 106 CONTINUE RETURN END SUBROUTINE DRADF4(IDO,L1,CC,CH,WA1,WA2,WA3) C***BEGIN PROLOGUE DRADF4 C***REFER TO DRFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADF4 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) , 1 WA1(*) ,WA2(*) ,WA3(*) C***FIRST EXECUTABLE STATEMENT DRADF4 HSQT2 = 0.5D0*SQRT(2.0D0) DO 101 K=1,L1 TR1 = CC(1,K,2)+CC(1,K,4) TR2 = CC(1,K,1)+CC(1,K,3) CH(1,1,K) = TR1+TR2 CH(IDO,4,K) = TR2-TR1 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3) CH(1,3,K) = CC(1,K,4)-CC(1,K,2) 101 CONTINUE IF (IDO-2) 107,105,102 102 IDP2 = IDO+2 IF((IDO-1)/2.LT.L1) GO TO 111 DO 104 K=1,L1 CDIR$ IVDEP DO 103 I=3,IDO,2 IC = IDP2-I CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) TR1 = CR2+CR4 TR4 = CR4-CR2 TI1 = CI2+CI4 TI4 = CI2-CI4 TI2 = CC(I,K,1)+CI3 TI3 = CC(I,K,1)-CI3 TR2 = CC(I-1,K,1)+CR3 TR3 = CC(I-1,K,1)-CR3 CH(I-1,1,K) = TR1+TR2 CH(IC-1,4,K) = TR2-TR1 CH(I,1,K) = TI1+TI2 CH(IC,4,K) = TI1-TI2 CH(I-1,3,K) = TI4+TR3 CH(IC-1,2,K) = TR3-TI4 CH(I,3,K) = TR4+TI3 CH(IC,2,K) = TR4-TI3 103 CONTINUE 104 CONTINUE GO TO 110 111 DO 109 I=3,IDO,2 IC = IDP2-I CDIR$ IVDEP DO 108 K=1,L1 CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) TR1 = CR2+CR4 TR4 = CR4-CR2 TI1 = CI2+CI4 TI4 = CI2-CI4 TI2 = CC(I,K,1)+CI3 TI3 = CC(I,K,1)-CI3 TR2 = CC(I-1,K,1)+CR3 TR3 = CC(I-1,K,1)-CR3 CH(I-1,1,K) = TR1+TR2 CH(IC-1,4,K) = TR2-TR1 CH(I,1,K) = TI1+TI2 CH(IC,4,K) = TI1-TI2 CH(I-1,3,K) = TI4+TR3 CH(IC-1,2,K) = TR3-TI4 CH(I,3,K) = TR4+TI3 CH(IC,2,K) = TR4-TI3 108 CONTINUE 109 CONTINUE 110 IF (MOD(IDO,2) .EQ. 1) RETURN 105 DO 106 K=1,L1 TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4)) TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4)) CH(IDO,1,K) = TR1+CC(IDO,K,1) CH(IDO,3,K) = CC(IDO,K,1)-TR1 CH(1,2,K) = TI1-CC(IDO,K,3) CH(1,4,K) = TI1+CC(IDO,K,3) 106 CONTINUE 107 RETURN END SUBROUTINE DRADF5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4) C***BEGIN PROLOGUE DRADF5 C***REFER TO DRFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADF5 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) , 1 WA1(*) ,WA2(*) ,WA3(*) ,WA4(*) C***FIRST EXECUTABLE STATEMENT DRADF5 PI = 4.0D0*ATAN(1.0D0) TR11 = SIN(0.1D0*PI) TI11 = SIN(0.4D0*PI) TR12 = -SIN(0.3D0*PI) TI12 = SIN(0.2D0*PI) DO 101 K=1,L1 CR2 = CC(1,K,5)+CC(1,K,2) CI5 = CC(1,K,5)-CC(1,K,2) CR3 = CC(1,K,4)+CC(1,K,3) CI4 = CC(1,K,4)-CC(1,K,3) CH(1,1,K) = CC(1,K,1)+CR2+CR3 CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3 CH(1,3,K) = TI11*CI5+TI12*CI4 CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3 CH(1,5,K) = TI12*CI5-TI11*CI4 101 CONTINUE IF (IDO .EQ. 1) RETURN IDP2 = IDO+2 IF((IDO-1)/2.LT.L1) GO TO 104 DO 103 K=1,L1 CDIR$ IVDEP DO 102 I=3,IDO,2 IC = IDP2-I DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5) DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5) CR2 = DR2+DR5 CI5 = DR5-DR2 CR5 = DI2-DI5 CI2 = DI2+DI5 CR3 = DR3+DR4 CI4 = DR4-DR3 CR4 = DI3-DI4 CI3 = DI3+DI4 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3 CH(I,1,K) = CC(I,K,1)+CI2+CI3 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3 TR5 = TI11*CR5+TI12*CR4 TI5 = TI11*CI5+TI12*CI4 TR4 = TI12*CR5-TI11*CR4 TI4 = TI12*CI5-TI11*CI4 CH(I-1,3,K) = TR2+TR5 CH(IC-1,2,K) = TR2-TR5 CH(I,3,K) = TI2+TI5 CH(IC,2,K) = TI5-TI2 CH(I-1,5,K) = TR3+TR4 CH(IC-1,4,K) = TR3-TR4 CH(I,5,K) = TI3+TI4 CH(IC,4,K) = TI4-TI3 102 CONTINUE 103 CONTINUE RETURN 104 DO 106 I=3,IDO,2 IC = IDP2-I CDIR$ IVDEP DO 105 K=1,L1 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2) DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2) DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3) DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3) DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4) DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4) DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5) DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5) CR2 = DR2+DR5 CI5 = DR5-DR2 CR5 = DI2-DI5 CI2 = DI2+DI5 CR3 = DR3+DR4 CI4 = DR4-DR3 CR4 = DI3-DI4 CI3 = DI3+DI4 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3 CH(I,1,K) = CC(I,K,1)+CI2+CI3 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3 TR5 = TI11*CR5+TI12*CR4 TI5 = TI11*CI5+TI12*CI4 TR4 = TI12*CR5-TI11*CR4 TI4 = TI12*CI5-TI11*CI4 CH(I-1,3,K) = TR2+TR5 CH(IC-1,2,K) = TR2-TR5 CH(I,3,K) = TI2+TI5 CH(IC,2,K) = TI5-TI2 CH(I-1,5,K) = TR3+TR4 CH(IC-1,4,K) = TR3-TR4 CH(I,5,K) = TI3+TI4 CH(IC,4,K) = TI4-TI3 105 CONTINUE 106 CONTINUE RETURN END SUBROUTINE DRADFG(IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA) C***BEGIN PROLOGUE DRADFG C***REFER TO DRFFTF C***ROUTINES CALLED (NONE) C***END PROLOGUE DRADFG IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) , 1 C1(IDO,L1,IP) ,C2(IDL1,IP), 2 CH2(IDL1,IP) ,WA(*) C***FIRST EXECUTABLE STATEMENT DRADFG TPI = 8.0D0*ATAN(1.0D0) ARG = TPI/DBLE(IP) DCP = COS(ARG) DSP = SIN(ARG) IPPH = (IP+1)/2 IPP2 = IP+2 IDP2 = IDO+2 NBD = (IDO-1)/2 IF (IDO .EQ. 1) GO TO 119 DO 101 IK=1,IDL1 CH2(IK,1) = C2(IK,1) 101 CONTINUE DO 103 J=2,IP DO 102 K=1,L1 CH(1,K,J) = C1(1,K,J) 102 CONTINUE 103 CONTINUE IF (NBD .GT. L1) GO TO 107 IS = -IDO DO 106 J=2,IP IS = IS+IDO IDIJ = IS DO 105 I=3,IDO,2 IDIJ = IDIJ+2 DO 104 K=1,L1 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J) CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J) 104 CONTINUE 105 CONTINUE 106 CONTINUE GO TO 111 107 IS = -IDO DO 110 J=2,IP IS = IS+IDO DO 109 K=1,L1 IDIJ = IS CDIR$ IVDEP DO 108 I=3,IDO,2 IDIJ = IDIJ+2 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J) CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J) 108 CONTINUE 109 CONTINUE 110 CONTINUE 111 IF (NBD .LT. L1) GO TO 115 DO 114 J=2,IPPH JC = IPP2-J DO 113 K=1,L1 CDIR$ IVDEP DO 112 I=3,IDO,2 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC) C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC) C1(I,K,J) = CH(I,K,J)+CH(I,K,JC) C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J) 112 CONTINUE 113 CONTINUE 114 CONTINUE GO TO 121 115 DO 118 J=2,IPPH JC = IPP2-J DO 117 I=3,IDO,2 DO 116 K=1,L1 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC) C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC) C1(I,K,J) = CH(I,K,J)+CH(I,K,JC) C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J) 116 CONTINUE 117 CONTINUE 118 CONTINUE GO TO 121 119 DO 120 IK=1,IDL1 C2(IK,1) = CH2(IK,1) 120 CONTINUE 121 DO 123 J=2,IPPH JC = IPP2-J DO 122 K=1,L1 C1(1,K,J) = CH(1,K,J)+CH(1,K,JC) C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J) 122 CONTINUE 123 CONTINUE C AR1 = 1.0D0 AI1 = 0.0D0 DO 127 L=2,IPPH LC = IPP2-L AR1H = DCP*AR1-DSP*AI1 AI1 = DCP*AI1+DSP*AR1 AR1 = AR1H DO 124 IK=1,IDL1 CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2) CH2(IK,LC) = AI1*C2(IK,IP) 124 CONTINUE DC2 = AR1 DS2 = AI1 AR2 = AR1 AI2 = AI1 DO 126 J=3,IPPH JC = IPP2-J AR2H = DC2*AR2-DS2*AI2 AI2 = DC2*AI2+DS2*AR2 AR2 = AR2H DO 125 IK=1,IDL1 CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J) CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC) 125 CONTINUE 126 CONTINUE 127 CONTINUE DO 129 J=2,IPPH DO 128 IK=1,IDL1 CH2(IK,1) = CH2(IK,1)+C2(IK,J) 128 CONTINUE 129 CONTINUE C IF (IDO .LT. L1) GO TO 132 DO 131 K=1,L1 DO 130 I=1,IDO CC(I,1,K) = CH(I,K,1) 130 CONTINUE 131 CONTINUE GO TO 135 132 DO 134 I=1,IDO DO 133 K=1,L1 CC(I,1,K) = CH(I,K,1) 133 CONTINUE 134 CONTINUE 135 DO 137 J=2,IPPH JC = IPP2-J J2 = J+J DO 136 K=1,L1 CC(IDO,J2-2,K) = CH(1,K,J) CC(1,J2-1,K) = CH(1,K,JC) 136 CONTINUE 137 CONTINUE IF (IDO .EQ. 1) RETURN IF (NBD .LT. L1) GO TO 141 DO 140 J=2,IPPH JC = IPP2-J J2 = J+J DO 139 K=1,L1 CDIR$ IVDEP DO 138 I=3,IDO,2 IC = IDP2-I CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC) CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC) CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC) CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J) 138 CONTINUE 139 CONTINUE 140 CONTINUE RETURN 141 DO 144 J=2,IPPH JC = IPP2-J J2 = J+J DO 143 I=3,IDO,2 IC = IDP2-I DO 142 K=1,L1 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC) CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC) CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC) CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J) 142 CONTINUE 143 CONTINUE 144 CONTINUE RETURN END SUBROUTINE DEZFTI(N,WSAVE) C***BEGIN PROLOGUE DEZFTI C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 860115 (YYMMDD) C***CATEGORY NO. J1A1 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE Initialize DEZFTF and DEZFTB C***DESCRIPTION C C Subroutine DEZFTI initializes the array WSAVE which is used in C both DEZFTF and DEZFTB. The prime factorization of N together with C a tabulation of the trigonometric functions are computed and C stored in WSAVE. C C Input Parameter C C N the length of the sequence to be transformed. C C Output Parameter C C WSAVE a work array which must be dimensioned at least 3*N+15. C The same work array can be used for both DEZFTF and DEZFTB C as long as N remains unchanged. Different WSAVE arrays C are required for different values of N. C***REFERENCES (NONE) C***ROUTINES CALLED DEZFT1 C***END PROLOGUE DEZFTI IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WSAVE(*) C***FIRST EXECUTABLE STATEMENT DEZFTI IF (N .EQ. 1) RETURN CALL DEZFT1 (N,WSAVE(2*N+1),WSAVE(3*N+1)) RETURN END SUBROUTINE DEZFT1(N,WA,IFAC) C***BEGIN PROLOGUE DEZFT1 C***REFER TO DEZFTI C***ROUTINES CALLED (NONE) C***END PROLOGUE DEZFT1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ C***FIRST EXECUTABLE STATEMENT DEZFT1 TPI = 8.0D0*ATAN(1.0D0) NL = N NF = 0 J = 0 101 J = J+1 IF (J-4) 102,102,103 102 NTRY = NTRYH(J) GO TO 104 103 NTRY = NTRY+2 104 NQ = NL/NTRY NR = NL-NTRY*NQ IF (NR) 101,105,101 105 NF = NF+1 IFAC(NF+2) = NTRY NL = NQ IF (NTRY .NE. 2) GO TO 107 IF (NF .EQ. 1) GO TO 107 DO 106 I=2,NF IB = NF-I+2 IFAC(IB+2) = IFAC(IB+1) 106 CONTINUE IFAC(3) = 2 107 IF (NL .NE. 1) GO TO 104 IFAC(1) = N IFAC(2) = NF ARGH = TPI/DBLE(N) IS = 0 NFM1 = NF-1 L1 = 1 IF (NFM1 .EQ. 0) RETURN DO 111 K1=1,NFM1 IP = IFAC(K1+2) L2 = L1*IP IDO = N/L2 IPM = IP-1 ARG1 = DBLE(L1)*ARGH CH1 = 1.0D0 SH1 = 0.0D0 DCH1 = COS(ARG1) DSH1 = SIN(ARG1) DO 110 J=1,IPM CH1H = DCH1*CH1-DSH1*SH1 SH1 = DCH1*SH1+DSH1*CH1 CH1 = CH1H I = IS+2 WA(I-1) = CH1 WA(I) = SH1 IF (IDO .LT. 5) GO TO 109 DO 108 II=5,IDO,2 I = I+2 WA(I-1) = CH1*WA(I-3)-SH1*WA(I-2) WA(I) = CH1*WA(I-2)+SH1*WA(I-3) 108 CONTINUE 109 IS = IS+IDO 110 CONTINUE L1 = L2 111 CONTINUE RETURN END SUBROUTINE DRFFTI(N,WSAVE) C***BEGIN PROLOGUE DRFFTI C***DATE WRITTEN 790601 (YYMMDD) C***REVISION DATE 860115 (YYMMDD) C***CATEGORY NO. J1A1 C***KEYWORDS FOURIER TRANSFORM C***AUTHOR SWARZTRAUBER, P. N., (NCAR) C***PURPOSE Initialize for DRFFTF and DRFFTB. C***DESCRIPTION C C From the book, "Numerical Methods and Software" by C D. Kahaner, C. Moler, S. Nash C Prentice Hall, 1989 C C Subroutine DRFFTI initializes the array WSAVE which is used in C both DRFFTF and DRFFTB. The prime factorization of N together with C a tabulation of the trigonometric functions are computed and C stored in WSAVE. C C Input Parameter C C N the length of the sequence to be transformed. C C Output Parameter C C WSAVE a DP work array which must be dimensioned at least 2*N+15. C The same work array can be used for both DRFFTF and DRFFTB C as long as N remains unchanged. Different WSAVE arrays C are required for different values of N. The contents of C WSAVE must not be changed between calls of DRFFTF or DRFFTB. C C ********************************************************************* C * * C * SUBPROGRAM REVISION HISTORY * C * * C * 06/01/79 - Original version by Paul Swarztrauber. * C * Distributed by NCAR (ref. 1). * C * 04/01/83 - SLATEC Common Math Library Subcommittee. * C * Modified to use SLATEC library source file format. * C * Distributed in the SLATEC library (ref. 2). * C * 01/15/86 - Ron Boisvert, National Bureau of Standards. * C * Modified to convert to portable Fortran 77. * C * * C * The changes introduced in the most recent modification are * C * * C * (a) Dummy array size declarations (1) changed to (*) * C * (b) References to intrinsic function FLOAT changed to DP * C * (c) Mathematical constants previously coded in DATA state- * C * ments now computed at runtime using Fortran intrinsic * C * functions. The affected variables are * C * * C * PI SQRT2 SQRT3 TAUR TR11 TR12 * C * PIH TSQRT2 TAUI TI11 TI12 * C * TPI HSQT2 * C * * C * References * C * * C * 1. P.N. Swarztrauber, Vectorizing the FFTs, in Parallel * C * Computations (G. Rodrigue, ed.), Academic Press, 1982, * C * pp. 51-83. * C * 2. B.L. Buzbee, The SLATEC Common Math Library, in Sources * C * and Development of Mathematical Software (W. Cowell, ed.), * C * Prentice-Hall, 1984, pp. 302-318. * C * * C ********************************************************************* C C***REFERENCES (NONE) C***ROUTINES CALLED DRFTI1 C***END PROLOGUE DRFFTI IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WSAVE(*) C***FIRST EXECUTABLE STATEMENT DRFFTI IF (N .EQ. 1) RETURN CALL DRFTI1 (N,WSAVE(N+1),WSAVE(2*N+1)) RETURN END SUBROUTINE DRFTI1(N,WA,IFAC) C***BEGIN PROLOGUE DRFTI1 C***REFER TO DRFFTI C***ROUTINES CALLED (NONE) C***END PROLOGUE DRFTI1 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION WA(*) ,IFAC(*) ,NTRYH(4) DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/ C***FIRST EXECUTABLE STATEMENT DRFTI1 NL = N NF = 0 J = 0 101 J = J+1 IF (J-4) 102,102,103 102 NTRY = NTRYH(J) GO TO 104 103 NTRY = NTRY+2 104 NQ = NL/NTRY NR = NL-NTRY*NQ IF (NR) 101,105,101 105 NF = NF+1 IFAC(NF+2) = NTRY NL = NQ IF (NTRY .NE. 2) GO TO 107 IF (NF .EQ. 1) GO TO 107 DO 106 I=2,NF IB = NF-I+2 IFAC(IB+2) = IFAC(IB+1) 106 CONTINUE IFAC(3) = 2 107 IF (NL .NE. 1) GO TO 104 IFAC(1) = N IFAC(2) = NF TPI = 8.D0*ATAN(1.D0) ARGH = TPI/N IS = 0 NFM1 = NF-1 L1 = 1 IF (NFM1 .EQ. 0) RETURN DO 110 K1=1,NFM1 IP = IFAC(K1+2) LD = 0 L2 = L1*IP IDO = N/L2 IPM = IP-1 DO 109 J=1,IPM LD = LD+L1 I = IS ARGLD = LD*ARGH FI = 0.D0 DO 108 II=3,IDO,2 I = I+2 FI = FI+1.D0 ARG = FI*ARGLD WA(I-1) = COS(ARG) WA(I) = SIN(ARG) 108 CONTINUE IS = IS+IDO 109 CONTINUE L1 = L2 110 CONTINUE RETURN END