SUBROUTINE fft99(a,work,trigs,ifax,inc,jump,n,lot,ISIGN),4
!
! PURPOSE PERFORMS MULTIPLE FAST FOURIER TRANSFORMS. THIS PACKAGE
! WILL PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX
! PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE
! TRANSFORMS, I.E. GIVEN A SET OF REAL DATA VECTORS, THE
! PACKAGE RETURNS A SET OF 'HALF-COMPLEX' FOURIER
! COEFFICIENT VECTORS, OR VICE VERSA. THE LENGTH OF THE
! TRANSFORMS MUST BE AN EVEN NUMBER GREATER THAN 4 THAT HAS
! NO OTHER FACTORS EXCEPT POSSIBLY POWERS OF 2, 3, AND 5.
! THIS IS AN ALL FORTRAN VERSION OF THE CRAYLIB PACKAGE
! THAT IS MOSTLY WRITTEN IN CAL.
!
! THE PACKAGE FFT99F CONTAINS SEVERAL USER-LEVEL ROUTINES:
!
! SUBROUTINE SET99
! AN INITIALIZATION ROUTINE THAT MUST BE CALLED ONCE
! BEFORE A SEQUENCE OF CALLS TO THE FFT ROUTINES
! (PROVIDED THAT N IS NOT CHANGED).
!
! SUBROUTINES FFT99 AND FFT991
! TWO FFT ROUTINES THAT RETURN SLIGHTLY DIFFERENT
! ARRANGEMENTS OF THE DATA IN GRIDPOINT SPACE.
!
!
! ACCESS THIS FORTRAN VERSION MAY BE ACCESSED WITH
!
! *FORTRAN,P=XLIB,SN=FFT99F
!
! TO ACCESS THE CRAY OBJECT CODE, CALLING THE USER ENTRY
! POINTS FROM A CRAY PROGRAM IS SUFFICIENT. THE SOURCE
! FORTRAN AND CAL CODE FOR THE CRAYLIB VERSION MAY BE
! ACCESSED USING
!
! FETCH P=CRAYLIB,SN=FFT99
! FETCH P=CRAYLIB,SN=CAL99
!
! USAGE LET N BE OF THE FORM 2**P * 3**Q * 5**R, WHERE P .GE. 1,
! Q .GE. 0, AND R .GE. 0. THEN A TYPICAL SEQUENCE OF
! CALLS TO TRANSFORM A GIVEN SET OF REAL VECTORS OF LENGTH
! N TO A SET OF 'HALF-COMPLEX' FOURIER COEFFICIENT VECTORS
! OF LENGTH N IS
!
! DIMENSION IFAX(13),TRIGS(3*N/2+1),A(M*(N+2)),
! + WORK(M*(N+1))
!
! CALL SET99 (TRIGS, IFAX, N)
! CALL FFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN)
!
! SEE THE INDIVIDUAL WRITE-UPS FOR SET99, FFT99, AND
! FFT991 BELOW, FOR A DETAILED DESCRIPTION OF THE
! ARGUMENTS.
!
! HISTORY THE PACKAGE WAS WRITTEN BY CLIVE TEMPERTON AT ECMWF IN
! NOVEMBER, 1978. IT WAS MODIFIED, DOCUMENTED, AND TESTED
! FOR NCAR BY RUSS REW IN SEPTEMBER, 1980.
!
!-----------------------------------------------------------------------
!
! SUBROUTINE SET99 (TRIGS, IFAX, N)
!
! PURPOSE A SET-UP ROUTINE FOR FFT99 AND FFT991. IT NEED ONLY BE
! CALLED ONCE BEFORE A SEQUENCE OF CALLS TO THE FFT
! ROUTINES (PROVIDED THAT N IS NOT CHANGED).
!
! ARGUMENT IFAX(13),TRIGS(3*N/2+1)
! DIMENSIONS
!
! ARGUMENTS
!
! ON INPUT TRIGS
! A FLOATING POINT ARRAY OF DIMENSION 3*N/2 IF N/2 IS
! EVEN, OR 3*N/2+1 IF N/2 IS ODD.
!
! IFAX
! AN INTEGER ARRAY. THE NUMBER OF ELEMENTS ACTUALLY USED
! WILL DEPEND ON THE FACTORIZATION OF N. DIMENSIONING
! IFAX FOR 13 SUFFICES FOR ALL N LESS THAN A MILLION.
!
! N
! AN EVEN NUMBER GREATER THAN 4 THAT HAS NO PRIME FACTOR
! GREATER THAN 5. N IS THE LENGTH OF THE TRANSFORMS (SEE
! THE DOCUMENTATION FOR FFT99 AND FFT991 FOR THE
! DEFINITIONS OF THE TRANSFORMS).
!
! ON OUTPUT IFAX
! CONTAINS THE FACTORIZATION OF N/2. IFAX(1) IS THE
! NUMBER OF FACTORS, AND THE FACTORS THEMSELVES ARE STORED
! IN IFAX(2),IFAX(3),... IF SET99 IS CALLED WITH N ODD,
! OR IF N HAS ANY PRIME FACTORS GREATER THAN 5, IFAX(1)
! IS SET TO -99.
!
! TRIGS
! AN ARRAY OF TRIGNOMENTRIC FUNCTION VALUES SUBSEQUENTLY
! USED BY THE FFT ROUTINES.
!
!-----------------------------------------------------------------------
!
! SUBROUTINE FFT991 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN)
! AND
! SUBROUTINE FFT99 (A,WORK,TRIGS,IFAX,INC,JUMP,N,M,ISIGN)
!
! PURPOSE PERFORM A NUMBER OF SIMULTANEOUS REAL/HALF-COMPLEX
! PERIODIC FOURIER TRANSFORMS OR CORRESPONDING INVERSE
! TRANSFORMS, USING ORDINARY SPATIAL ORDER OF GRIDPOINT
! VALUES (FFT991) OR EXPLICIT CYCLIC CONTINUITY IN THE
! GRIDPOINT VALUES (FFT99). GIVEN A SET
! OF REAL DATA VECTORS, THE PACKAGE RETURNS A SET OF
! 'HALF-COMPLEX' FOURIER COEFFICIENT VECTORS, OR VICE
! VERSA. THE LENGTH OF THE TRANSFORMS MUST BE AN EVEN
! NUMBER THAT HAS NO OTHER FACTORS EXCEPT POSSIBLY POWERS
! OF 2, 3, AND 5. THESE VERSION OF FFT991 AND FFT99 ARE
! OPTIMIZED FOR USE ON THE CRAY-1.
!
! ARGUMENT A(M*(N+2)), WORK(M*(N+1)), TRIGS(3*N/2+1), IFAX(13)
! DIMENSIONS
!
! ARGUMENTS
!
! ON INPUT A
! AN ARRAY OF LENGTH M*(N+2) CONTAINING THE INPUT DATA
! OR COEFFICIENT VECTORS. THIS ARRAY IS OVERWRITTEN BY
! THE RESULTS.
!
! WORK
! A WORK ARRAY OF DIMENSION M*(N+1)
!
! TRIGS
! AN ARRAY SET UP BY SET99, WHICH MUST BE CALLED FIRST.
!
! IFAX
! AN ARRAY SET UP BY SET99, WHICH MUST BE CALLED FIRST.
!
! INC
! THE INCREMENT (IN WORDS) BETWEEN SUCCESSIVE ELEMENTS OF
! EACH DATA OR COEFFICIENT VECTOR (E.G. INC=1 FOR
! CONSECUTIVELY STORED DATA).
!
! JUMP
! THE INCREMENT (IN WORDS) BETWEEN THE FIRST ELEMENTS OF
! SUCCESSIVE DATA OR COEFFICIENT VECTORS. ON THE CRAY-1,
! TRY TO ARRANGE DATA SO THAT JUMP IS NOT A MULTIPLE OF 8
! (TO AVOID MEMORY BANK CONFLICTS). FOR CLARIFICATION OF
! INC AND JUMP, SEE THE EXAMPLES BELOW.
!
! N
! THE LENGTH OF EACH TRANSFORM (SEE DEFINITION OF
! TRANSFORMS, BELOW).
!
! M
! THE NUMBER OF TRANSFORMS TO BE DONE SIMULTANEOUSLY.
!
! ISIGN
! = +1 FOR A TRANSFORM FROM FOURIER COEFFICIENTS TO
! GRIDPOINT VALUES.
! = -1 FOR A TRANSFORM FROM GRIDPOINT VALUES TO FOURIER
! COEFFICIENTS.
!
! ON OUTPUT A
! IF ISIGN = +1, AND M COEFFICIENT VECTORS ARE SUPPLIED
! EACH CONTAINING THE SEQUENCE:
!
! A(0),B(0),A(1),B(1),...,A(N/2),B(N/2) (N+2 VALUES)
!
! THEN THE RESULT CONSISTS OF M DATA VECTORS EACH
! CONTAINING THE CORRESPONDING N+2 GRIDPOINT VALUES:
!
! FOR FFT991, X(0), X(1), X(2),...,X(N-1),0,0.
! FOR FFT99, X(N-1),X(0),X(1),X(2),...,X(N-1),X(0).
! (EXPLICIT CYCLIC CONTINUITY)
!
! WHEN ISIGN = +1, THE TRANSFORM IS DEFINED BY:
! X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
! AND I=SQRT (-1)
!
! IF ISIGN = -1, AND M DATA VECTORS ARE SUPPLIED EACH
! CONTAINING A SEQUENCE OF GRIDPOINT VALUES X(J) AS
! DEFINED ABOVE, THEN THE RESULT CONSISTS OF M VECTORS
! EACH CONTAINING THE CORRESPONDING FOURIER COFFICIENTS
! A(K), B(K), 0 .LE. K .LE N/2.
!
! WHEN ISIGN = -1, THE INVERSE TRANSFORM IS DEFINED BY:
! C(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*EXP(-2*I*J*K*PI/N))
! WHERE C(K)=A(K)+I*B(K) AND I=SQRT(-1)
!
! A CALL WITH ISIGN=+1 FOLLOWED BY A CALL WITH ISIGN=-1
! (OR VICE VERSA) RETURNS THE ORIGINAL DATA.
!
! NOTE: THE FACT THAT THE GRIDPOINT VALUES X(J) ARE REAL
! IMPLIES THAT B(0)=B(N/2)=0. FOR A CALL WITH ISIGN=+1,
! IT IS NOT ACTUALLY NECESSARY TO SUPPLY THESE ZEROS.
!
! EXAMPLES GIVEN 19 DATA VECTORS EACH OF LENGTH 64 (+2 FOR EXPLICIT
! CYCLIC CONTINUITY), COMPUTE THE CORRESPONDING VECTORS OF
! FOURIER COEFFICIENTS. THE DATA MAY, FOR EXAMPLE, BE
! ARRANGED LIKE THIS:
!
! FIRST DATA A(1)= . . . A(66)= A(70)
! VECTOR X(63) X(0) X(1) X(2) ... X(63) X(0) (4 EMPTY LOCATIONS)
!
! SECOND DATA A(71)= . . . A(140)
! VECTOR X(63) X(0) X(1) X(2) ... X(63) X(0) (4 EMPTY LOCATIONS)
!
! AND SO ON. HERE INC=1, JUMP=70, N=64, M=19, ISIGN=-1,
! AND FFT99 SHOULD BE USED (BECAUSE OF THE EXPLICIT CYCLIC
! CONTINUITY).
!
! ALTERNATIVELY THE DATA MAY BE ARRANGED LIKE THIS:
!
! FIRST SECOND LAST
! DATA DATA DATA
! VECTOR VECTOR VECTOR
!
! A(1)= A(2)= A(19)=
!
! X(63) X(63) . . . X(63)
! A(20)= X(0) X(0) . . . X(0)
! A(39)= X(1) X(1) . . . X(1)
! . . .
! . . .
! . . .
!
! IN WHICH CASE WE HAVE INC=19, JUMP=1, AND THE REMAINING
! PARAMETERS ARE THE SAME AS BEFORE. IN EITHER CASE, EACH
! COEFFICIENT VECTOR OVERWRITES THE CORRESPONDING INPUT
! DATA VECTOR.
!
!-----------------------------------------------------------------------
DIMENSION a((n+2)*lot),work((n+1)*lot),trigs(3*n/2+1),ifax(13)
!
! SUBROUTINE "FFT99" - MULTIPLE FAST REAL PERIODIC TRANSFORM
! CORRESPONDING TO OLD SCALAR ROUTINE FFT9
! PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM
! IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12
! (1970), 315-337)
!
! A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA
! WORK IS AN AREA OF SIZE (N+1)*LOT
! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2
! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'
! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
! N IS THE LENGTH OF THE DATA VECTORS
! LOT IS THE NUMBER OF DATA VECTORS
! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
!
! ORDERING OF COEFFICIENTS:
! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED
!
! ORDERING OF DATA:
! X(N-1),X(0),X(1),X(2),...,X(N),X(0)
! I.E. EXPLICIT CYCLIC CONTINUITY; (N+2) LOCATIONS REQUIRED
!
! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN
! PARALLEL
!
! *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER
!
! DEFINITION OF TRANSFORMS:
! -------------------------
!
! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
!
! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
!
!
!
!
nfax=ifax(1)
nx=n+1
nh=n/2
ink=inc+inc
IF (ISIGN == +1) GO TO 30
!
! IF NECESSARY, TRANSFER DATA TO WORK AREA
igo=50
IF (MOD(nfax,2) == 1) GO TO 40
ibase=inc+1
jbase=1
DO l=1,lot
i=ibase
j=jbase
!DIR$ IVDEP
DO m=1,n
work(j)=a(i)
i=i+inc
j=j+1
END DO
ibase=ibase+jump
jbase=jbase+nx
END DO
!
igo=60
GO TO 40
!
! PREPROCESSING (ISIGN=+1)
! ------------------------
!
30 CONTINUE
CALL fft99a
(a,work,trigs,inc,jump,n,lot)
igo=60
!
! COMPLEX TRANSFORM
! -----------------
!
40 CONTINUE
ia=inc+1
la=1
DO k=1,nfax
IF (igo == 60) GO TO 60
! 50 CONTINUE
CALL vpassm( ink,2,jump,nx,lot,nh,ifax(k+1),la)
igo=60
GO TO 70
60 CONTINUE
CALL vpassm(work(1),work(2),a(ia),a(ia+inc),trigs, &
2,ink,nx,jump,lot,nh,ifax(k+1),la)
igo=50
70 CONTINUE
la=la*ifax(k+1)
END DO
!
IF (ISIGN == -1) GO TO 130
!
! IF NECESSARY, TRANSFER DATA FROM WORK AREA
IF (MOD(nfax,2) == 1) GO TO 110
ibase=1
jbase=ia
DO l=1,lot
i=ibase
j=jbase
!DIR$ IVDEP
DO m=1,n
a(j)=work(i)
i=i+1
j=j+inc
END DO
ibase=ibase+nx
jbase=jbase+jump
END DO
!
! FILL IN CYCLIC BOUNDARY POINTS
110 CONTINUE
ia=1
ib=n*inc+1
!DIR$ IVDEP
DO l=1,lot
a(ia)=a(ib)
a(ib+inc)=a(ia+inc)
ia=ia+jump
ib=ib+jump
END DO
GO TO 140
!
! POSTPROCESSING (ISIGN=-1):
! --------------------------
!
130 CONTINUE
CALL fft99b(work,a,trigs,inc,jump,n,lot)
!
140 CONTINUE
RETURN
END SUBROUTINE fft99
SUBROUTINE fft99a(a,work,trigs,inc,jump,n,lot) 2
DIMENSION a(n),work(n),trigs(n)
!
! SUBROUTINE FFT99A - PREPROCESSING STEP FOR FFT99, ISIGN=+1
! (SPECTRAL TO GRIDPOINT TRANSFORM)
!
nh=n/2
nx=n+1
ink=inc+inc
!
! A(0) AND A(N/2)
ia=1
ib=n*inc+1
ja=1
jb=2
!DIR$ IVDEP
DO l=1,lot
work(ja)=a(ia)+a(ib)
work(jb)=a(ia)-a(ib)
ia=ia+jump
ib=ib+jump
ja=ja+nx
jb=jb+nx
END DO
!
! REMAINING WAVENUMBERS
iabase=2*inc+1
ibbase=(n-2)*inc+1
jabase=3
jbbase=n-1
!
DO k=3,nh,2
ia=iabase
ib=ibbase
ja=jabase
jb=jbbase
c=trigs(n+k)
s=trigs(n+k+1)
!DIR$ IVDEP
DO l=1,lot
work(ja)=(a(ia)+a(ib))- &
(s*(a(ia)-a(ib))+c*(a(ia+inc)+a(ib+inc)))
work(jb)=(a(ia)+a(ib))+ &
(s*(a(ia)-a(ib))+c*(a(ia+inc)+a(ib+inc)))
work(ja+1)=(c*(a(ia)-a(ib))-s*(a(ia+inc)+a(ib+inc)))+ &
(a(ia+inc)-a(ib+inc))
work(jb+1)=(c*(a(ia)-a(ib))-s*(a(ia+inc)+a(ib+inc)))- &
(a(ia+inc)-a(ib+inc))
ia=ia+jump
ib=ib+jump
ja=ja+nx
jb=jb+nx
END DO
iabase=iabase+ink
ibbase=ibbase-ink
jabase=jabase+2
jbbase=jbbase-2
END DO
!
IF (iabase /= ibbase) GO TO 50
! WAVENUMBER N/4 (IF IT EXISTS)
ia=iabase
ja=jabase
!DIR$ IVDEP
DO l=1,lot
work(ja)=2.0*a(ia)
work(ja+1)=-2.0*a(ia+inc)
ia=ia+jump
ja=ja+nx
END DO
!
50 CONTINUE
RETURN
END SUBROUTINE fft99a
SUBROUTINE fft99b(work,a,trigs,inc,jump,n,lot) 2
DIMENSION work(n),a(n),trigs(n)
!
! SUBROUTINE FFT99B - POSTPROCESSING STEP FOR FFT99, ISIGN=-1
! (GRIDPOINT TO SPECTRAL TRANSFORM)
!
nh=n/2
nx=n+1
ink=inc+inc
!
! A(0) AND A(N/2)
scale=1.0/FLOAT(n)
ia=1
ib=2
ja=1
jb=n*inc+1
!DIR$ IVDEP
DO l=1,lot
a(ja)=scale*(work(ia)+work(ib))
a(jb)=scale*(work(ia)-work(ib))
a(ja+inc)=0.0
a(jb+inc)=0.0
ia=ia+nx
ib=ib+nx
ja=ja+jump
jb=jb+jump
END DO
!
! REMAINING WAVENUMBERS
scale=0.5*scale
iabase=3
ibbase=n-1
jabase=2*inc+1
jbbase=(n-2)*inc+1
!
DO k=3,nh,2
ia=iabase
ib=ibbase
ja=jabase
jb=jbbase
c=trigs(n+k)
s=trigs(n+k+1)
!DIR$ IVDEP
DO l=1,lot
a(ja)=scale*((work(ia)+work(ib)) &
+(c*(work(ia+1)+work(ib+1))+s*(work(ia)-work(ib))))
a(jb)=scale*((work(ia)+work(ib)) &
-(c*(work(ia+1)+work(ib+1))+s*(work(ia)-work(ib))))
a(ja+inc)=scale*((c*(work(ia)-work(ib))-s*(work(ia+1)+work(ib+1))) &
+(work(ib+1)-work(ia+1)))
a(jb+inc)=scale*((c*(work(ia)-work(ib))-s*(work(ia+1)+work(ib+1))) &
-(work(ib+1)-work(ia+1)))
ia=ia+nx
ib=ib+nx
ja=ja+jump
jb=jb+jump
END DO
iabase=iabase+2
ibbase=ibbase-2
jabase=jabase+ink
jbbase=jbbase-ink
END DO
!
IF (iabase /= ibbase) GO TO 50
! WAVENUMBER N/4 (IF IT EXISTS)
ia=iabase
ja=jabase
scale=2.0*scale
!DIR$ IVDEP
DO l=1,lot
a(ja)=scale*work(ia)
a(ja+inc)=-scale*work(ia+1)
ia=ia+nx
ja=ja+jump
END DO
!
50 CONTINUE
RETURN
END SUBROUTINE fft99b
SUBROUTINE fft991(a,work,trigs,ifax,inc,jump,n,lot,ISIGN) 8,4
DIMENSION a(n),work(n),trigs(n),ifax(1)
!
! SUBROUTINE "FFT991" - MULTIPLE REAL/HALF-COMPLEX PERIODIC
! FAST FOURIER TRANSFORM
!
! SAME AS FFT99 EXCEPT THAT ORDERING OF DATA CORRESPONDS TO
! THAT IN MRFFT2
!
! PROCEDURE USED TO CONVERT TO HALF-LENGTH COMPLEX TRANSFORM
! IS GIVEN BY COOLEY, LEWIS AND WELCH (J. SOUND VIB., VOL. 12
! (1970), 315-337)
!
! A IS THE ARRAY CONTAINING INPUT AND OUTPUT DATA
! WORK IS AN AREA OF SIZE (N+1)*LOT
! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES
! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N/2
! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR'
! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA)
! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR
! N IS THE LENGTH OF THE DATA VECTORS
! LOT IS THE NUMBER OF DATA VECTORS
! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT
! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL
!
! ORDERING OF COEFFICIENTS:
! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED
!
! ORDERING OF DATA:
! X(0),X(1),X(2),...,X(N-1)
!
! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS IN
! PARALLEL
!
! *** N.B. N IS ASSUMED TO BE AN EVEN NUMBER
!
! DEFINITION OF TRANSFORMS:
! -------------------------
!
! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K)
!
! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
!
!
!
nfax=ifax(1)
nx=n+1
nh=n/2
ink=inc+inc
IF (ISIGN == +1) GO TO 30
!
! IF NECESSARY, TRANSFER DATA TO WORK AREA
igo=50
IF (MOD(nfax,2) == 1) GO TO 40
ibase=1
jbase=1
DO l=1,lot
i=ibase
j=jbase
!DIR$ IVDEP
DO m=1,n
work(j)=a(i)
i=i+inc
j=j+1
END DO
ibase=ibase+jump
jbase=jbase+nx
END DO
!
igo=60
GO TO 40
!
! PREPROCESSING (ISIGN=+1)
! ------------------------
!
30 CONTINUE
CALL fft99a(a,work,trigs,inc,jump,n,lot)
igo=60
!
! COMPLEX TRANSFORM
! -----------------
!
40 CONTINUE
ia=1
la=1
DO k=1,nfax
IF (igo == 60) GO TO 60
! 50 CONTINUE
CALL vpassm( ink,2,jump,nx,lot,nh,ifax(k+1),la)
igo=60
GO TO 70
60 CONTINUE
CALL vpassm(work(1),work(2),a(ia),a(ia+inc),trigs, &
2,ink,nx,jump,lot,nh,ifax(k+1),la)
igo=50
70 CONTINUE
la=la*ifax(k+1)
END DO
!
IF (ISIGN == -1) GO TO 130
!
! IF NECESSARY, TRANSFER DATA FROM WORK AREA
IF (MOD(nfax,2) == 1) GO TO 110
ibase=1
jbase=1
DO l=1,lot
i=ibase
j=jbase
!DIR$ IVDEP
DO m=1,n
a(j)=work(i)
i=i+1
j=j+inc
END DO
ibase=ibase+nx
jbase=jbase+jump
END DO
!
! FILL IN ZEROS AT END
110 CONTINUE
ib=n*inc+1
!DIR$ IVDEP
DO l=1,lot
a(ib)=0.0
a(ib+inc)=0.0
ib=ib+jump
END DO
GO TO 140
!
! POSTPROCESSING (ISIGN=-1):
! --------------------------
!
130 CONTINUE
CALL fft99b(work,a,trigs,inc,jump,n,lot)
!
140 CONTINUE
RETURN
END SUBROUTINE fft991
SUBROUTINE set99 (trigs, ifax, n) 8,3
DIMENSION ifax(13),trigs(3*n/2+1)
!
! MODE 3 IS USED FOR REAL/HALF-COMPLEX TRANSFORMS. IT IS POSSIBLE
! TO DO COMPLEX/COMPLEX TRANSFORMS WITH OTHER VALUES OF MODE, BUT
! DOCUMENTATION OF THE DETAILS WERE NOT AVAILABLE WHEN THIS ROUTINE
! WAS WRITTEN.
!
DATA mode /3/
CALL fax (ifax, n, mode)
i = ifax(1)
IF (ifax(i+1) > 5 .OR. n <= 4) ifax(1) = -99
IF (ifax(1) <= 0 ) THEN
WRITE(6,*) ' SET99 -- INVALID N'
CALL arpsstop('arpsstop called from SET99',1)
END IF
CALL fftrig (trigs, n, mode)
RETURN
END SUBROUTINE set99
SUBROUTINE fax(ifax,n,mode) 1
DIMENSION ifax(10)
nn=n
IF (IABS(mode) == 1) GO TO 10
IF (IABS(mode) == 8) GO TO 10
nn=n/2
IF ((nn+nn) == n) GO TO 10
ifax(1)=-99
RETURN
10 k=1
! TEST FOR FACTORS OF 4
20 IF (MOD(nn,4) /= 0) GO TO 30
k=k+1
ifax(k)=4
nn=nn/4
IF (nn == 1) GO TO 80
GO TO 20
! TEST FOR EXTRA FACTOR OF 2
30 IF (MOD(nn,2) /= 0) GO TO 40
k=k+1
ifax(k)=2
nn=nn/2
IF (nn == 1) GO TO 80
! TEST FOR FACTORS OF 3
40 IF (MOD(nn,3) /= 0) GO TO 50
k=k+1
ifax(k)=3
nn=nn/3
IF (nn == 1) GO TO 80
GO TO 40
! NOW FIND REMAINING FACTORS
50 l=5
inc=2
! INC ALTERNATELY TAKES ON VALUES 2 AND 4
60 IF (MOD(nn,l) /= 0) GO TO 70
k=k+1
ifax(k)=l
nn=nn/l
IF (nn == 1) GO TO 80
GO TO 60
70 l=l+inc
inc=6-inc
GO TO 60
80 ifax(1)=k-1
! IFAX(1) CONTAINS NUMBER OF FACTORS
nfax=ifax(1)
! SORT FACTORS INTO ASCENDING ORDER
IF (nfax == 1) GO TO 110
DO ii=2,nfax
istop=nfax+2-ii
DO i=2,istop
IF (ifax(i+1) >= ifax(i)) CYCLE
item=ifax(i)
ifax(i)=ifax(i+1)
ifax(i+1)=item
END DO
END DO
110 CONTINUE
RETURN
END SUBROUTINE fax
SUBROUTINE fftrig(trigs,n,mode) 1
DIMENSION trigs(1)
pi=2.0*ASIN(1.0)
imode=IABS(mode)
nn=n
IF (imode > 1.AND.imode < 6) nn=n/2
del=(pi+pi)/FLOAT(nn)
l=nn+nn
DO i=1,l,2
angle=0.5*FLOAT(i-1)*del
trigs(i)=COS(angle)
trigs(i+1)=SIN(angle)
END DO
IF (imode == 1) RETURN
IF (imode == 8) RETURN
del=0.5*del
nh=(nn+1)/2
l=nh+nh
la=nn+nn
DO i=1,l,2
angle=0.5*FLOAT(i-1)*del
trigs(la+i)=COS(angle)
trigs(la+i+1)=SIN(angle)
END DO
IF (imode <= 3) RETURN
del=0.5*del
la=la+nn
IF (mode == 5) GO TO 40
DO i=2,nn
angle=FLOAT(i-1)*del
trigs(la+i)=2.0*SIN(angle)
END DO
RETURN
40 CONTINUE
del=0.5*del
DO i=2,n
angle=FLOAT(i-1)*del
trigs(la+i)=SIN(angle)
END DO
RETURN
END SUBROUTINE fftrig
SUBROUTINE vpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la) 4
DIMENSION a(n),b(n),c(n),d(n),trigs(n)
!
! SUBROUTINE "VPASSM" - MULTIPLE VERSION OF "VPASSA"
! PERFORMS ONE PASS THROUGH DATA
! AS PART OF MULTIPLE COMPLEX FFT ROUTINE
! A IS FIRST REAL INPUT VECTOR
! B IS FIRST IMAGINARY INPUT VECTOR
! C IS FIRST REAL OUTPUT VECTOR
! D IS FIRST IMAGINARY OUTPUT VECTOR
! TRIGS IS PRECALCULATED TABLE OF SINES " COSINES
! INC1 IS ADDRESSING INCREMENT FOR A AND B
! INC2 IS ADDRESSING INCREMENT FOR C AND D
! INC3 IS ADDRESSING INCREMENT BETWEEN A"S & B"S
! INC4 IS ADDRESSING INCREMENT BETWEEN C"S & D"S
! LOT IS THE NUMBER OF VECTORS
! N IS LENGTH OF VECTORS
! IFAC IS CURRENT FACTOR OF N
! LA IS PRODUCT OF PREVIOUS FACTORS
!
DATA sin36/0.587785252292473/,cos36/0.809016994374947/, &
sin72/0.951056516295154/,cos72/0.309016994374947/, &
sin60/0.866025403784437/
!
m=n/ifac
iink=m*inc1
jink=la*inc2
jump=(ifac-1)*jink
ibase=0
jbase=0
igo=ifac-1
IF (igo > 4) RETURN
! GO TO (10,50,90,130),igo
!
! obsolescent feature correction temporarily by WYH.
!
SELECT CASE (igo)
CASE (1)
GO TO 10
CASE (2)
GO TO 50
CASE (3)
GO TO 90
CASE (4)
GO TO 130
CASE DEFAULT
! Do nothing
END SELECT
!
! end of correction by WYH.
!
!
! CODING FOR FACTOR 2
!
10 ia=1
ja=1
ib=ia+iink
jb=ja+jink
DO l=1,la
i=ibase
j=jbase
!DIR$ IVDEP
DO ijk=1,lot
c(ja+j)=a(ia+i)+a(ib+i)
d(ja+j)=b(ia+i)+b(ib+i)
c(jb+j)=a(ia+i)-a(ib+i)
d(jb+j)=b(ia+i)-b(ib+i)
i=i+inc3
j=j+inc4
END DO
ibase=ibase+inc1
jbase=jbase+inc2
END DO
IF (la == m) RETURN
la1=la+1
jbase=jbase+jump
DO k=la1,m,la
kb=k+k-2
c1=trigs(kb+1)
s1=trigs(kb+2)
DO l=1,la
i=ibase
j=jbase
!DIR$ IVDEP
DO ijk=1,lot
c(ja+j)=a(ia+i)+a(ib+i)
d(ja+j)=b(ia+i)+b(ib+i)
c(jb+j)=c1*(a(ia+i)-a(ib+i))-s1*(b(ia+i)-b(ib+i))
d(jb+j)=s1*(a(ia+i)-a(ib+i))+c1*(b(ia+i)-b(ib+i))
i=i+inc3
j=j+inc4
END DO
ibase=ibase+inc1
jbase=jbase+inc2
END DO
jbase=jbase+jump
END DO
RETURN
!
! CODING FOR FACTOR 3
!
50 ia=1
ja=1
ib=ia+iink
jb=ja+jink
ic=ib+iink
jc=jb+jink
DO l=1,la
i=ibase
j=jbase
!DIR$ IVDEP
DO ijk=1,lot
c(ja+j)=a(ia+i)+(a(ib+i)+a(ic+i))
d(ja+j)=b(ia+i)+(b(ib+i)+b(ic+i))
c(jb+j)=(a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)-b(ic+i)))
c(jc+j)=(a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)-b(ic+i)))
d(jb+j)=(b(ia+i)-0.5*(b(ib+i)+b(ic+i)))+(sin60*(a(ib+i)-a(ic+i)))
d(jc+j)=(b(ia+i)-0.5*(b(ib+i)+b(ic+i)))-(sin60*(a(ib+i)-a(ic+i)))
i=i+inc3
j=j+inc4
END DO
ibase=ibase+inc1
jbase=jbase+inc2
END DO
IF (la == m) RETURN
la1=la+1
jbase=jbase+jump
DO k=la1,m,la
kb=k+k-2
kc=kb+kb
c1=trigs(kb+1)
s1=trigs(kb+2)
c2=trigs(kc+1)
s2=trigs(kc+2)
DO l=1,la
i=ibase
j=jbase
!DIR$ IVDEP
DO ijk=1,lot
c(ja+j)=a(ia+i)+(a(ib+i)+a(ic+i))
d(ja+j)=b(ia+i)+(b(ib+i)+b(ic+i))
c(jb+j)= &
c1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)-b(ic+i)))) &
-s1*((b(ia+i)-0.5*(b(ib+i)+b(ic+i)))+(sin60*(a(ib+i)-a(ic+i))))
d(jb+j)= &
s1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)-b(ic+i)))) &
+c1*((b(ia+i)-0.5*(b(ib+i)+b(ic+i)))+(sin60*(a(ib+i)-a(ic+i))))
c(jc+j)= &
c2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)-b(ic+i)))) &
-s2*((b(ia+i)-0.5*(b(ib+i)+b(ic+i)))-(sin60*(a(ib+i)-a(ic+i))))
d(jc+j)= &
s2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)-b(ic+i)))) &
+c2*((b(ia+i)-0.5*(b(ib+i)+b(ic+i)))-(sin60*(a(ib+i)-a(ic+i))))
i=i+inc3
j=j+inc4
END DO
ibase=ibase+inc1
jbase=jbase+inc2
END DO
jbase=jbase+jump
END DO
RETURN
!
! CODING FOR FACTOR 4
!
90 ia=1
ja=1
ib=ia+iink
jb=ja+jink
ic=ib+iink
jc=jb+jink
id=ic+iink
jd=jc+jink
DO l=1,la
i=ibase
j=jbase
!DIR$ IVDEP
DO ijk=1,lot
c(ja+j)=(a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i))
c(jc+j)=(a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))
d(ja+j)=(b(ia+i)+b(ic+i))+(b(ib+i)+b(id+i))
d(jc+j)=(b(ia+i)+b(ic+i))-(b(ib+i)+b(id+i))
c(jb+j)=(a(ia+i)-a(ic+i))-(b(ib+i)-b(id+i))
c(jd+j)=(a(ia+i)-a(ic+i))+(b(ib+i)-b(id+i))
d(jb+j)=(b(ia+i)-b(ic+i))+(a(ib+i)-a(id+i))
d(jd+j)=(b(ia+i)-b(ic+i))-(a(ib+i)-a(id+i))
i=i+inc3
j=j+inc4
END DO
ibase=ibase+inc1
jbase=jbase+inc2
END DO
IF (la == m) RETURN
la1=la+1
jbase=jbase+jump
DO k=la1,m,la
kb=k+k-2
kc=kb+kb
kd=kc+kb
c1=trigs(kb+1)
s1=trigs(kb+2)
c2=trigs(kc+1)
s2=trigs(kc+2)
c3=trigs(kd+1)
s3=trigs(kd+2)
DO l=1,la
i=ibase
j=jbase
!DIR$ IVDEP
DO ijk=1,lot
c(ja+j)=(a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i))
d(ja+j)=(b(ia+i)+b(ic+i))+(b(ib+i)+b(id+i))
c(jc+j)= &
c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) &
-s2*((b(ia+i)+b(ic+i))-(b(ib+i)+b(id+i)))
d(jc+j)= &
s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) &
+c2*((b(ia+i)+b(ic+i))-(b(ib+i)+b(id+i)))
c(jb+j)= &
c1*((a(ia+i)-a(ic+i))-(b(ib+i)-b(id+i))) &
-s1*((b(ia+i)-b(ic+i))+(a(ib+i)-a(id+i)))
d(jb+j)= &
s1*((a(ia+i)-a(ic+i))-(b(ib+i)-b(id+i))) &
+c1*((b(ia+i)-b(ic+i))+(a(ib+i)-a(id+i)))
c(jd+j)= &
c3*((a(ia+i)-a(ic+i))+(b(ib+i)-b(id+i))) &
-s3*((b(ia+i)-b(ic+i))-(a(ib+i)-a(id+i)))
d(jd+j)= &
s3*((a(ia+i)-a(ic+i))+(b(ib+i)-b(id+i))) &
+c3*((b(ia+i)-b(ic+i))-(a(ib+i)-a(id+i)))
i=i+inc3
j=j+inc4
END DO
ibase=ibase+inc1
jbase=jbase+inc2
END DO
jbase=jbase+jump
END DO
RETURN
!
! CODING FOR FACTOR 5
!
130 ia=1
ja=1
ib=ia+iink
jb=ja+jink
ic=ib+iink
jc=jb+jink
id=ic+iink
jd=jc+jink
ie=id+iink
je=jd+jink
DO l=1,la
i=ibase
j=jbase
!DIR$ IVDEP
DO ijk=1,lot
c(ja+j)=a(ia+i)+(a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))
d(ja+j)=b(ia+i)+(b(ib+i)+b(ie+i))+(b(ic+i)+b(id+i))
c(jb+j)=(a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i))) &
-(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))
c(je+j)=(a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i))) &
+(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))
d(jb+j)=(b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i))) &
+(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i)))
d(je+j)=(b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i))) &
-(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i)))
c(jc+j)=(a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i))) &
-(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))
c(jd+j)=(a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i))) &
+(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))
d(jc+j)=(b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i))) &
+(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i)))
d(jd+j)=(b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i))) &
-(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i)))
i=i+inc3
j=j+inc4
END DO
ibase=ibase+inc1
jbase=jbase+inc2
END DO
IF (la == m) RETURN
la1=la+1
jbase=jbase+jump
DO k=la1,m,la
kb=k+k-2
kc=kb+kb
kd=kc+kb
ke=kd+kb
c1=trigs(kb+1)
s1=trigs(kb+2)
c2=trigs(kc+1)
s2=trigs(kc+2)
c3=trigs(kd+1)
s3=trigs(kd+2)
c4=trigs(ke+1)
s4=trigs(ke+2)
DO l=1,la
i=ibase
j=jbase
!DIR$ IVDEP
DO ijk=1,lot
c(ja+j)=a(ia+i)+(a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))
d(ja+j)=b(ia+i)+(b(ib+i)+b(ie+i))+(b(ic+i)+b(id+i))
c(jb+j)= &
c1*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i))) &
-(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))) &
-s1*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i))) &
+(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
d(jb+j)= &
s1*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i))) &
-(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))) &
+c1*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i))) &
+(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
c(je+j)= &
c4*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i))) &
+(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))) &
-s4*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i))) &
-(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
d(je+j)= &
s4*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i))) &
+(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))) &
+c4*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i))) &
-(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
c(jc+j)= &
c2*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i))) &
-(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))) &
-s2*((b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i))) &
+(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
d(jc+j)= &
s2*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i))) &
-(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))) &
+c2*((b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i))) &
+(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
c(jd+j)= &
c3*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i))) &
+(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))) &
-s3*((b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i))) &
-(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
d(jd+j)= &
s3*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i))) &
+(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))) &
+c3*((b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i))) &
-(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
i=i+inc3
j=j+inc4
END DO
ibase=ibase+inc1
jbase=jbase+inc2
END DO
jbase=jbase+jump
END DO
RETURN
END SUBROUTINE vpassm