1 C    -------------------------------------------------------------------------
  2 C    |  File name:   fftd.for              Sean O'Connor             9-22-87 |
  3 C    |                                                               7-01-17 |
  4 C    ------------------------------------------------------------------------- 
  5 
  6 C    -------------------------------------------------------------------------
  7 C    |  VAX 11/785       VMS 4.5                   VAX EXTENDED FORTRAN 77   |
  8 C    |  MacBook Pro      Ubuntu LTS 16.04 Linux    gfortran                  |
  9 C    -------------------------------------------------------------------------
 10 
 11 
 12 C  ****************************************************************************
 13 C  *                                                                          *
 14 C  *                            FFTD                                          *
 15 C  *                                                                          *
 16 C  ****************************************************************************
 17 C
 18 C
 19 C  DESCRIPTION
 20 C 
 21 C      This program tests the subroutine FFT in the file FFT.FOR  It does a 
 22 C      forward discrete Fourier transform (DFT).  It then does an inverse
 23 C      transform of the forward transform to check the roundoff error.
 24 C
 25 C
 26 C  INPUT 
 27 C
 28 C      The number of complex points in the transform.
 29 C      Whether to do a forward (1) or inverse (-1) transform.
 30 C      The data points (complex numbers).
 31 C
 32 C
 33 C  OUTPUT
 34 C
 35 C      The DFT of the input.
 36 C      The inverse DFT of the first DFT.
 37 C
 38 C
 39 C  EXAMPLE 
 40 C
 41 C      $ ./FFTD
 42 C
 43 C      Enter number of complex points in transform: 4
 44 C
 45 C      Enter 1 (forward) or -1 (inverse) transform:  1
 46 C
 47 C      Enter point no.   0 (Real, Imag) >  (1,1)
 48 C      Enter point no.   1 (Real, Imag) >  (2,2)
 49 C      Enter point no.   2 (Real, Imag) >  (3,3)
 50 C      Enter point no.   3 (Real, Imag) >  (4,4)
 51 C
 52 C      num_points =          4
 53 C      direction  =          1
 54 C
 55 C      Point No.   0 = (   1.0000000000000       1.0000000000000     ) 
 56 C      Point No.   1 = (   2.0000000000000       2.0000000000000     ) 
 57 C      Point No.   2 = (   3.0000000000000       3.0000000000000     ) 
 58 C      Point No.   3 = (   4.0000000000000       4.0000000000000     ) 
 59 C
 60 C      The transform of the input is . . .
 61 C 
 62 C      Point No.   0= (  5.0000000000000       5.0000000000000     ) 
 63 C      Point No.   1= ( -2.0000000000000      0.00000000000000E+000) 
 64 C      Point No.   2= ( -1.0000000000000      -1.0000000000000     ) 
 65 C      Point No.   3= (-0.59604644775391E-007 -2.0000000000000     ) 
 66 C
 67 C      The inverse transform of the transform above is . . .
 68 C  
 69 C      Point No.   0= (  1.0000000000000       1.0000000000000     ) 
 70 C      Point No.   1= (  2.0000000000000       2.0000000000000     ) 
 71 C      Point No.   2= (  3.0000000000000       3.0000000000000     ) 
 72 C      Point No.   3= (  4.0000000000000       4.0000000000000     ) 
 73 C
 74 C
 75 C  INSTALLATION
 76 C
 77 C      To compile, type
 78 C
 79 C          $ gfortran FFTD.FOR FORTRAN FFT.FOR -o FFT
 80 C
 81 C      To run, type
 82 C
 83 C          $ ./FFTD
 84 C       
 85 C 
 86 C  ****************************************************************************
 87 C  *                                                                          *
 88 C  ****************************************************************************
 89 
 90       PROGRAM fftd
 91 
 92 C  -------------------------------------------------------------------------
 93 C  |                            MACRO DEFINITIONS                          |
 94 C  ------------------------------------------------------------------------- 
 95 
 96       PARAMETER ( MAX_NUM_POINTS = 255 )   !  Maximum number of points in FFT.
 97 
 98 
 99 C -------------------------------------------------------------------------
100 C |                          LOCAL VARIABLES                              |
101 C -------------------------------------------------------------------------
102 
103         INTEGER i,             !  Loop counter
104      $          num_points,    !  Number of points in the transform
105      $          direction      !  Direction of transform 1 for forward,
106                                ! -1 for inverse
107 
108         COMPLEX X( 0 : MAX_NUM_POINTS )  !  Array of complex data points
109                                          !  to be transformed.
110 
111 
112 C -------------------------------------------------------------------------
113 C |                          PROGRAM CODE                                 |
114 C -------------------------------------------------------------------------
115 
116 C
117 C    Read the number of points in the FFT
118 C
119         WRITE( 6, 100 )
120 100     FORMAT ( ' Enter number of complex points in transform:  ' )
121 
122         READ ( 5 , * ) num_points      ! Free format input
123 
124 
125 
126 C
127 C       Decide whether to do the forward or inverse FFT
128 C
129 
130         WRITE( 6, 200 )
131 200     FORMAT ( //, ' Enter 1 for a forward transform ', /,
132      $               ' or -1 for an inverse transform: ' )
133 
134         READ ( 5, * ) direction     ! Free format input
135 
136 
137 C
138 C       Check to see that the number of points in the FFT
139 C       is valid.  If not, halt.
140 C
141 
142         IF ( (num_points .LE. 0) .OR.
143      $       (num_points .GT. MAX_NUM_POINTS) ) THEN
144 
145                 WRITE( 6, 300 )
146 300             FORMAT( ' The number of points is out of range ')
147                 STOP
148 
149         END IF
150 
151 
152 
153 C
154 C       Read in the complex data points, one by one 
155 C
156         DO i = 0, num_points - 1
157 
158             WRITE( 6, 400 ) i
159 400         FORMAT ( ' Enter point no. ', I3,' (Real, Imag) > ')
160             READ ( 5 , * ) X( i )      ! Free format input
161 
162         END DO
163 
164 
165 C
166 C       Print out the points just read in.
167 C
168 
169         WRITE( 6, 500 ) num_points
170 500     FORMAT ( //, ' num_points = ', I10, / )
171 
172         WRITE( 6, 600 ) direction
173 600     FORMAT ( ' direction =  ', I10, // )
174 
175         DO i = 0, num_points - 1
176 
177             WRITE( 6 , 700 ) i , X ( i )
178 700         FORMAT ( ' Point no. ', I3, ' = ( ', 2G22.14E3, ') ' )
179 
180         END DO
181 
182 
183 C
184 C       Call the FFT subroutine.
185 C
186         CALL FFT( num_points, X, direction )
187 
188 
189 C
190 C       Print the transform results.
191 C
192         WRITE( 6, 800 )
193 800     FORMAT ( //, ' The transform of the input data is . . . ', // )
194 
195 
196         DO i = 0, num_points - 1
197 
198             WRITE( 6, 900 ) i , X( i )
199 900         FORMAT (' Point No. ', I3, '= (', 2G22.14E3, ') ' )
200 
201         END DO
202 
203 
204 C
205 C       Call the FFT subroutine with -direction to do the inverse transform.
206 C
207         CALL FFT( num_points, X, -direction )
208 
209 
210 
211 C
212 C       Print the inverse transform coefficients.
213 C
214         WRITE ( 6, 1000 )
215 1000    FORMAT ( //,' The inverse transform ',
216      $               'of the transform above is . . .', // )
217 
218 
219         DO i = 0, num_points - 1
220 
221             WRITE ( 6, 1100 ) i , X( i )
222 1100        FORMAT (' Point No. ', I3, '= (', 2G22.14E3, ') ' )
223 
224         END DO
225 
226 
227 
228 C
229 C  Test the error checking in the FFT routine.
230 C
231 
232 C  Number of points is too large.  Currently, MAX_TABLE_SIZE in
233 C  subroutine FFT is 20, so we can transform up to 2 ** 21 points.
234 
235         PRINT *, ' Calling FFT with too many points.'
236 
237         CALL FFT( 2 ** 22, X, 1 )
238 
239 
240 C  Number of points is less than 2.
241 
242         PRINT *, ' Calling FFT with 1 point.'
243 
244         CALL FFT( 1, X, 1 )
245 
246 
247 C  Calling FFT with 2^ 15 + 1 points (not a power of 2)
248 
249         PRINT *, ' Calling FFT with 2 ^ 15 + 1 points'
250 
251         CALL FFT( 2 * 15 + 1, X, 1 )
252 
253 
254 C  Direction is not 1 or -1.
255 
256         PRINT *, ' Calling FFT with direction = 2'
257 
258         CALL FFT( 16, X, 2)
259 
260 
261         STOP
262         END   ! of program fftd