C C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C * * C * F F T P A C K * C * * C * * C * A PACKAGE OF FORTRAN SUBPROGRAMS FOR THE SOLUTION OF * C * * C * SEPARABLE ELLIPTIC PARTIAL DIFFERENTIAL EQUATIONS * C * * C * (VERSION 4.1 , NOVEMBER 1988) * C * * C * BY * C * * C * PAUL SWARZTRAUBER * C * * C * OF * C * * C * THE NATIONAL CENTER FOR ATMOSPHERIC RESEARCH * C * * C * BOULDER, COLORADO (80307) U.S.A. * C * * C * WHICH IS SPONSORED BY * C * * C * THE NATIONAL SCIENCE FOUNDATION * C * * C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * C C C LATEST REVISION C --------------- C NOVEMBER 1988 (VERSION 4.1) C C PURPOSE C ------- C THIS PACKAGE CONSISTS OF PROGRAMS WHICH PERFORM FAST FOURIER C TRANSFORMS FOR BOTH COMPLEX AND REAL PERIODIC SEQUENCES AND C CERTAIN OTHER SYMMETRIC SEQUENCES THAT ARE LISTED BELOW. C C USAGE C ----- C 1. RFFTI INITIALIZE RFFTF AND RFFTB C 2. RFFTF FORWARD TRANSFORM OF A REAL PERIODIC SEQUENCE C 3. RFFTB BACKWARD TRANSFORM OF A REAL COEFFICIENT ARRAY C C 4. EZFFTI INITIALIZE EZFFTF AND EZFFTB C 5. EZFFTF A SIMPLIFIED REAL PERIODIC FORWARD TRANSFORM C 6. EZFFTB A SIMPLIFIED REAL PERIODIC BACKWARD TRANSFORM C C 7. SINTI INITIALIZE SINT C 8. SINT SINE TRANSFORM OF A REAL ODD SEQUENCE C C 9. COSTI INITIALIZE COST C 10. COST COSINE TRANSFORM OF A REAL EVEN SEQUENCE C C 11. SINQI INITIALIZE SINQF AND SINQB C 12. SINQF FORWARD SINE TRANSFORM WITH ODD WAVE NUMBERS C 13. SINQB UNNORMALIZED INVERSE OF SINQF C C 14. COSQI INITIALIZE COSQF AND COSQB C 15. COSQF FORWARD COSINE TRANSFORM WITH ODD WAVE NUMBERS C 16. COSQB UNNORMALIZED INVERSE OF COSQF C C 17. CFFTI INITIALIZE CFFTF AND CFFTB C 18. CFFTF FORWARD TRANSFORM OF A COMPLEX PERIODIC SEQUENCE C 19. CFFTB UNNORMALIZED INVERSE OF CFFTF C C SPECIAL CONDITIONS C ------------------ C BEFORE CALLING ROUTINES EZFFTB AND EZFFTF FOR THE FIRST TIME, C OR BEFORE CALLING EZFFTB AND EZFFTF WITH A DIFFERENT LENGTH, C USERS MUST INITIALIZE BY CALLING ROUTINE EZFFTI. C C I/O C --- C NONE C C PRECISION C --------- C NONE C C REQUIRED LIBRARY FILES C ---------------------- C NONE C C LANGUAGE C -------- C FORTRAN C C HISTORY C ------- C DEVELOPED AT NCAR IN BOULDER, COLORADO BY PAUL N. SWARZTRAUBER C OF THE SCIENTIFIC COMPUTING DIVISION. RELEASED ON NCAR'S PUBLIC C SOFTWARE LIBRARIES IN JANUARY 1980. MODIFIED MAY 29, 1985 TO C INCREASE EFFICIENCY. C C PORTABILITY C ----------- C FORTRAN 77 C C ********************************************************************** C C SUBROUTINE RFFTI(N,WSAVE) C C SUBROUTINE RFFTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN C BOTH RFFTF AND RFFTB. 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 2*N+15. C THE SAME WORK ARRAY CAN BE USED FOR BOTH RFFTF AND RFFTB 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 RFFTF OR RFFTB. C C ********************************************************************** C C SUBROUTINE RFFTF(N,R,WSAVE) C C SUBROUTINE RFFTF COMPUTES THE FOURIER COEFFICIENTS OF A REAL 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 REAL ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE C TO BE TRANSFORMED C C WSAVE A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 2*N+15. C IN THE PROGRAM THAT CALLS RFFTF. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE RFFTI(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 RFFTF AND RFFTB. 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 RFFTF C FOLLOWED BY A CALL OF RFFTB WILL MULTIPLY THE INPUT C SEQUENCE BY N. C C WSAVE CONTAINS RESULTS WHICH MUST NOT BE DESTROYED BETWEEN C CALLS OF RFFTF OR RFFTB. C C C ********************************************************************** C C SUBROUTINE RFFTB(N,R,WSAVE) C C SUBROUTINE RFFTB COMPUTES THE REAL 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 REAL ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE C TO BE TRANSFORMED C C WSAVE A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 2*N+15. C IN THE PROGRAM THAT CALLS RFFTB. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE RFFTI(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 RFFTF AND RFFTB. 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 RFFTF C FOLLOWED BY A CALL OF RFFTB WILL MULTIPLY THE INPUT C SEQUENCE BY N. C C WSAVE CONTAINS RESULTS WHICH MUST NOT BE DESTROYED BETWEEN C CALLS OF RFFTB OR RFFTF. C C C ********************************************************************** C C SUBROUTINE EZFFTI(N,WSAVE) C C SUBROUTINE EZFFTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN C BOTH EZFFTF AND EZFFTB. 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 EZFFTF AND EZFFTB C AS LONG AS N REMAINS UNCHANGED. DIFFERENT WSAVE ARRAYS C ARE REQUIRED FOR DIFFERENT VALUES OF N. C C C ********************************************************************** C C SUBROUTINE EZFFTF(N,R,AZERO,A,B,WSAVE) C C SUBROUTINE EZFFTF COMPUTES THE FOURIER COEFFICIENTS OF A REAL C PERODIC SEQUENCE (FOURIER ANALYSIS). THE TRANSFORM IS DEFINED C BELOW AT OUTPUT PARAMETERS AZERO,A AND B. EZFFTF IS A SIMPLIFIED C BUT SLOWER VERSION OF RFFTF. 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 REAL 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 EZFFTF. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE EZFFTI(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 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. 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./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./N*R(I)*SIN(K*(I-1)*2*PI/N) C C C ********************************************************************** C C SUBROUTINE EZFFTB(N,R,AZERO,A,B,WSAVE) C C SUBROUTINE EZFFTB COMPUTES A REAL PERODIC SEQUENCE FROM ITS C FOURIER COEFFICIENTS (FOURIER SYNTHESIS). THE TRANSFORM IS C DEFINED BELOW AT OUTPUT PARAMETER R. EZFFTB IS A SIMPLIFIED C BUT SLOWER VERSION OF RFFTB. 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 EZFFTB. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE EZFFTI(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 EZFFTB. 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) = .5*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 C SUBROUTINE SINTI(N,WSAVE) C C SUBROUTINE SINTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN C SUBROUTINE SINT. 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. THE METHOD C IS MOST EFFICIENT WHEN N+1 IS A PRODUCT OF SMALL PRIMES. C C OUTPUT PARAMETER C C WSAVE A WORK ARRAY WITH AT LEAST INT(2.5*N+15) LOCATIONS. C DIFFERENT WSAVE ARRAYS ARE REQUIRED FOR DIFFERENT VALUES C OF N. THE CONTENTS OF WSAVE MUST NOT BE CHANGED BETWEEN C CALLS OF SINT. C C ********************************************************************** C C SUBROUTINE SINT(N,X,WSAVE) C C SUBROUTINE SINT COMPUTES THE DISCRETE FOURIER SINE TRANSFORM C OF AN ODD SEQUENCE X(I). THE TRANSFORM IS DEFINED BELOW AT C OUTPUT PARAMETER X. C C SINT IS THE UNNORMALIZED INVERSE OF ITSELF SINCE A CALL OF SINT C FOLLOWED BY ANOTHER CALL OF SINT WILL MULTIPLY THE INPUT SEQUENCE C X BY 2*(N+1). C C THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE SINT MUST BE C INITIALIZED BY CALLING SUBROUTINE SINTI(N,WSAVE). C C INPUT PARAMETERS C C N THE LENGTH OF THE SEQUENCE TO BE TRANSFORMED. THE METHOD C IS MOST EFFICIENT WHEN N+1 IS THE PRODUCT OF SMALL PRIMES. C C X AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED C C C WSAVE A WORK ARRAY WITH DIMENSION AT LEAST INT(2.5*N+15) C IN THE PROGRAM THAT CALLS SINT. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE SINTI(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 C OUTPUT PARAMETERS C C X FOR I=1,...,N C C X(I)= THE SUM FROM K=1 TO K=N C C 2*X(K)*SIN(K*I*PI/(N+1)) C C A CALL OF SINT FOLLOWED BY ANOTHER CALL OF C SINT WILL MULTIPLY THE SEQUENCE X BY 2*(N+1). C HENCE SINT IS THE UNNORMALIZED INVERSE C OF ITSELF. C C WSAVE CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT BE C DESTROYED BETWEEN CALLS OF SINT. C C ********************************************************************** C C SUBROUTINE COSTI(N,WSAVE) C C SUBROUTINE COSTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN C SUBROUTINE COST. 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. THE METHOD C IS MOST EFFICIENT WHEN N-1 IS A PRODUCT OF SMALL PRIMES. C C OUTPUT PARAMETER C C WSAVE A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3*N+15. C DIFFERENT WSAVE ARRAYS ARE REQUIRED FOR DIFFERENT VALUES C OF N. THE CONTENTS OF WSAVE MUST NOT BE CHANGED BETWEEN C CALLS OF COST. C C ********************************************************************** C C SUBROUTINE COST(N,X,WSAVE) C C SUBROUTINE COST COMPUTES THE DISCRETE FOURIER COSINE TRANSFORM C OF AN EVEN SEQUENCE X(I). THE TRANSFORM IS DEFINED BELOW AT OUTPUT C PARAMETER X. C C COST IS THE UNNORMALIZED INVERSE OF ITSELF SINCE A CALL OF COST C FOLLOWED BY ANOTHER CALL OF COST WILL MULTIPLY THE INPUT SEQUENCE C X BY 2*(N-1). THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER X C C THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE COST MUST BE C INITIALIZED BY CALLING SUBROUTINE COSTI(N,WSAVE). C C INPUT PARAMETERS C C N THE LENGTH OF THE SEQUENCE X. N MUST BE GREATER THAN 1. C THE METHOD IS MOST EFFICIENT WHEN N-1 IS A PRODUCT OF C SMALL PRIMES. C C X AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED C C WSAVE A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3*N+15 C IN THE PROGRAM THAT CALLS COST. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE COSTI(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 C OUTPUT PARAMETERS C C X FOR I=1,...,N C C X(I) = X(1)+(-1)**(I-1)*X(N) C C + THE SUM FROM K=2 TO K=N-1 C C 2*X(K)*COS((K-1)*(I-1)*PI/(N-1)) C C A CALL OF COST FOLLOWED BY ANOTHER CALL OF C COST WILL MULTIPLY THE SEQUENCE X BY 2*(N-1) C HENCE COST IS THE UNNORMALIZED INVERSE C OF ITSELF. C C WSAVE CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT BE C DESTROYED BETWEEN CALLS OF COST. C C ********************************************************************** C C SUBROUTINE SINQI(N,WSAVE) C C SUBROUTINE SINQI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN C BOTH SINQF AND SINQB. 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. THE METHOD C IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES. 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 SINQF AND SINQB 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 SINQF OR SINQB. C C ********************************************************************** C C SUBROUTINE SINQF(N,X,WSAVE) C C SUBROUTINE SINQF COMPUTES THE FAST FOURIER TRANSFORM OF QUARTER C WAVE DATA. THAT IS , SINQF COMPUTES THE COEFFICIENTS IN A SINE C SERIES REPRESENTATION WITH ONLY ODD WAVE NUMBERS. THE TRANSFORM C IS DEFINED BELOW AT OUTPUT PARAMETER X. C C SINQB IS THE UNNORMALIZED INVERSE OF SINQF SINCE A CALL OF SINQF C FOLLOWED BY A CALL OF SINQB WILL MULTIPLY THE INPUT SEQUENCE X C BY 4*N. C C THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE SINQF MUST BE C INITIALIZED BY CALLING SUBROUTINE SINQI(N,WSAVE). C C C INPUT PARAMETERS C C N THE LENGTH OF THE ARRAY X TO BE TRANSFORMED. THE METHOD C IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES. C C X AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED C C WSAVE A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3*N+15. C IN THE PROGRAM THAT CALLS SINQF. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE SINQI(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 C OUTPUT PARAMETERS C C X FOR I=1,...,N C C X(I) = (-1)**(I-1)*X(N) C C + THE SUM FROM K=1 TO K=N-1 OF C C 2*X(K)*SIN((2*I-1)*K*PI/(2*N)) C C A CALL OF SINQF FOLLOWED BY A CALL OF C SINQB WILL MULTIPLY THE SEQUENCE X BY 4*N. C THEREFORE SINQB IS THE UNNORMALIZED INVERSE C OF SINQF. C C WSAVE CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT C BE DESTROYED BETWEEN CALLS OF SINQF OR SINQB. C C ********************************************************************** C C SUBROUTINE SINQB(N,X,WSAVE) C C SUBROUTINE SINQB COMPUTES THE FAST FOURIER TRANSFORM OF QUARTER C WAVE DATA. THAT IS , SINQB COMPUTES A SEQUENCE FROM ITS C REPRESENTATION IN TERMS OF A SINE SERIES WITH ODD WAVE NUMBERS. C THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER X. C C SINQF IS THE UNNORMALIZED INVERSE OF SINQB SINCE A CALL OF SINQB C FOLLOWED BY A CALL OF SINQF WILL MULTIPLY THE INPUT SEQUENCE X C BY 4*N. C C THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE SINQB MUST BE C INITIALIZED BY CALLING SUBROUTINE SINQI(N,WSAVE). C C C INPUT PARAMETERS C C N THE LENGTH OF THE ARRAY X TO BE TRANSFORMED. THE METHOD C IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES. C C X AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED C C WSAVE A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3*N+15. C IN THE PROGRAM THAT CALLS SINQB. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE SINQI(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 C OUTPUT PARAMETERS C C X FOR I=1,...,N C C X(I)= THE SUM FROM K=1 TO K=N OF C C 4*X(K)*SIN((2K-1)*I*PI/(2*N)) C C A CALL OF SINQB FOLLOWED BY A CALL OF C SINQF WILL MULTIPLY THE SEQUENCE X BY 4*N. C THEREFORE SINQF IS THE UNNORMALIZED INVERSE C OF SINQB. C C WSAVE CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT C BE DESTROYED BETWEEN CALLS OF SINQB OR SINQF. C C ********************************************************************** C C SUBROUTINE COSQI(N,WSAVE) C C SUBROUTINE COSQI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN C BOTH COSQF AND COSQB. 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 ARRAY TO BE TRANSFORMED. THE METHOD C IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES. 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 COSQF AND COSQB 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 COSQF OR COSQB. C C ********************************************************************** C C SUBROUTINE COSQF(N,X,WSAVE) C C SUBROUTINE COSQF COMPUTES THE FAST FOURIER TRANSFORM OF QUARTER C WAVE DATA. THAT IS , COSQF COMPUTES THE COEFFICIENTS IN A COSINE C SERIES REPRESENTATION WITH ONLY ODD WAVE NUMBERS. THE TRANSFORM C IS DEFINED BELOW AT OUTPUT PARAMETER X C C COSQF IS THE UNNORMALIZED INVERSE OF COSQB SINCE A CALL OF COSQF C FOLLOWED BY A CALL OF COSQB WILL MULTIPLY THE INPUT SEQUENCE X C BY 4*N. C C THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE COSQF MUST BE C INITIALIZED BY CALLING SUBROUTINE COSQI(N,WSAVE). C C C INPUT PARAMETERS C C N THE LENGTH OF THE ARRAY X TO BE TRANSFORMED. THE METHOD C IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES. C C X AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED C C WSAVE A WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 3*N+15 C IN THE PROGRAM THAT CALLS COSQF. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE COSQI(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 C OUTPUT PARAMETERS C C X FOR I=1,...,N C C X(I) = X(1) PLUS THE SUM FROM K=2 TO K=N OF C C 2*X(K)*COS((2*I-1)*(K-1)*PI/(2*N)) C C A CALL OF COSQF FOLLOWED BY A CALL OF C COSQB WILL MULTIPLY THE SEQUENCE X BY 4*N. C THEREFORE COSQB IS THE UNNORMALIZED INVERSE C OF COSQF. C C WSAVE CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT C BE DESTROYED BETWEEN CALLS OF COSQF OR COSQB. C C ********************************************************************** C C SUBROUTINE COSQB(N,X,WSAVE) C C SUBROUTINE COSQB COMPUTES THE FAST FOURIER TRANSFORM OF QUARTER C WAVE DATA. THAT IS , COSQB COMPUTES A SEQUENCE FROM ITS C REPRESENTATION IN TERMS OF A COSINE SERIES WITH ODD WAVE NUMBERS. C THE TRANSFORM IS DEFINED BELOW AT OUTPUT PARAMETER X. C C COSQB IS THE UNNORMALIZED INVERSE OF COSQF SINCE A CALL OF COSQB C FOLLOWED BY A CALL OF COSQF WILL MULTIPLY THE INPUT SEQUENCE X C BY 4*N. C C THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE COSQB MUST BE C INITIALIZED BY CALLING SUBROUTINE COSQI(N,WSAVE). C C C INPUT PARAMETERS C C N THE LENGTH OF THE ARRAY X TO BE TRANSFORMED. THE METHOD C IS MOST EFFICIENT WHEN N IS A PRODUCT OF SMALL PRIMES. C C X AN ARRAY WHICH CONTAINS THE SEQUENCE TO BE TRANSFORMED C C WSAVE A WORK ARRAY THAT MUST BE DIMENSIONED AT LEAST 3*N+15 C IN THE PROGRAM THAT CALLS COSQB. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE COSQI(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 C OUTPUT PARAMETERS C C X FOR I=1,...,N C C X(I)= THE SUM FROM K=1 TO K=N OF C C 4*X(K)*COS((2*K-1)*(I-1)*PI/(2*N)) C C A CALL OF COSQB FOLLOWED BY A CALL OF C COSQF WILL MULTIPLY THE SEQUENCE X BY 4*N. C THEREFORE COSQF IS THE UNNORMALIZED INVERSE C OF COSQB. C C WSAVE CONTAINS INITIALIZATION CALCULATIONS WHICH MUST NOT C BE DESTROYED BETWEEN CALLS OF COSQB OR COSQF. C C ********************************************************************** C C SUBROUTINE CFFTI(N,WSAVE) C C SUBROUTINE CFFTI INITIALIZES THE ARRAY WSAVE WHICH IS USED IN C BOTH CFFTF AND CFFTB. 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 CFFTF AND CFFTB 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 CFFTF OR CFFTB. C C ********************************************************************** C C SUBROUTINE CFFTF(N,C,WSAVE) C C SUBROUTINE CFFTF COMPUTES THE FORWARD COMPLEX DISCRETE FOURIER C TRANSFORM (THE FOURIER ANALYSIS). EQUIVALENTLY , CFFTF 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 CFFTF C FOLLOWED BY A CALL OF CFFTB WILL MULTIPLY THE SEQUENCE BY N. C C THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE CFFTF MUST BE C INITIALIZED BY CALLING SUBROUTINE CFFTI(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. N C C C A COMPLEX ARRAY OF LENGTH N WHICH CONTAINS THE SEQUENCE C C WSAVE A REAL WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4N+15 C IN THE PROGRAM THAT CALLS CFFTF. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE CFFTI(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 CFFTF AND CFFTB. 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-1)*(K-1)*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 CFFTF OR CFFTB C C ********************************************************************** C C SUBROUTINE CFFTB(N,C,WSAVE) C C SUBROUTINE CFFTB COMPUTES THE BACKWARD COMPLEX DISCRETE FOURIER C TRANSFORM (THE FOURIER SYNTHESIS). EQUIVALENTLY , CFFTB 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 CFFTF FOLLOWED BY A CALL OF CFFTB WILL MULTIPLY THE C SEQUENCE BY N. C C THE ARRAY WSAVE WHICH IS USED BY SUBROUTINE CFFTB MUST BE C INITIALIZED BY CALLING SUBROUTINE CFFTI(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 REAL WORK ARRAY WHICH MUST BE DIMENSIONED AT LEAST 4N+15 C IN THE PROGRAM THAT CALLS CFFTB. THE WSAVE ARRAY MUST BE C INITIALIZED BY CALLING SUBROUTINE CFFTI(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 CFFTF AND CFFTB. 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-1)*(K-1)*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 CFFTF OR CFFTB C **********************************************************************