SUBROUTINE CFFT2D(N,F,LDF,W,FORWD)
C***BEGIN PROLOGUE CFFT2D
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) 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 CFFT2D. For example,
C if you declare F by either
C COMPLEX F(0:15,0:20) OR F(16,21) then
C set LDF=16. You must have LDF >= N, NOT CHECKED.
C W: (COMPLEX) Array for internal use as work storage.
C Must be dimensioned by calling program to be
C at least 6N+15 REAL words or 3N+8 COMPLEX 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) Forward or reverse transformed input matrix.
C Output is unscaled, that is, a call to CFFT2D
C with FORWD=.TRUE. followed by a call to CFFT2D
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.
C
C
C***REFERENCES (NONE)
C***ROUTINES CALLED CFFTI, CFFTF, CFFTB
C***END PROLOGUE CFFT2D
COMPLEX F(0:LDF-1,0:*), W(0:*)
LOGICAL FORWD
C***FIRST EXECUTABLE STATEMENT CFFT2D
C Find transform of each row, then of each column
C First row transforms
CALL CFFTI(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 CFFTF(N,W,W(N))
ELSE
CALL CFFTB(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 CFFTF or CFFTB by passing F(0,J).
C Similarly CFFTF or CFFTB places results back there.
IF(FORWD) THEN
CALL CFFTF(N, F(0,J), W(N))
ELSE
CALL CFFTB(N, F(0,J), W(N))
ENDIF
20 CONTINUE
RETURN
END