Tweet
SUBROUTINE SPCFFT(U,N,ISIGN,WORK,INTERP)
* VARIABLES
* ---------
IMPLICIT NONE
LOGICAL*1
| INU ! Flag for calling SUBROUTINE SPCPFT( arguments ).
| ,SCALE ! .TRUE.=inverse transform -- .FALSE.=forward transform
INTEGER*4
| A ! After \
| ,B ! Before > Factors of N.
| ,C ! Current /
| ,N ! Length of the array to be transformed.
| ,I ! DO loop index.
| ,ISIGN ! sign of transform
REAL*4
| INTERP ! interpolation factor
COMPLEX*8
| U(*) ! Vector to be transformed
| ,WORK(*) ! Working storage.
* Initialize parameters.
A = 1
B = N
C = 1
INU = .TRUE.
IF (ISIGN.EQ.1) THEN
SCALE = .TRUE.
ELSE IF (ISIGN.EQ.-1) THEN
SCALE = .FALSE.
END IF
* Calculate Fourier transform by means of Glassman's algorithm
DO WHILE ( B .GT. 1 )
A = C * A
* Start of routine to get factors of N
C = 2
* Increment C until it is an integer factor of B
DO WHILE ( MOD(B,C) .NE. 0 )
C = C + 1
END DO
* Calculate largest factor of B
B = B / C
* Call Glassman's Fast Fourier Transform routine
IF ( INU ) THEN
CALL SPCPFT (A,B,C,U,WORK,ISIGN)
ELSE
CALL SPCPFT (A,B,C,WORK,U,ISIGN)
END IF
* Set flag to swap input & output arrays to SPCPFT
INU = ( .NOT. INU )
END DO
* If odd number of calls to SPCPFT swap WORK into U
IF ( .NOT. INU ) THEN
DO I = 1, N
U(I) = WORK(I)
END DO
END IF
* Scale the output for inverse Fourier transforms.
IF ( SCALE ) THEN
DO I = 1, N
U(I) = U(I) / (N/INTERP)
END DO
END IF
* TERMINATION
* -----------
RETURN
END
*=======================================================================
* Single Precision Complex Prime Factor Transform
*
* REFERENCES
* ----------
*
* Calling routines:
* SPCFFT
*
* Subroutines called:
* -none-
*
* Functions called:
* CMLPX
* COS
* SIN
* FLOAT
*=======================================================================
SUBROUTINE SPCPFT( A, B, C, UIN, UOUT, ISIGN )
* VARIABLES
* ---------
IMPLICIT NONE
INTEGER*4
| ISIGN ! ISIGN of the Fourier transform.
| ,A ! After \
| ,B ! Before > Factors of N.
| ,C ! Current /
| ,IA ! \
| ,IB ! \ DO loop indicies.
| ,IC ! /
| ,JCR ! /
| ,JC ! Dummy index.
REAL*8
| ANGLE
COMPLEX*8
| UIN(B,C,A) ! Input vector.
| ,UOUT(B,A,C) ! Output vector.
| ,DELTA ! Fourier transform kernel.
| ,OMEGA ! Multiples of exp( i TWOPI ).
| ,SUM ! Dummy register for addition for UOUT(B,A,C)
* ALGORITHM
* ---------
* Initialize run time parameters.
ANGLE =6.28318530717958 / FLOAT( A * C )
OMEGA = CMPLX( 1.0, 0.0 )
* Check the ISIGN of the transform.
IF( ISIGN .EQ. 1 ) THEN
DELTA = CMPLX( DCOS(ANGLE), DSIN(ANGLE) )
ELSE
DELTA = CMPLX( DCOS(ANGLE), -DSIN(ANGLE) )
END IF
* Do the computations.
DO IC = 1, C
DO IA = 1, A
DO IB = 1, B
SUM = UIN( IB, C, IA )
DO JCR = 2, C
JC = C + 1 - JCR
SUM = UIN( IB, JC, IA ) + OMEGA * SUM
END DO
UOUT( IB, IA, IC ) = SUM
END DO
OMEGA = DELTA * OMEGA
END DO
END DO
* TERMINATION
* -----------
RETURN
END
#ifdef FFT_TEST
program test
implicit none
integer*4 nmax
parameter( nmax=10000 )
complex*8 input(nmax)
complex*8 scratch(nmax)
real*4 a , b
integer*4 i
integer*4 n
n = 0
1000 format( 2e25.15 )
do i=1, nmax
read( 5 ,1000 ,end=10 ) input(i)
n = n + 1
end do
10 continue
call spcfft( input ,n ,0 ,scratch ,1.0 )
do i=1 ,n
write( 6 ,1000 ) real(input(i)) ,imag(input(i))
end do
call spcfft( input ,n ,1 ,scratch ,1.0 )
do i=1 ,n
write( 6 ,1000 ) real(input(i)) ,imag(input(i))
end do
end
#endif
#ifdef TEST
program test
implicit none
real*4 data_in(10000)
real*4 data_out(10000)
integer*4 i
integer*4 m
integer*4 n
read( 5 ,* ) m ,n
do i=1,m
read( 5 ,* ,end=10 ) data_in(i)
end do
10 continue
write( 6 ,* ) '#input'
do i=1 ,m
write( 6 ,'(2f25.10)' ) real(i-1)/(m) ,data_in(i)
end do
call resample( data_in ,m ,data_out ,n )
write( 6 ,* ) '#output'
do i=1 ,n
write( 6 ,'(2f25.10)' ) real(i-1)/(n) ,data_out(i)
end do
end
#endif
). I still maintain it is impossible to use something from the GoF book in C.
). I still maintain it is impossible to use something from the GoF book in C.begin
integer i;
real procedure sum (i, lo, hi, term);
value lo, hi;
integer i, lo, hi;
real term;
comment term is passed by-name, and so is i;
begin
real temp;
temp := 0;
for i := lo step 1 until hi do
temp := temp + term;
sum := temp
end;
print (sum (i, 1, 100, 1/i))
end
Comment