module fft_work ! ! Wrapper/interface for fft routines ! ! (c) Copyright 1991 to 1998 by Michael A. Beer, William D. Dorland, ! P. B. Snyder, Q. P. Liu, and Gregory W. Hammett. ALL RIGHTS RESERVED. ! implicit none ! ccfft ! ! Generic: ! table(ccfft_table0 + ccfft_table1*n) ! work( ccfft_work0 + ccfft_work1 *n) ! work(2*n) ! UNICOS: ! table(100 + 8*n) ! work(8*n) ! UNICOS/mk: ! table(2*n) ! work(4*n) ! Origin: ! table(30 + 2*n) ! work(2*n) integer, parameter :: ccfft_table0 = 1 integer, parameter :: ccfft_table1 = 0 integer, parameter :: ccfft_work0 = 0 integer, parameter :: ccfft_work1 = 1 ! csfft/scfft ! UNICOS: ! table(100 + 4*n) ! work(4 + 4*n) ! UNICOS/mk: ! table(2*n) ! work(2*n) ! Origin: ! table(15 + n) ! work(n) integer, parameter :: csfft_table0 = 1 integer, parameter :: csfft_table1 = 0 integer, parameter :: csfft_work0 = 0 integer, parameter :: csfft_work1 = 1 contains ! translate calls to Cray FFT subroutines to NAG subroutines. ! On a 450 MHz DEC Alpha, these are about 1.5 times slower than the ! T3E or Origin2000. It may be possible to speed this up by switching ! to NAG routines which precompute and save some of the trig work ! arrays, or by switching to the native fft routines in the DEC dxml ! library (type "man dxml"). subroutine ccfft(isign,n,scale,x,y,table,work,isys) !!$DESCRIPTION (From Cray Scientific Library for ccfft): !!$ CCFFT computes the Fast Fourier Transform (FFT) of the complex vector !!$ x, and it stores the result in vector y. !!$ !!$ In FFT applications, it is customary to use zero-based subscripts; the !!$ formulas are simpler that way. Suppose that the arrays are !!$ dimensioned as follows: !!$ !!$ COMPLEX X(0:N-1), Y(0:N-1) !!$ !!$ The output array is the FFT of the input array, using the following !!$ formula for the FFT: !!$ !!$ n-1 !!$ Y(k) = scale * Sum [ X(j)*w**(isign*j*k) ] for k = 0, ..., n-1 !!$ j=0 !!$ !!$ where: !!$ w = exp(2*pi*i/n), !!$ i = + sqrt(-1), !!$ pi = 3.14159..., !!$ isign = +1 or -1 integer :: isign, n, isys real :: scale complex, dimension(0:n-1) :: x, y real, dimension(*) :: table real, dimension(*) :: work real, dimension(n) :: rx, ry integer ifail if(isign == 0) return ! NAG doesn't need initialization ifail = 0 rx(:) = real(x(:)) ry(:) = -isign*aimag(x(:)) call c06fce(rx,ry,n,work,ifail) y(:) = scale*sqrt(n*1.0d00)*cmplx(rx(:),-isign*ry(:)) return end subroutine ccfft subroutine csfft(isign,n,scale,z,x,table,work,isys) integer :: isign, n, isys real :: scale real, dimension(0:n-1) :: x complex, dimension(0:n/2) :: z real, dimension(*) :: table real, dimension(*) :: work integer ifail, k if(isign == 0) return ! NAG doesn't need initialization ifail = 0 do k=0,n/2 x(k) = real(z(k))*scale*sqrt(n*1.0d00) enddo do k=1,(n-1)/2 x(n-k) = -isign*aimag(z(k))*scale*sqrt(n*1.0d00) enddo call c06fbe(x, n, work, ifail) return end subroutine csfft subroutine scfft(isign,n,scale,x,z,table,work,isys) integer :: isign, n, isys real :: scale real, dimension(0:n-1) :: x complex, dimension(0:n/2) :: z real, dimension(*) :: table real, dimension(*) :: work real rx(0:n-1) integer ifail, k if(isign == 0) return ! NAG doesn't need initialization ifail = 0 rx(:) = x(:)*scale*sqrt(n*1.0) call c06fae(rx, n, work, ifail) z(0)=cmplx(rx(0),0.0) do k=1,(n-1)/2 z(k)=cmplx(rx(k),-isign*rx(n-k)) enddo if( n/2 > (n-1)/2 ) z(n/2) = cmplx(rx(n-k),0.0) return end subroutine scfft end module fft_work