Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_rfft
8 : use m_juDFT
9 : PRIVATE
10 : ! module also contains routines vrffti,vrfftf&vrfftb as private routines below
11 : PUBLIC rfft
12 : CONTAINS
13 0 : SUBROUTINE rfft( &
14 : isn,n1d,n2d,n3d,n1,n2,n3, &
15 0 : nw1,nw2,nw3,wsave,b, &
16 0 : a)
17 : ! **********************************************************************
18 :
19 : ! isn = +1 : FFT a real, 3dim array "a" of dimensions n1*n2*n3 to a complex
20 : ! array of size n1*n2*(n3+1)/2 [if n3 is odd] or n1*n2*(n3/2+1)
21 : ! [n3 is even].
22 : ! isn = -1 : back-FFT of a complex array "a" of size n1*n2*(n3+1)/2 [odd]
23 : ! or n1*n2*(n3/2+1) [even] to a real, 3dim array
24 :
25 : ! the actual array is assumed to be located between
26 : ! 1 ... nw1+1 & n1-nw1+1 ... n1
27 : ! 1 ... nw2+1 & n2-nw2+1 ... n2
28 : ! 1 ... nw3+1 & n3-nw3+1 ... n3
29 : ! and padded with zeros in between.
30 : ! if nw1 >= (n1-1)/2, no padding is assumed (-> aliasing errors!)
31 :
32 : ! G.Bihlmayer (UniWien)
33 :
34 : ! **********************************************************************
35 : USE m_constants
36 : USE m_cfft
37 : IMPLICIT NONE
38 :
39 : INTEGER :: n1d,n2d,n3d,n1,n2,n3,nw1,nw2,nw3,isn
40 : REAL :: a(n1d,n2d,0:n3d),b(n1d,n2d,n3d),wsave(n3d+15)
41 :
42 : INTEGER :: i1,i2,i3,nup
43 : REAL :: factor
44 : LOGICAL :: l_nopad
45 :
46 : ! a ... array for FFT
47 : ! b ... work array
48 : ! wsave ... stores tables for r-2-c FFT
49 : ! n1,n2,n3 ... dimensions to be transformed
50 : ! nw1,nw2,nw3 ... actual dimensions of a before FFT
51 : ! n1d,n2d,n3d ... dimensions of a,b
52 :
53 : ! check for input errors
54 :
55 0 : IF ((isn/=-1) .AND. (isn /= 1)) CALL juDFT_error("choose isn=+/- 1" &
56 0 : ,calledby ="rfft")
57 0 : IF ((n1d < n1) .OR. (n2d < n2) .OR. (n3d < n3)) THEN
58 0 : WRITE (oUnit,*) 'n1d,n2d,n3d =',n1d,n2d,n3d
59 0 : WRITE (oUnit,*) 'n1 ,n2 ,n3 =',n1 ,n2 ,n3
60 0 : CALL juDFT_error("n(i) > n(i)d",calledby ="rfft")
61 : ENDIF
62 : IF ((n1 <= 2*nw1+1) .OR. &
63 0 : (n2 <= 2*nw2+1) .OR. &
64 : (n3 <= 2*nw3+1)) THEN
65 : ! WRITE (oUnit,*) 'n1 ,n2 ,n3 =',n1 ,n2 ,n3
66 : ! WRITE (oUnit,*) 'nw1,nw2,nw3 =',nw1,nw2,nw3
67 : l_nopad= .TRUE.
68 : ELSE
69 0 : l_nopad= .FALSE.
70 : ENDIF
71 :
72 : ! ******************** f o r w a r d - t r a n s f o r m *******************
73 :
74 0 : IF (isn == 1) THEN
75 :
76 : ! array a is assumed to be zero from (1,1,0) to (n1d,n2d,0) and the array
77 : ! to be FFT'd starts at (1,1,1) as n1*n2*n3 real numbers.
78 : ! first transform n1*n2 real sequences of lenghth n3 to n3/2 complex values
79 :
80 0 : CALL vrffti(n3,wsave)
81 0 : IF (l_nopad) THEN
82 0 : CALL vrfftf(n1*n2,n3,a(1,1,1),b,n1d*n2d,wsave)
83 : ELSE
84 0 : DO i2=1,nw2+1
85 0 : CALL vrfftf(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
86 0 : CALL vrfftf(nw1 ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
87 : ENDDO
88 0 : DO i2=n2-nw2+1,n2
89 0 : CALL vrfftf(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
90 0 : CALL vrfftf(nw1 ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
91 : ENDDO
92 : ENDIF
93 :
94 : ! now we have the FFT'd array stored as described in vrfftf
95 : ! (mixed real & compex data)
96 : ! remove the norm 1/sqrt(n3) (to be compatible with cfft)
97 :
98 0 : factor = sqrt(1.0*n3)
99 : ! ALL CPP_BLAS_sscal(n1d*n2d*n3,factor,a(1,1,1),1)
100 0 : a=a*factor
101 :
102 : ! now, the real part of f(0) has to be moved to a(n1,n2,0) to get a purely
103 : ! complex array starting at a(1,1,0)
104 :
105 0 : DO i1=1,n1
106 0 : DO i2=1,n2
107 0 : a(i1,i2,0)=a(i1,i2,1)
108 0 : a(i1,i2,1)=0.0
109 : ENDDO
110 : ENDDO
111 :
112 : ENDIF
113 :
114 : ! ******************** g e n e r a l p a r t *******************
115 :
116 : ! now perform n2*n3/2 and n1*n3/2 complex FFT's; a is assumed to be
117 : ! complex and starting at (1,1,0)
118 :
119 0 : IF (ABS((n3/2.)-NINT(n3/2.)) > 0.1) THEN
120 : nup = n3
121 : ELSE
122 0 : nup = n3+1
123 0 : IF (n3+1>n3d) CALL juDFT_error("n3 even & n3+1 > n3d" ,calledby &
124 0 : ="rfft")
125 0 : a(:,:,n3+1)=0.0
126 : ENDIF
127 :
128 0 : IF (l_nopad) THEN
129 0 : DO i3=1,nup,2
130 0 : CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n1,n1,isn)
131 0 : CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n2,n1*n2,isn)
132 : ENDDO
133 : ELSE
134 0 : DO i3=1,nup,2
135 0 : CALL cfft(a(1,1,i3-1),a(1,1,i3),(nw2+1)*n1,n1,n1,isn)
136 : CALL cfft(a(1,n2-nw2+1,i3-1),a(1,n2-nw2+1,i3), &
137 0 : nw2*n1,n1,n1,isn)
138 0 : CALL cfft(a(1,1,i3-1),a(1,1,i3),n1*n2,n2,n1*n2,isn)
139 : ENDDO
140 : ENDIF
141 :
142 : ! ******************** b a c k w a r d - t r a n s f o r m *******************
143 :
144 0 : IF (isn == -1) THEN
145 :
146 : ! the real part of f(0) has to be moved to a(n1,n2,1) for a correct
147 : ! setup for vrfftb (see comments therein)
148 :
149 0 : DO i1=1,n1
150 0 : DO i2=1,n2
151 0 : a(i1,i2,1)=a(i1,i2,0)
152 0 : a(i1,i2,0)=0.0
153 : ENDDO
154 : ENDDO
155 :
156 : ! transform n1*n2 mixed real and complex sequences of lenghth n3/2
157 : ! to n3 real values
158 :
159 0 : CALL vrffti(n3,wsave)
160 0 : IF (l_nopad) THEN
161 0 : CALL vrfftb(n1*n2,n3,a(1,1,1),b,n1d*n2d,wsave)
162 : ELSE
163 0 : DO i2=1,nw2+1
164 0 : CALL vrfftb(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
165 0 : CALL vrfftb(nw1 ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
166 : ENDDO
167 0 : DO i2=n2-nw2+1,n2
168 0 : CALL vrfftb(nw1+1,n3,a(1,i2,1),b,n1d*n2d,wsave)
169 0 : CALL vrfftb(nw1 ,n3,a(n1-nw1+1,i2,1),b,n1d*n2d,wsave)
170 : ENDDO
171 : ENDIF
172 :
173 : ! remove the norm 1/sqrt(n3) (compatibility with cfft)
174 :
175 0 : factor = sqrt(1.0*n3)
176 : ! ALL CPP_BLAS_sscal(n1d*n2d*n3,factor,a(1,1,1),1)
177 0 : a=a*factor
178 :
179 : ENDIF
180 :
181 0 : END SUBROUTINE rfft
182 :
183 0 : SUBROUTINE vrffti(n,wsave)
184 : !***BEGIN PROLOGUE VRFFTI
185 : !***DATE WRITTEN 860701 (YYMMDD)
186 : !***REVISION DATE 900509 (YYMMDD)
187 : !***CATEGORY NO. J1A1
188 : !***KEYWORDS FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
189 : ! MULTIPLE SEQUENCES
190 : !***AUTHOR SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
191 : !***PURPOSE Initialization for VRFFTF and VRFFTB.
192 : !***DESCRIPTION
193 :
194 : ! Subroutine VRFFTI initializes the array WSAVE which is used in
195 : ! both VRFFTF and VRFFTB. The prime factorization of N together with
196 : ! a tabulation of certain trigonometric functions are computed and
197 : ! stored in the array WSAVE.
198 :
199 : ! Input Parameter
200 :
201 : ! N the length of the sequence to be transformed. There is no
202 : ! restriction on N.
203 :
204 : ! Output Parameter
205 :
206 : ! WSAVE a work array which must be dimensioned at least N+15.
207 : ! The same work array can be used for both VRFFTF and VRFFTB
208 : ! as long as N remains unchanged. Different WSAVE arrays
209 : ! are required for different values of N. The contents of
210 : ! WSAVE must not be changed between calls of VRFFTF or VRFFTB.
211 :
212 : ! * * * * * * * * * * * * * * * * * * * * *
213 : ! * *
214 : ! * PROGRAM SPECIFICATIONS *
215 : ! * *
216 : ! * * * * * * * * * * * * * * * * * * * * *
217 :
218 : ! Dimension of R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
219 : ! Arguments
220 :
221 : ! Latest AUGUST 1, 1985
222 : ! Revision
223 :
224 : ! Subprograms VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
225 : ! Required VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
226 : ! VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
227 :
228 : ! Special NONE
229 : ! Conditions
230 :
231 : ! Common NONE
232 : ! blocks
233 :
234 : ! I/O NONE
235 :
236 : ! Precision SINGLE
237 :
238 : ! Specialist ROLAND SWEET
239 :
240 : ! Language FORTRAN
241 :
242 : ! History WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
243 : ! NATIONAL BUREAU OF STANDARDS (BOULDER).
244 :
245 : ! Algorithm A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
246 : ! OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
247 :
248 : ! Portability AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
249 : ! THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
250 : ! THE FUNCTION PIMACH.
251 :
252 : ! Required COS,SIN
253 : ! resident
254 : ! Routines
255 :
256 : !***REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
257 : ! Computations, (G. Rodrigue, ed.), Academic Press, 1982,
258 : ! pp. 51-83.
259 : !***ROUTINES CALLED VRFTI1
260 : !***END PROLOGUE VRFFTI
261 :
262 : ! VRFFTPK, VERSION 1, AUGUST 1985
263 :
264 : DIMENSION wsave(n+15)
265 : !***FIRST EXECUTABLE STATEMENT VRFFTI
266 0 : IF (n == 1) RETURN
267 0 : CALL vrfti1(n,wsave(1),wsave(n+1))
268 0 : RETURN
269 : END subroutine
270 0 : SUBROUTINE vrfti1(n,wa,fac)
271 :
272 : ! VRFFTPK, VERSION 1, AUGUST 1985
273 :
274 : USE m_constants, ONLY : pimach
275 : DIMENSION wa(n),fac(15),ntryh(4)
276 : DATA ntryh(1),ntryh(2),ntryh(3),ntryh(4)/4,2,3,5/
277 :
278 0 : nl = n
279 0 : nf = 0
280 0 : j = 0
281 0 : 10 j = j + 1
282 0 : IF ( j <=4 ) THEN
283 0 : ntry = ntryh(j)
284 : ELSE
285 0 : ntry = ntry + 2
286 : ENDIF
287 0 : 40 nq = nl/ntry
288 0 : nr = nl - ntry*nq
289 0 : IF ( nr /= 0 ) THEN
290 : GOTO 10
291 : ENDIF
292 0 : nf = nf + 1
293 0 : fac(nf+2) = ntry
294 0 : nl = nq
295 0 : IF (ntry /= 2) GO TO 70
296 0 : IF (nf == 1) GO TO 70
297 0 : DO 60 i = 2,nf
298 0 : ib = nf - i + 2
299 0 : fac(ib+2) = fac(ib+1)
300 0 : 60 END DO
301 0 : fac(3) = 2
302 0 : 70 IF (nl /= 1) GO TO 40
303 0 : fac(1) = n
304 0 : fac(2) = nf
305 0 : tpi = 2.*pimach()
306 0 : argh = tpi/real(n)
307 0 : is = 0
308 0 : nfm1 = nf - 1
309 0 : l1 = 1
310 0 : IF (nfm1 == 0) RETURN
311 0 : DO 100 k1 = 1,nfm1
312 0 : ip = fac(k1+2)
313 0 : ld = 0
314 0 : l2 = l1*ip
315 0 : ido = n/l2
316 0 : ipm = ip - 1
317 0 : DO 90 j = 1,ipm
318 0 : ld = ld + l1
319 0 : i = is
320 0 : argld = real(ld)*argh
321 0 : fi = 0.
322 0 : DO 80 ii = 3,ido,2
323 0 : i = i + 2
324 0 : fi = fi + 1.
325 0 : arg = fi*argld
326 0 : wa(i-1) = cos(arg)
327 0 : wa(i) = sin(arg)
328 0 : 80 END DO
329 0 : is = is + ido
330 0 : 90 END DO
331 0 : l1 = l2
332 0 : 100 END DO
333 : RETURN
334 : END subroutine
335 0 : SUBROUTINE vrfftf(m,n,r,rt,mdimr,wsave)
336 :
337 : !***BEGIN PROLOGUE VRFFTF
338 : !***DATE WRITTEN 850801 (YYMMDD)
339 : !***REVISION DATE 900509 (YYMMDD)
340 : !***CATEGORY NO. J1A1
341 : !***KEYWORDS FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
342 : ! FOURIER ANALYSIS, FORWARD TRANSFORM, MULTIPLE SEQUENCES
343 : !***AUTHOR SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
344 : !***PURPOSE Forward real periodic transform, M sequences.
345 : !***DESCRIPTION
346 :
347 : ! Subroutine VRFFTF computes the Fourier coefficients (forward
348 : ! transform) of a number of real periodic sequences. Specifically,
349 : ! for each sequence the subroutine claculates the independent
350 : ! Fourier coefficients described below at output parameter R.
351 :
352 : ! The array WSAVE which is used by subroutine VRFFTF must be
353 : ! initialized by calling subroutine VRFFTI(N,WSAVE).
354 :
355 : ! Input Parameters
356 :
357 : ! M the number of sequences to be transformed.
358 :
359 : ! N the length of the sequences to be transformed. The method
360 : ! is most efficient when N is a product of small primes,
361 : ! however n may be any positive integer.
362 :
363 : ! R areal two-dimensional array of size MDIMX x N containing the
364 : ! the sequences to be transformed. The sequences are stored
365 : ! in the ROWS of R. Thus, the I-th sequence to be transformed,
366 : ! X(I,J), J=0,1,...,N-1, is stored as
367 :
368 : ! R(I,J) = X(I,J-1) , J=1, 2, . . . , N.
369 :
370 : ! RT a real two-dimensional work array of size MDIMX x N.
371 :
372 : ! MDIMR the row (or first) dimension of the arrays R and RT exactly
373 : ! as they appear in the calling program. This parameter is
374 : ! used to specify the variable dimension of these arrays.
375 :
376 : ! WSAVE a real one-dimensional work array which must be dimensioned
377 : ! at least N+15. The WSAVE array must be initialized by
378 : ! calling subroutine VRFFTI. A different WSAVE array must be
379 : ! used for each different value of N. This initialization does
380 : ! not have to be repeated so long as N remains unchanged. The
381 : ! same WSAVE array may be used by VRFFTF and VRFFTB.
382 :
383 : ! Output Parameters
384 :
385 : ! R contains the Fourier coefficients F(K) for each of the M
386 : ! input sequences. Specifically, row I of R, R(I,J),
387 : ! J=1,2,..,N, contains the independent Fourier coefficients
388 : ! F(I,K), for the I-th input sequence stored as
389 :
390 : ! R(I,1) = REAL( F(I,0) ),
391 : ! = SQRT(1/N)*SUM(J=0,N-1)[ X(I,J) ],
392 :
393 : ! R(I,2*K) = REAL( F(I,K) )
394 : ! = SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*COS(2J*K*PI/N)]
395 :
396 : ! R(I,2*K+1) = IMAG( F(I,K) )
397 : ! =-SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*SIN(2J*K*PI/N)]
398 :
399 : ! for K = 1, 2, . . . , M-1,
400 :
401 : ! and, when N is even,
402 :
403 : ! R(I,N) = REAL( F(I,N/2) ).
404 : ! = SQRT(1/N)*SUM(J=0,N-1)[ (-1)**J*X(I,J) ].
405 :
406 : ! WSAVE contains results which must not be destroyed between calls
407 : ! to VRFFTF or VRFFTB.
408 :
409 : ! -----------------------------------------------------------------
410 :
411 : ! NOTE - A call of VRFFTF followed immediately by a call of
412 : ! of VRFFTB will return the original sequences R. Thus,
413 : ! VRFFTB is the correctly normalized inverse of VRFFTF.
414 :
415 : ! -----------------------------------------------------------------
416 :
417 : ! VRFFTF is a straightforward extension of the subprogram RFFTF to
418 : ! handle M simultaneous sequences. RFFTF was originally developed
419 : ! by P. N. Swarztrauber of NCAR.
420 :
421 : ! * * * * * * * * * * * * * * * * * * * * *
422 : ! * *
423 : ! * PROGRAM SPECIFICATIONS *
424 : ! * *
425 : ! * * * * * * * * * * * * * * * * * * * * *
426 :
427 : ! Dimension of R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
428 : ! Arguments
429 :
430 : ! Latest AUGUST 1, 1985
431 : ! Revision
432 :
433 : ! Subprograms VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
434 : ! Required VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
435 : ! VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
436 :
437 : ! Special NONE
438 : ! Conditions
439 :
440 : ! Common NONE
441 : ! blocks
442 :
443 : ! I/O NONE
444 :
445 : ! Precision SINGLE
446 :
447 : ! Specialist ROLAND SWEET
448 :
449 : ! Language FORTRAN
450 :
451 : ! History WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
452 : ! NATIONAL BUREAU OF STANDARDS (BOULDER).
453 :
454 : ! Algorithm A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
455 : ! OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
456 :
457 : ! Portability AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
458 : ! THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
459 : ! THE FUNCTION PIMACH.
460 :
461 : ! Required COS,SIN
462 : ! resident
463 : ! Routines
464 :
465 : !***REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
466 : ! Computations, (G. Rodrigue, ed.), Academic Press, 1982,
467 : ! pp. 51-83.
468 : !***ROUTINES CALLED VRFTF1
469 : !***END PROLOGUE VRFFTF
470 :
471 : ! VRFFTPK, VERSION 1, AUGUST 1985
472 :
473 : DIMENSION r(mdimr,n),rt(mdimr,n),wsave(n+15)
474 : !***FIRST EXECUTABLE STATEMENT VRFFTF
475 0 : IF (n == 1) RETURN
476 0 : CALL vrftf1(m,n,r,rt,mdimr,wsave(1),wsave(n+1))
477 0 : RETURN
478 : END subroutine
479 0 : SUBROUTINE vradf2(mp,ido,l1,cc,ch,mdimc,wa1)
480 :
481 : ! VRFFTPK, VERSION 1, AUGUST 1985
482 :
483 : DIMENSION ch(mdimc,ido,2,l1),cc(mdimc,ido,l1,2),wa1(ido)
484 :
485 0 : DO 20 k = 1,l1
486 0 : DO 10 m = 1,mp
487 0 : ch(m,1,1,k) = cc(m,1,k,1) + cc(m,1,k,2)
488 0 : ch(m,ido,2,k) = cc(m,1,k,1) - cc(m,1,k,2)
489 0 : 10 END DO
490 0 : 20 END DO
491 0 : IF ( ido == 2 ) THEN
492 : GOTO 70
493 0 : ELSEIF ( ido < 2 ) THEN
494 : GOTO 100
495 : ENDIF
496 0 : idp2 = ido + 2
497 0 : DO 60 k = 1,l1
498 0 : DO 50 i = 3,ido,2
499 0 : ic = idp2 - i
500 0 : DO 40 m = 1,mp
501 : ch(m,i,1,k) = cc(m,i,k,1) + &
502 : (wa1(i-2)*cc(m,i,k,2)-wa1(i-1)* &
503 0 : cc(m,i-1,k,2))
504 : ch(m,ic,2,k) = (wa1(i-2)*cc(m,i,k,2)- &
505 0 : wa1(i-1)*cc(m,i-1,k,2)) - cc(m,i,k,1)
506 : ch(m,i-1,1,k) = cc(m,i-1,k,1) + &
507 : (wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)* &
508 0 : cc(m,i,k,2))
509 : ch(m,ic-1,2,k) = cc(m,i-1,k,1) - &
510 : (wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)* &
511 0 : cc(m,i,k,2))
512 0 : 40 END DO
513 0 : 50 END DO
514 0 : 60 END DO
515 0 : IF (mod(ido,2) == 1) RETURN
516 0 : 70 DO 90 k = 1,l1
517 0 : DO 80 m = 1,mp
518 0 : ch(m,1,2,k) = -cc(m,ido,k,2)
519 0 : ch(m,ido,1,k) = cc(m,ido,k,1)
520 0 : 80 END DO
521 0 : 90 END DO
522 : 100 RETURN
523 : END subroutine
524 0 : SUBROUTINE vradf3(mp,ido,l1,cc,ch,mdimc,wa1,wa2)
525 :
526 : ! VRFFTPK, VERSION 1, AUGUST 1985
527 :
528 : USE m_constants, ONLY : pimach
529 : DIMENSION ch(mdimc,ido,3,l1),cc(mdimc,ido,l1,3),wa1(ido),wa2(ido)
530 :
531 0 : arg = 2.*pimach()/3.
532 0 : taur = cos(arg)
533 0 : taui = sin(arg)
534 0 : DO 20 k = 1,l1
535 0 : DO 10 m = 1,mp
536 0 : ch(m,1,1,k) = cc(m,1,k,1) + (cc(m,1,k,2)+cc(m,1,k,3))
537 0 : ch(m,1,3,k) = taui* (cc(m,1,k,3)-cc(m,1,k,2))
538 : ch(m,ido,2,k) = cc(m,1,k,1) + &
539 0 : taur* (cc(m,1,k,2)+cc(m,1,k,3))
540 0 : 10 END DO
541 0 : 20 END DO
542 0 : IF (ido == 1) RETURN
543 0 : idp2 = ido + 2
544 0 : DO 50 k = 1,l1
545 0 : DO 40 i = 3,ido,2
546 0 : ic = idp2 - i
547 0 : DO 30 m = 1,mp
548 : ch(m,i-1,1,k) = cc(m,i-1,k,1) + &
549 : ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i, &
550 : k,2))+ (wa2(i-2)*cc(m,i-1,k, &
551 0 : & 3)+wa2(i-1)*cc(m,i,k,3)))
552 : ch(m,i,1,k) = cc(m,i,k,1) + &
553 : ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
554 : & 2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m, &
555 0 : i-1,k,3)))
556 : ch(m,i-1,3,k) = (cc(m,i-1,k,1)+ &
557 : taur* ((wa1(i-2)*cc(m,i-1,k, &
558 : & 2)+wa1(i-1)*cc(m,i,k,2))+ (wa2(i-2)*cc(m, &
559 : i-1,k,3)+wa2(i-1)*cc(m,i,k,3)))) + &
560 : (taui* ((wa1(i-2)*cc(m,i,k, &
561 : & 2)-wa1(i-1)*cc(m,i-1,k, &
562 : & 2))- (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m, &
563 0 : i-1,k,3))))
564 : ch(m,ic-1,2,k) = (cc(m,i-1,k,1)+ &
565 : taur* ((wa1(i-2)*cc(m,i-1,k, &
566 : & 2)+wa1(i-1)*cc(m,i,k, &
567 : & 2))+ (wa2(i-2)*cc(m,i-1,k, &
568 : & 3)+wa2(i-1)*cc(m,i,k,3)))) - &
569 : (taui* ((wa1(i-2)*cc(m,i,k, &
570 : & 2)-wa1(i-1)*cc(m,i-1,k, &
571 : & 2))- (wa2(i-2)*cc(m,i,k, &
572 0 : & 3)-wa2(i-1)*cc(m,i-1,k,3))))
573 : ch(m,i,3,k) = (cc(m,i,k,1)+taur* &
574 : ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
575 : & 2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m, &
576 : i-1,k,3)))) + (taui* ((wa2(i-2)*cc(m,i-1,k, &
577 : & 3)+wa2(i-1)*cc(m,i,k,3))- (wa1(i-2)*cc(m, &
578 0 : i-1,k,2)+wa1(i-1)*cc(m,i,k,2))))
579 : ch(m,ic,2,k) = (taui* ((wa2(i-2)*cc(m,i-1,k, &
580 : & 3)+wa2(i-1)*cc(m,i,k,3))- (wa1(i-2)*cc(m, &
581 : i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))) - &
582 : (cc(m,i,k,1)+taur* ((wa1(i-2)*cc(m,i,k, &
583 : & 2)-wa1(i-1)*cc(m,i-1,k, &
584 : & 2))+ (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m, &
585 0 : i-1,k,3))))
586 0 : 30 END DO
587 0 : 40 END DO
588 0 : 50 END DO
589 : RETURN
590 : END subroutine
591 0 : SUBROUTINE vradf4(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3)
592 :
593 : ! VRFFTPK, VERSION 1, AUGUST 1985
594 :
595 : DIMENSION cc(mdimc,ido,l1,4),ch(mdimc,ido,4,l1),wa1(ido),wa2(ido), &
596 : wa3(ido)
597 :
598 0 : hsqt2 = sqrt(2.)/2.
599 0 : DO 20 k = 1,l1
600 0 : DO 10 m = 1,mp
601 : ch(m,1,1,k) = (cc(m,1,k,2)+cc(m,1,k,4)) + &
602 0 : (cc(m,1,k,1)+cc(m,1,k,3))
603 : ch(m,ido,4,k) = (cc(m,1,k,1)+cc(m,1,k,3)) - &
604 0 : (cc(m,1,k,2)+cc(m,1,k,4))
605 0 : ch(m,ido,2,k) = cc(m,1,k,1) - cc(m,1,k,3)
606 0 : ch(m,1,3,k) = cc(m,1,k,4) - cc(m,1,k,2)
607 0 : 10 END DO
608 0 : 20 END DO
609 0 : IF ( ido == 2 ) THEN
610 : GOTO 70
611 0 : ELSEIF ( ido < 2 ) THEN
612 : GOTO 100
613 : ENDIF
614 0 : idp2 = ido + 2
615 0 : DO 60 k = 1,l1
616 0 : DO 50 i = 3,ido,2
617 0 : ic = idp2 - i
618 0 : DO 40 m = 1,mp
619 : ch(m,i-1,1,k) = ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i, &
620 : k,2))+ (wa3(i-2)*cc(m,i-1,k, &
621 : & 4)+wa3(i-1)*cc(m,i,k,4))) + &
622 : (cc(m,i-1,k,1)+ (wa2(i-2)*cc(m,i-1,k, &
623 0 : & 3)+wa2(i-1)*cc(m,i,k,3)))
624 : ch(m,ic-1,4,k) = (cc(m,i-1,k,1)+ &
625 : (wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i, &
626 : k,3))) - ((wa1(i-2)*cc(m,i-1,k, &
627 : & 2)+wa1(i-1)*cc(m,i,k,2))+ &
628 : (wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i, &
629 0 : k,4)))
630 : ch(m,i,1,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
631 : & 2))+ (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
632 : i-1,k,4))) + (cc(m,i,k,1)+ &
633 : (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k, &
634 0 : & 3)))
635 : ch(m,ic,4,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1, &
636 : k,2))+ (wa3(i-2)*cc(m,i,k, &
637 : & 4)-wa3(i-1)*cc(m,i-1,k,4))) - &
638 : (cc(m,i,k,1)+ (wa2(i-2)*cc(m,i,k, &
639 0 : & 3)-wa2(i-1)*cc(m,i-1,k,3)))
640 : ch(m,i-1,3,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1, &
641 : k,2))- (wa3(i-2)*cc(m,i,k, &
642 : & 4)-wa3(i-1)*cc(m,i-1,k,4))) + &
643 : (cc(m,i-1,k,1)- (wa2(i-2)*cc(m,i-1,k, &
644 0 : & 3)+wa2(i-1)*cc(m,i,k,3)))
645 : ch(m,ic-1,2,k) = (cc(m,i-1,k,1)- &
646 : (wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i, &
647 : k,3))) - ((wa1(i-2)*cc(m,i,k, &
648 : & 2)-wa1(i-1)*cc(m,i-1,k,2))- &
649 : (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,i-1, &
650 0 : k,4)))
651 : ch(m,i,3,k) = ((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i,k, &
652 : & 4))- (wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m, &
653 : i,k,2))) + (cc(m,i,k,1)- &
654 : (wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k, &
655 0 : & 3)))
656 : ch(m,ic,2,k) = ((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i, &
657 : k,4))- (wa1(i-2)*cc(m,i-1,k, &
658 : & 2)+wa1(i-1)*cc(m,i,k,2))) - &
659 : (cc(m,i,k,1)- (wa2(i-2)*cc(m,i,k, &
660 0 : & 3)-wa2(i-1)*cc(m,i-1,k,3)))
661 0 : 40 END DO
662 0 : 50 END DO
663 0 : 60 END DO
664 0 : IF (mod(ido,2) == 1) RETURN
665 : 70 CONTINUE
666 0 : DO 90 k = 1,l1
667 0 : DO 80 m = 1,mp
668 : ch(m,ido,1,k) = (hsqt2* (cc(m,ido,k,2)-cc(m,ido,k,4))) + &
669 0 : cc(m,ido,k,1)
670 : ch(m,ido,3,k) = cc(m,ido,k,1) - &
671 0 : (hsqt2* (cc(m,ido,k,2)-cc(m,ido,k,4)))
672 : ch(m,1,2,k) = (-hsqt2* (cc(m,ido,k,2)+cc(m,ido,k,4))) - &
673 0 : cc(m,ido,k,3)
674 : ch(m,1,4,k) = (-hsqt2* (cc(m,ido,k,2)+cc(m,ido,k,4))) + &
675 0 : cc(m,ido,k,3)
676 0 : 80 END DO
677 0 : 90 END DO
678 : 100 RETURN
679 : END subroutine
680 0 : SUBROUTINE vradf5(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3,wa4)
681 :
682 : ! VRFFTPK, VERSION 1, AUGUST 1985
683 :
684 : USE m_constants, ONLY : pimach
685 : DIMENSION cc(mdimc,ido,l1,5),ch(mdimc,ido,5,l1),wa1(ido),wa2(ido), &
686 : wa3(ido),wa4(ido)
687 :
688 0 : arg = 2.*pimach()/5.
689 0 : tr11 = cos(arg)
690 0 : ti11 = sin(arg)
691 0 : tr12 = cos(2.*arg)
692 0 : ti12 = sin(2.*arg)
693 0 : DO 20 k = 1,l1
694 0 : DO 10 m = 1,mp
695 : ch(m,1,1,k) = cc(m,1,k,1) + (cc(m,1,k,5)+cc(m,1,k,2)) + &
696 0 : (cc(m,1,k,4)+cc(m,1,k,3))
697 : ch(m,ido,2,k) = cc(m,1,k,1) + &
698 : tr11* (cc(m,1,k,5)+cc(m,1,k,2)) + &
699 0 : tr12* (cc(m,1,k,4)+cc(m,1,k,3))
700 : ch(m,1,3,k) = ti11* (cc(m,1,k,5)-cc(m,1,k,2)) + &
701 0 : ti12* (cc(m,1,k,4)-cc(m,1,k,3))
702 : ch(m,ido,4,k) = cc(m,1,k,1) + &
703 : tr12* (cc(m,1,k,5)+cc(m,1,k,2)) + &
704 0 : tr11* (cc(m,1,k,4)+cc(m,1,k,3))
705 : ch(m,1,5,k) = ti12* (cc(m,1,k,5)-cc(m,1,k,2)) - &
706 0 : ti11* (cc(m,1,k,4)-cc(m,1,k,3))
707 0 : 10 END DO
708 0 : 20 END DO
709 0 : IF (ido == 1) RETURN
710 0 : idp2 = ido + 2
711 0 : DO 50 k = 1,l1
712 0 : DO 40 i = 3,ido,2
713 0 : ic = idp2 - i
714 0 : DO 30 m = 1,mp
715 : ch(m,i-1,1,k) = cc(m,i-1,k,1) + &
716 : ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i, &
717 : k,2))+ (wa4(i-2)*cc(m,i-1,k, &
718 : & 5)+wa4(i-1)*cc(m,i,k,5))) + &
719 : ((wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i, &
720 : k,3))+ (wa3(i-2)*cc(m,i-1,k, &
721 0 : & 4)+wa3(i-1)*cc(m,i,k,4)))
722 : ch(m,i,1,k) = cc(m,i,k,1) + &
723 : ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
724 : & 2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
725 : i-1,k,5))) + ((wa2(i-2)*cc(m,i,k, &
726 : & 3)-wa2(i-1)*cc(m,i-1,k,3))+ &
727 : (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,i-1,k, &
728 0 : & 4)))
729 : ch(m,i-1,3,k) = cc(m,i-1,k,1) + &
730 : tr11* (wa1(i-2)*cc(m,i-1,k,2)+ &
731 : wa1(i-1)*cc(m,i,k,2)+ &
732 : wa4(i-2)*cc(m,i-1,k,5)+ &
733 : wa4(i-1)*cc(m,i,k,5)) + &
734 : tr12* (wa2(i-2)*cc(m,i-1,k,3)+ &
735 : wa2(i-1)*cc(m,i,k,3)+ &
736 : wa3(i-2)*cc(m,i-1,k,4)+ &
737 : wa3(i-1)*cc(m,i,k,4)) + &
738 : ti11* (wa1(i-2)*cc(m,i,k,2)- &
739 : wa1(i-1)*cc(m,i-1,k,2)- &
740 : (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m,i-1, &
741 : k,5))) + ti12* (wa2(i-2)*cc(m,i,k,3)- &
742 : wa2(i-1)*cc(m,i-1,k,3)- &
743 : (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,i-1, &
744 0 : k,4)))
745 : ch(m,ic-1,2,k) = cc(m,i-1,k,1) + &
746 : tr11* (wa1(i-2)*cc(m,i-1,k,2)+ &
747 : wa1(i-1)*cc(m,i,k,2)+ &
748 : wa4(i-2)*cc(m,i-1,k,5)+ &
749 : wa4(i-1)*cc(m,i,k,5)) + &
750 : tr12* (wa2(i-2)*cc(m,i-1,k,3)+ &
751 : wa2(i-1)*cc(m,i,k,3)+ &
752 : wa3(i-2)*cc(m,i-1,k,4)+ &
753 : wa3(i-1)*cc(m,i,k,4)) - &
754 : (ti11* (wa1(i-2)*cc(m,i,k, &
755 : & 2)-wa1(i-1)*cc(m,i-1,k, &
756 : & 2)- (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
757 : i-1,k,5)))+ti12* (wa2(i-2)*cc(m,i,k, &
758 : & 3)-wa2(i-1)*cc(m,i-1,k, &
759 : & 3)- (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
760 0 : i-1,k,4))))
761 : ch(m,i,3,k) = (cc(m,i,k,1)+tr11* &
762 : ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
763 : & 2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
764 : i-1,k,5)))+tr12* ((wa2(i-2)*cc(m,i,k, &
765 : & 3)-wa2(i-1)*cc(m,i-1,k,3))+ (wa3(i-2)*cc(m, &
766 : i,k,4)-wa3(i-1)*cc(m,i-1,k,4)))) + &
767 : (ti11* ((wa4(i-2)*cc(m,i-1,k, &
768 : & 5)+wa4(i-1)*cc(m,i,k,5))- (wa1(i-2)*cc(m, &
769 : i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))+ &
770 : ti12* ((wa3(i-2)*cc(m,i-1,k, &
771 : & 4)+wa3(i-1)*cc(m,i,k,4))- (wa2(i-2)*cc(m, &
772 0 : i-1,k,3)+wa2(i-1)*cc(m,i,k,3))))
773 : ch(m,ic,2,k) = (ti11* ((wa4(i-2)*cc(m,i-1,k, &
774 : & 5)+wa4(i-1)*cc(m,i,k,5))- (wa1(i-2)*cc(m, &
775 : i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))+ &
776 : ti12* ((wa3(i-2)*cc(m,i-1,k, &
777 : & 4)+wa3(i-1)*cc(m,i,k,4))- (wa2(i-2)*cc(m, &
778 : i-1,k,3)+wa2(i-1)*cc(m,i,k,3)))) - &
779 : (cc(m,i,k,1)+tr11* ((wa1(i-2)*cc(m,i,k, &
780 : & 2)-wa1(i-1)*cc(m,i-1,k, &
781 : & 2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
782 : i-1,k,5)))+tr12* ((wa2(i-2)*cc(m,i,k, &
783 : & 3)-wa2(i-1)*cc(m,i-1,k, &
784 : & 3))+ (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
785 0 : i-1,k,4))))
786 : ch(m,i-1,5,k) = (cc(m,i-1,k,1)+ &
787 : tr12* ((wa1(i-2)*cc(m,i-1,k, &
788 : & 2)+wa1(i-1)*cc(m,i,k,2))+ (wa4(i-2)*cc(m, &
789 : i-1,k,5)+wa4(i-1)*cc(m,i,k,5)))+ &
790 : tr11* ((wa2(i-2)*cc(m,i-1,k, &
791 : & 3)+wa2(i-1)*cc(m,i,k,3))+ (wa3(i-2)*cc(m, &
792 : i-1,k,4)+wa3(i-1)*cc(m,i,k,4)))) + &
793 : (ti12* ((wa1(i-2)*cc(m,i,k, &
794 : & 2)-wa1(i-1)*cc(m,i-1,k, &
795 : & 2))- (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
796 : i-1,k,5)))-ti11* ((wa2(i-2)*cc(m,i,k, &
797 : & 3)-wa2(i-1)*cc(m,i-1,k, &
798 : & 3))- (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
799 0 : i-1,k,4))))
800 : ch(m,ic-1,4,k) = (cc(m,i-1,k,1)+ &
801 : tr12* ((wa1(i-2)*cc(m,i-1,k, &
802 : & 2)+wa1(i-1)*cc(m,i,k, &
803 : & 2))+ (wa4(i-2)*cc(m,i-1,k, &
804 : & 5)+wa4(i-1)*cc(m,i,k,5)))+ &
805 : tr11* ((wa2(i-2)*cc(m,i-1,k, &
806 : & 3)+wa2(i-1)*cc(m,i,k, &
807 : & 3))+ (wa3(i-2)*cc(m,i-1,k, &
808 : & 4)+wa3(i-1)*cc(m,i,k,4)))) - &
809 : (ti12* ((wa1(i-2)*cc(m,i,k, &
810 : & 2)-wa1(i-1)*cc(m,i-1,k, &
811 : & 2))- (wa4(i-2)*cc(m,i,k, &
812 : & 5)-wa4(i-1)*cc(m,i-1,k,5)))- &
813 : ti11* ((wa2(i-2)*cc(m,i,k, &
814 : & 3)-wa2(i-1)*cc(m,i-1,k, &
815 : & 3))- (wa3(i-2)*cc(m,i,k, &
816 0 : & 4)-wa3(i-1)*cc(m,i-1,k,4))))
817 : ch(m,i,5,k) = (cc(m,i,k,1)+tr12* &
818 : ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k, &
819 : & 2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
820 : i-1,k,5)))+tr11* ((wa2(i-2)*cc(m,i,k, &
821 : & 3)-wa2(i-1)*cc(m,i-1,k,3))+ (wa3(i-2)*cc(m, &
822 : i,k,4)-wa3(i-1)*cc(m,i-1,k,4)))) + &
823 : (ti12* ((wa4(i-2)*cc(m,i-1,k, &
824 : & 5)+wa4(i-1)*cc(m,i,k,5))- (wa1(i-2)*cc(m, &
825 : i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))- &
826 : ti11* ((wa3(i-2)*cc(m,i-1,k, &
827 : & 4)+wa3(i-1)*cc(m,i,k,4))- (wa2(i-2)*cc(m, &
828 0 : i-1,k,3)+wa2(i-1)*cc(m,i,k,3))))
829 : ch(m,ic,4,k) = (ti12* ((wa4(i-2)*cc(m,i-1,k, &
830 : & 5)+wa4(i-1)*cc(m,i,k,5))- (wa1(i-2)*cc(m, &
831 : i-1,k,2)+wa1(i-1)*cc(m,i,k,2)))- &
832 : ti11* ((wa3(i-2)*cc(m,i-1,k, &
833 : & 4)+wa3(i-1)*cc(m,i,k,4))- (wa2(i-2)*cc(m, &
834 : i-1,k,3)+wa2(i-1)*cc(m,i,k,3)))) - &
835 : (cc(m,i,k,1)+tr12* ((wa1(i-2)*cc(m,i,k, &
836 : & 2)-wa1(i-1)*cc(m,i-1,k, &
837 : & 2))+ (wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m, &
838 : i-1,k,5)))+tr11* ((wa2(i-2)*cc(m,i,k, &
839 : & 3)-wa2(i-1)*cc(m,i-1,k, &
840 : & 3))+ (wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m, &
841 0 : i-1,k,4))))
842 0 : 30 END DO
843 0 : 40 END DO
844 0 : 50 END DO
845 : RETURN
846 : END subroutine
847 0 : SUBROUTINE vradfg(mp,ido,ip,l1,idl1,cc,c1,c2,ch,ch2,mdimc,wa)
848 :
849 : ! VRFFTPK, VERSION 1, AUGUST 1985
850 :
851 : USE m_constants, ONLY : pimach
852 : DIMENSION ch(mdimc,ido,l1,ip),cc(mdimc,ido,ip,l1), &
853 : c1(mdimc,ido,l1,ip),c2(mdimc,idl1,ip), &
854 : ch2(mdimc,idl1,ip),wa(ido)
855 :
856 0 : tpi = 2.*pimach()
857 0 : arg = tpi/real(ip)
858 0 : dcp = cos(arg)
859 0 : dsp = sin(arg)
860 0 : ipph = (ip+1)/2
861 0 : ipp2 = ip + 2
862 0 : idp2 = ido + 2
863 0 : nbd = (ido-1)/2
864 0 : IF (ido == 1) GO TO 250
865 0 : DO 20 ik = 1,idl1
866 0 : DO 10 m = 1,mp
867 0 : ch2(m,ik,1) = c2(m,ik,1)
868 0 : 10 END DO
869 0 : 20 END DO
870 0 : DO 50 j = 2,ip
871 0 : DO 40 k = 1,l1
872 0 : DO 30 m = 1,mp
873 0 : ch(m,1,k,j) = c1(m,1,k,j)
874 0 : 30 END DO
875 0 : 40 END DO
876 0 : 50 END DO
877 0 : IF (nbd > l1) GO TO 100
878 0 : is = -ido
879 0 : DO 90 j = 2,ip
880 0 : is = is + ido
881 0 : idij = is
882 0 : DO 80 i = 3,ido,2
883 0 : idij = idij + 2
884 0 : DO 70 k = 1,l1
885 0 : DO 60 m = 1,mp
886 : ch(m,i-1,k,j) = wa(idij-1)*c1(m,i-1,k,j) + &
887 0 : wa(idij)*c1(m,i,k,j)
888 : ch(m,i,k,j) = wa(idij-1)*c1(m,i,k,j) - &
889 0 : wa(idij)*c1(m,i-1,k,j)
890 0 : 60 END DO
891 0 : 70 END DO
892 0 : 80 END DO
893 0 : 90 END DO
894 0 : GO TO 150
895 0 : 100 is = -ido
896 0 : DO 140 j = 2,ip
897 0 : is = is + ido
898 0 : DO 130 k = 1,l1
899 0 : idij = is
900 0 : DO 120 i = 3,ido,2
901 0 : idij = idij + 2
902 0 : DO 110 m = 1,mp
903 : ch(m,i-1,k,j) = wa(idij-1)*c1(m,i-1,k,j) + &
904 0 : wa(idij)*c1(m,i,k,j)
905 : ch(m,i,k,j) = wa(idij-1)*c1(m,i,k,j) - &
906 0 : wa(idij)*c1(m,i-1,k,j)
907 0 : 110 END DO
908 0 : 120 END DO
909 0 : 130 END DO
910 0 : 140 END DO
911 0 : 150 IF (nbd < l1) GO TO 200
912 0 : DO 190 j = 2,ipph
913 0 : jc = ipp2 - j
914 0 : DO 180 k = 1,l1
915 0 : DO 170 i = 3,ido,2
916 0 : DO 160 m = 1,mp
917 0 : c1(m,i-1,k,j) = ch(m,i-1,k,j) + ch(m,i-1,k,jc)
918 0 : c1(m,i-1,k,jc) = ch(m,i,k,j) - ch(m,i,k,jc)
919 0 : c1(m,i,k,j) = ch(m,i,k,j) + ch(m,i,k,jc)
920 0 : c1(m,i,k,jc) = ch(m,i-1,k,jc) - ch(m,i-1,k,j)
921 0 : 160 END DO
922 0 : 170 END DO
923 0 : 180 END DO
924 0 : 190 END DO
925 0 : GO TO 280
926 0 : 200 DO 240 j = 2,ipph
927 0 : jc = ipp2 - j
928 0 : DO 230 i = 3,ido,2
929 0 : DO 220 k = 1,l1
930 0 : DO 210 m = 1,mp
931 0 : c1(m,i-1,k,j) = ch(m,i-1,k,j) + ch(m,i-1,k,jc)
932 0 : c1(m,i-1,k,jc) = ch(m,i,k,j) - ch(m,i,k,jc)
933 0 : c1(m,i,k,j) = ch(m,i,k,j) + ch(m,i,k,jc)
934 0 : c1(m,i,k,jc) = ch(m,i-1,k,jc) - ch(m,i-1,k,j)
935 0 : 210 END DO
936 0 : 220 END DO
937 0 : 230 END DO
938 0 : 240 END DO
939 0 : GO TO 280
940 0 : 250 DO 270 ik = 1,idl1
941 0 : DO 260 m = 1,mp
942 0 : c2(m,ik,1) = ch2(m,ik,1)
943 0 : 260 END DO
944 0 : 270 END DO
945 0 : 280 DO 310 j = 2,ipph
946 0 : jc = ipp2 - j
947 0 : DO 300 k = 1,l1
948 0 : DO 290 m = 1,mp
949 0 : c1(m,1,k,j) = ch(m,1,k,j) + ch(m,1,k,jc)
950 0 : c1(m,1,k,jc) = ch(m,1,k,jc) - ch(m,1,k,j)
951 0 : 290 END DO
952 0 : 300 END DO
953 0 : 310 END DO
954 :
955 : ar1 = 1.
956 : ai1 = 0.
957 0 : DO 370 l = 2,ipph
958 0 : lc = ipp2 - l
959 0 : ar1h = dcp*ar1 - dsp*ai1
960 0 : ai1 = dcp*ai1 + dsp*ar1
961 0 : ar1 = ar1h
962 0 : DO 330 ik = 1,idl1
963 0 : DO 320 m = 1,mp
964 0 : ch2(m,ik,l) = c2(m,ik,1) + ar1*c2(m,ik,2)
965 0 : ch2(m,ik,lc) = ai1*c2(m,ik,ip)
966 0 : 320 END DO
967 0 : 330 END DO
968 : dc2 = ar1
969 : ds2 = ai1
970 : ar2 = ar1
971 : ai2 = ai1
972 0 : DO 360 j = 3,ipph
973 0 : jc = ipp2 - j
974 0 : ar2h = dc2*ar2 - ds2*ai2
975 0 : ai2 = dc2*ai2 + ds2*ar2
976 0 : ar2 = ar2h
977 0 : DO 350 ik = 1,idl1
978 0 : DO 340 m = 1,mp
979 0 : ch2(m,ik,l) = ch2(m,ik,l) + ar2*c2(m,ik,j)
980 0 : ch2(m,ik,lc) = ch2(m,ik,lc) + ai2*c2(m,ik,jc)
981 0 : 340 END DO
982 0 : 350 END DO
983 0 : 360 END DO
984 0 : 370 END DO
985 0 : DO 400 j = 2,ipph
986 0 : DO 390 ik = 1,idl1
987 0 : DO 380 m = 1,mp
988 0 : ch2(m,ik,1) = ch2(m,ik,1) + c2(m,ik,j)
989 0 : 380 END DO
990 0 : 390 END DO
991 0 : 400 END DO
992 :
993 0 : IF (ido < l1) GO TO 440
994 0 : DO 430 k = 1,l1
995 0 : DO 420 i = 1,ido
996 0 : DO 410 m = 1,mp
997 0 : cc(m,i,1,k) = ch(m,i,k,1)
998 0 : 410 END DO
999 0 : 420 END DO
1000 0 : 430 END DO
1001 0 : GO TO 480
1002 0 : 440 DO 470 i = 1,ido
1003 0 : DO 460 k = 1,l1
1004 0 : DO 450 m = 1,mp
1005 0 : cc(m,i,1,k) = ch(m,i,k,1)
1006 0 : 450 END DO
1007 0 : 460 END DO
1008 0 : 470 END DO
1009 0 : 480 DO 510 j = 2,ipph
1010 0 : jc = ipp2 - j
1011 0 : j2 = j + j
1012 0 : DO 500 k = 1,l1
1013 0 : DO 490 m = 1,mp
1014 0 : cc(m,ido,j2-2,k) = ch(m,1,k,j)
1015 0 : cc(m,1,j2-1,k) = ch(m,1,k,jc)
1016 0 : 490 END DO
1017 0 : 500 END DO
1018 0 : 510 END DO
1019 0 : IF (ido == 1) RETURN
1020 0 : IF (nbd < l1) GO TO 560
1021 0 : DO 550 j = 2,ipph
1022 0 : jc = ipp2 - j
1023 0 : j2 = j + j
1024 0 : DO 540 k = 1,l1
1025 0 : DO 530 i = 3,ido,2
1026 0 : ic = idp2 - i
1027 0 : DO 520 m = 1,mp
1028 0 : cc(m,i-1,j2-1,k) = ch(m,i-1,k,j) + ch(m,i-1,k,jc)
1029 0 : cc(m,ic-1,j2-2,k) = ch(m,i-1,k,j) - ch(m,i-1,k,jc)
1030 0 : cc(m,i,j2-1,k) = ch(m,i,k,j) + ch(m,i,k,jc)
1031 0 : cc(m,ic,j2-2,k) = ch(m,i,k,jc) - ch(m,i,k,j)
1032 0 : 520 END DO
1033 0 : 530 END DO
1034 0 : 540 END DO
1035 0 : 550 END DO
1036 0 : RETURN
1037 0 : 560 DO 600 j = 2,ipph
1038 0 : jc = ipp2 - j
1039 0 : j2 = j + j
1040 0 : DO 590 i = 3,ido,2
1041 0 : ic = idp2 - i
1042 0 : DO 580 k = 1,l1
1043 0 : DO 570 m = 1,mp
1044 0 : cc(m,i-1,j2-1,k) = ch(m,i-1,k,j) + ch(m,i-1,k,jc)
1045 0 : cc(m,ic-1,j2-2,k) = ch(m,i-1,k,j) - ch(m,i-1,k,jc)
1046 0 : cc(m,i,j2-1,k) = ch(m,i,k,j) + ch(m,i,k,jc)
1047 0 : cc(m,ic,j2-2,k) = ch(m,i,k,jc) - ch(m,i,k,j)
1048 0 : 570 END DO
1049 0 : 580 END DO
1050 0 : 590 END DO
1051 0 : 600 END DO
1052 : RETURN
1053 : END subroutine
1054 0 : SUBROUTINE vrftf1(m,n,c,ch,mdimc,wa,fac)
1055 :
1056 : ! VRFFTPK, VERSION 1, AUGUST 1985
1057 :
1058 : DIMENSION ch(mdimc,n),c(mdimc,n),wa(n),fac(15)
1059 :
1060 0 : nf = fac(2)
1061 0 : na = 1
1062 0 : l2 = n
1063 0 : iw = n
1064 0 : DO 110 k1 = 1,nf
1065 0 : kh = nf - k1
1066 0 : ip = fac(kh+3)
1067 0 : l1 = l2/ip
1068 0 : ido = n/l2
1069 0 : idl1 = ido*l1
1070 0 : iw = iw - (ip-1)*ido
1071 0 : na = 1 - na
1072 0 : IF (ip /= 4) GO TO 20
1073 0 : ix2 = iw + ido
1074 0 : ix3 = ix2 + ido
1075 0 : IF (na /= 0) GO TO 10
1076 0 : CALL vradf4(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3))
1077 0 : GO TO 100
1078 0 : 10 CALL vradf4(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3))
1079 0 : GO TO 100
1080 0 : 20 IF (ip /= 2) GO TO 40
1081 0 : IF (na /= 0) GO TO 30
1082 0 : CALL vradf2(m,ido,l1,c,ch,mdimc,wa(iw))
1083 0 : GO TO 100
1084 0 : 30 CALL vradf2(m,ido,l1,ch,c,mdimc,wa(iw))
1085 0 : GO TO 100
1086 0 : 40 IF (ip /= 3) GO TO 60
1087 0 : ix2 = iw + ido
1088 0 : IF (na /= 0) GO TO 50
1089 0 : CALL vradf3(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2))
1090 0 : GO TO 100
1091 0 : 50 CALL vradf3(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2))
1092 0 : GO TO 100
1093 0 : 60 IF (ip /= 5) GO TO 80
1094 0 : ix2 = iw + ido
1095 0 : ix3 = ix2 + ido
1096 0 : ix4 = ix3 + ido
1097 0 : IF (na /= 0) GO TO 70
1098 0 : CALL vradf5(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
1099 0 : GO TO 100
1100 0 : 70 CALL vradf5(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
1101 0 : GO TO 100
1102 0 : 80 IF (ido == 1) na = 1 - na
1103 0 : IF (na /= 0) GO TO 90
1104 0 : CALL vradfg(m,ido,ip,l1,idl1,c,c,c,ch,ch,mdimc,wa(iw))
1105 0 : na = 1
1106 0 : GO TO 100
1107 0 : 90 CALL vradfg(m,ido,ip,l1,idl1,ch,ch,ch,c,c,mdimc,wa(iw))
1108 0 : na = 0
1109 0 : 100 l2 = l1
1110 0 : 110 END DO
1111 0 : scale = sqrt(1./n)
1112 0 : IF (na == 1) GO TO 140
1113 0 : DO 130 j = 1,n
1114 0 : DO 120 i = 1,m
1115 0 : c(i,j) = scale*ch(i,j)
1116 0 : 120 END DO
1117 0 : 130 END DO
1118 0 : RETURN
1119 0 : 140 DO 160 j = 1,n
1120 0 : DO 150 i = 1,m
1121 0 : c(i,j) = scale*c(i,j)
1122 0 : 150 END DO
1123 0 : 160 END DO
1124 : RETURN
1125 : END subroutine
1126 :
1127 0 : SUBROUTINE vrfftb(m,n,r,rt,mdimr,wsave)
1128 : !***BEGIN PROLOGUE VRFFTB
1129 : !***DATE WRITTEN 850801 (YYMMDD)
1130 : !***REVISION DATE 900509 (YYMMDD)
1131 : !***CATEGORY NO. J1A1
1132 : !***KEYWORDS FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
1133 : ! FOURIER SYNTHESIS, BACKWARD TRANSFORM, MULTIPLE SEQUENCES
1134 : !***AUTHOR SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
1135 : !***PURPOSE Backward real periodic transform, M sequences.
1136 : !***DESCRIPTION
1137 :
1138 : ! Subroutine VRFFTB computes the synthesis (backward transform) of a
1139 : ! number of real periodic sequences from their Fourier coefficients.
1140 : ! Specifically, for each set of independent Fourier coefficients
1141 : ! F(K), the corresponding real periodic sequence is computed.
1142 :
1143 : ! The array WSAVE which is used by subroutine VRFFTB must be
1144 : ! initialized by calling subroutine VRFFTI(N,WSAVE).
1145 :
1146 : ! Input Parameters
1147 :
1148 : ! M the number of sets of coefficients.
1149 :
1150 : ! N the length of the sequences of coefficients to be
1151 : ! transformed. The method is most efficient when N is a
1152 : ! product of small primes, however n may be any positive
1153 : ! integer.
1154 :
1155 : ! R areal two-dimensional array of size MDIMX x N containing the
1156 : ! coefficients to be transformed. Each set of coefficients
1157 : ! F(K), K\0,1,..,N-1, is stored as a ROW of R. Specifically,
1158 : ! the I-th set of independent Fourier coefficients is stored
1159 :
1160 : ! R(I,1) = REAL( F(I,0) ),
1161 :
1162 : ! R(I,2*K) = REAL( F(I,K) )
1163 :
1164 : ! R(I,2*K+1) = IMAG( F(I,K) )
1165 :
1166 : ! for K = 1, 2, . . . , M-1,
1167 :
1168 : ! and, when N is even,
1169 :
1170 : ! R(I,N) = REAL( F(I,N/2) ).
1171 :
1172 : ! RT a real two-dimensional work array of size MDIMX x N.
1173 :
1174 : ! MDIMR the row (or first) dimension of the arrays R and RT exactly
1175 : ! as they appear in the calling program. This parameter is
1176 : ! used to specify the variable dimension of these arrays.
1177 :
1178 : ! WSAVE a real one-dimensional work array which must be dimensioned
1179 : ! at least N+15. The WSAVE array must be initialized by
1180 : ! calling subroutine VRFFTI. A different WSAVE array must be
1181 : ! used for each different value of N. This initialization does
1182 : ! not have to be repeated so long as N remains unchanged. The
1183 : ! same WSAVE array may be used by VRFFTB and VRFFTB.
1184 :
1185 : ! Output Parameters
1186 :
1187 : ! R contains M real periodic sequences corresponding to the given
1188 : ! coefficients. Specifically, the I-th row of R contains the
1189 : ! real periodic sequence corresponding to the I-th set of
1190 : ! independent Fourier coefficients F(I,K) stored as
1191 :
1192 : ! R(I,J) = X(I,J-1) , J = 1, 2, . . . , N, where
1193 :
1194 : ! X(I,J) = SQRT(1/N)* F(I,0) + (-1)**J*F(I,N/2)
1195 : ! + 2*SUM(K=1,M)[ REAL(F(I,2K))*COS(2K*J*PI/N)
1196 : ! - IMAG(F(I,2K+1))*SIN(2K*J*PI/N) ] ,
1197 :
1198 : ! when N is even, and
1199 :
1200 : ! X(I,J) = SQRT(1/N)* F(I,0) +
1201 : ! 2*SUM(K=1,M)[ REAL(F(I,2K))*COS(2K*J*PI/N)
1202 : ! - IMAG(F(I,2K+1))*SIN(2K*J*PI/N) ] ,
1203 :
1204 : ! when N is odd.
1205 :
1206 : ! WSAVE contains results which must not be destroyed between calls
1207 : ! to VRFFTF or VRFFTB.
1208 :
1209 : ! -----------------------------------------------------------------
1210 :
1211 : ! NOTE - A call of VRFFTF followed immediately by a call of
1212 : ! of VRFFTB will return the original sequences R. Thus,
1213 : ! VRFFTB is the correctly normalized inverse of VRFFTF.
1214 :
1215 : ! -----------------------------------------------------------------
1216 :
1217 : ! VRFFTB is a straightforward extension of the subprogram RFFTB to
1218 : ! handle M simultaneous sequences. RFFTB was originally developed
1219 : ! by P. N. Swarztrauber of NCAR.
1220 :
1221 : ! * * * * * * * * * * * * * * * * * * * * *
1222 : ! * *
1223 : ! * PROGRAM SPECIFICATIONS *
1224 : ! * *
1225 : ! * * * * * * * * * * * * * * * * * * * * *
1226 :
1227 : ! Dimension of R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
1228 : ! Arguments
1229 :
1230 : ! Latest AUGUST 1, 1985
1231 : ! Revision
1232 :
1233 : ! Subprograms VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
1234 : ! required VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
1235 : ! VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
1236 :
1237 : ! Special NONE
1238 : ! Conditions
1239 :
1240 : ! Common NONE
1241 : ! blocks
1242 :
1243 : ! I/O NONE
1244 :
1245 : ! Precision SINGLE
1246 :
1247 : ! Specialist ROLAND SWEET
1248 :
1249 : ! Language FORTRAN
1250 :
1251 : ! History WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
1252 : ! NATIONAL BUREAU OF STANDARDS (BOULDER).
1253 :
1254 : ! Algorithm A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
1255 : ! OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
1256 :
1257 : ! Portability AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
1258 : ! THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
1259 : ! THE FUNCTION PIMACH.
1260 :
1261 : ! Required COS,SIN
1262 : ! resident
1263 : ! Routines
1264 :
1265 : !***REFERENCES P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
1266 : ! Computations, (G. Rodrigue, ed.), Academic Press, 1982,
1267 : ! pp. 51-83.
1268 : !***ROUTINES CALLED VRFTB1
1269 : !***END PROLOGUE VRFFTB
1270 :
1271 : ! VRFFTPK, VERSION 1, AUGUST 1985
1272 :
1273 : IMPLICIT NONE
1274 : INTEGER, INTENT(IN) :: n, m, mdimr
1275 : REAL :: r(mdimr,n),rt(mdimr,n),wsave(n+15)
1276 :
1277 0 : IF (n == 1) RETURN
1278 0 : CALL vrftb1(m,n,r,rt,mdimr,wsave(1),wsave(n+1))
1279 0 : RETURN
1280 :
1281 : CONTAINS
1282 :
1283 0 : SUBROUTINE vradb2(mp,ido,l1,cc,ch,mdimc,wa1)
1284 :
1285 : ! VRFFTPK, VERSION 1, AUGUST 1985
1286 :
1287 : IMPLICIT NONE
1288 : INTEGER, INTENT(IN) :: mp, ido, l1, mdimc
1289 : REAL, INTENT(IN) :: cc(mdimc,ido,2,l1), wa1(ido)
1290 : REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,2)
1291 :
1292 : INTEGER :: i, m, k, ic, idp2
1293 :
1294 0 : DO 20 k = 1,l1
1295 0 : DO 10 m = 1,mp
1296 0 : ch(m,1,k,1) = cc(m,1,1,k) + cc(m,ido,2,k)
1297 0 : ch(m,1,k,2) = cc(m,1,1,k) - cc(m,ido,2,k)
1298 0 : 10 END DO
1299 0 : 20 END DO
1300 0 : IF ( ido == 2 ) THEN
1301 : GOTO 70
1302 0 : ELSEIF ( ido < 2 ) THEN
1303 : GOTO 100
1304 : ENDIF
1305 0 : idp2 = ido + 2
1306 0 : DO 60 k = 1,l1
1307 0 : DO 50 i = 3,ido,2
1308 0 : ic = idp2 - i
1309 0 : DO 40 m = 1,mp
1310 0 : ch(m,i-1,k,1) = cc(m,i-1,1,k) + cc(m,ic-1,2,k)
1311 0 : ch(m,i,k,1) = cc(m,i,1,k) - cc(m,ic,2,k)
1312 : ch(m,i-1,k,2) = wa1(i-2)* (cc(m,i-1,1,k)- &
1313 : cc(m,ic-1,2,k)) - wa1(i-1)* &
1314 0 : (cc(m,i,1,k)+cc(m,ic,2,k))
1315 : ch(m,i,k,2) = wa1(i-2)* (cc(m,i,1,k)+cc(m,ic,2,k)) + &
1316 0 : wa1(i-1)* (cc(m,i-1,1,k)-cc(m,ic-1,2,k))
1317 0 : 40 END DO
1318 0 : 50 END DO
1319 0 : 60 END DO
1320 0 : IF (mod(ido,2) == 1) RETURN
1321 0 : 70 DO 90 k = 1,l1
1322 0 : DO 80 m = 1,mp
1323 0 : ch(m,ido,k,1) = cc(m,ido,1,k) + cc(m,ido,1,k)
1324 0 : ch(m,ido,k,2) = - (cc(m,1,2,k)+cc(m,1,2,k))
1325 0 : 80 END DO
1326 0 : 90 END DO
1327 : 100 RETURN
1328 : END SUBROUTINE vradb2
1329 :
1330 0 : SUBROUTINE vradb3(mp,ido,l1,cc,ch,mdimc,wa1,wa2)
1331 :
1332 : ! VRFFTPK, VERSION 1, AUGUST 1985
1333 :
1334 : USE m_constants, ONLY : pimach
1335 : IMPLICIT NONE
1336 : INTEGER, INTENT(IN) :: mp, ido, l1, mdimc
1337 : REAL, INTENT(IN) :: cc(mdimc,ido,3,l1), wa1(ido), wa2(ido)
1338 : REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,3)
1339 :
1340 : INTEGER :: i, k, idp2, ic, m
1341 : REAL :: arg, taur, taui
1342 :
1343 0 : arg = 2.*pimach()/3.
1344 0 : taur = cos(arg)
1345 0 : taui = sin(arg)
1346 0 : DO 20 k = 1,l1
1347 0 : DO 10 m = 1,mp
1348 0 : ch(m,1,k,1) = cc(m,1,1,k) + 2.*cc(m,ido,2,k)
1349 : ch(m,1,k,2) = cc(m,1,1,k) + (2.*taur)*cc(m,ido,2,k) - &
1350 0 : (2.*taui)*cc(m,1,3,k)
1351 : ch(m,1,k,3) = cc(m,1,1,k) + (2.*taur)*cc(m,ido,2,k) + &
1352 0 : & 2.*taui*cc(m,1,3,k)
1353 0 : 10 END DO
1354 0 : 20 END DO
1355 0 : IF (ido == 1) RETURN
1356 0 : idp2 = ido + 2
1357 0 : DO 50 k = 1,l1
1358 0 : DO 40 i = 3,ido,2
1359 0 : ic = idp2 - i
1360 0 : DO 30 m = 1,mp
1361 : ch(m,i-1,k,1) = cc(m,i-1,1,k) + &
1362 0 : (cc(m,i-1,3,k)+cc(m,ic-1,2,k))
1363 0 : ch(m,i,k,1) = cc(m,i,1,k) + (cc(m,i,3,k)-cc(m,ic,2,k))
1364 : ch(m,i-1,k,2) = wa1(i-2)* ((cc(m,i-1,1,k)+taur* (cc(m, &
1365 : i-1,3,k)+cc(m,ic-1,2,k)))- &
1366 : (taui* (cc(m,i,3,k)+cc(m,ic,2,k)))) - &
1367 : wa1(i-1)* ((cc(m,i,1,k)+taur* (cc(m,i,3, &
1368 : k)-cc(m,ic,2,k)))+ (taui* (cc(m,i-1,3, &
1369 0 : k)-cc(m,ic-1,2,k))))
1370 : ch(m,i,k,2) = wa1(i-2)* ((cc(m,i,1,k)+taur* (cc(m,i,3, &
1371 : k)-cc(m,ic,2,k)))+ (taui* (cc(m,i-1,3, &
1372 : k)-cc(m,ic-1,2,k)))) + &
1373 : wa1(i-1)* ((cc(m,i-1,1,k)+taur* (cc(m,i-1, &
1374 : & 3,k)+cc(m,ic-1,2,k)))- &
1375 0 : (taui* (cc(m,i,3,k)+cc(m,ic,2,k))))
1376 : ch(m,i-1,k,3) = wa2(i-2)* ((cc(m,i-1,1,k)+taur* (cc(m, &
1377 : i-1,3,k)+cc(m,ic-1,2,k)))+ &
1378 : (taui* (cc(m,i,3,k)+cc(m,ic,2,k)))) - &
1379 : wa2(i-1)* ((cc(m,i,1,k)+taur* (cc(m,i,3, &
1380 : k)-cc(m,ic,2,k)))- (taui* (cc(m,i-1,3, &
1381 0 : k)-cc(m,ic-1,2,k))))
1382 : ch(m,i,k,3) = wa2(i-2)* ((cc(m,i,1,k)+taur* (cc(m,i,3, &
1383 : k)-cc(m,ic,2,k)))- (taui* (cc(m,i-1,3, &
1384 : k)-cc(m,ic-1,2,k)))) + &
1385 : wa2(i-1)* ((cc(m,i-1,1,k)+taur* (cc(m,i-1, &
1386 : & 3,k)+cc(m,ic-1,2,k)))+ &
1387 0 : (taui* (cc(m,i,3,k)+cc(m,ic,2,k))))
1388 0 : 30 END DO
1389 0 : 40 END DO
1390 0 : 50 END DO
1391 : RETURN
1392 : END SUBROUTINE vradb3
1393 :
1394 0 : SUBROUTINE vradb4(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3)
1395 :
1396 : ! VRFFTPK, VERSION 1, AUGUST 1985
1397 :
1398 : IMPLICIT NONE
1399 : INTEGER, INTENT(IN) :: mp, ido, l1, mdimc
1400 : REAL, INTENT(IN) :: cc(mdimc,ido,4,l1), wa1(ido), wa2(ido), &
1401 : wa3(ido)
1402 : REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,4)
1403 :
1404 : INTEGER :: k, m, i, idp2, ic
1405 : REAL :: sqrt2
1406 :
1407 0 : sqrt2 = sqrt(2.)
1408 0 : DO 20 k = 1,l1
1409 0 : DO 10 m = 1,mp
1410 : ch(m,1,k,3) = (cc(m,1,1,k)+cc(m,ido,4,k)) - &
1411 0 : (cc(m,ido,2,k)+cc(m,ido,2,k))
1412 : ch(m,1,k,1) = (cc(m,1,1,k)+cc(m,ido,4,k)) + &
1413 0 : (cc(m,ido,2,k)+cc(m,ido,2,k))
1414 : ch(m,1,k,4) = (cc(m,1,1,k)-cc(m,ido,4,k)) + &
1415 0 : (cc(m,1,3,k)+cc(m,1,3,k))
1416 : ch(m,1,k,2) = (cc(m,1,1,k)-cc(m,ido,4,k)) - &
1417 0 : (cc(m,1,3,k)+cc(m,1,3,k))
1418 0 : 10 END DO
1419 0 : 20 END DO
1420 0 : IF ( ido == 2 ) THEN
1421 : GOTO 70
1422 0 : ELSEIF ( ido < 2 ) THEN
1423 : GOTO 100
1424 : ENDIF
1425 0 : idp2 = ido + 2
1426 0 : DO 60 k = 1,l1
1427 0 : DO 50 i = 3,ido,2
1428 0 : ic = idp2 - i
1429 0 : DO 40 m = 1,mp
1430 : ch(m,i-1,k,1) = (cc(m,i-1,1,k)+cc(m,ic-1,4,k)) + &
1431 0 : (cc(m,i-1,3,k)+cc(m,ic-1,2,k))
1432 : ch(m,i,k,1) = (cc(m,i,1,k)-cc(m,ic,4,k)) + &
1433 0 : (cc(m,i,3,k)-cc(m,ic,2,k))
1434 : ch(m,i-1,k,2) = wa1(i-2)* ((cc(m,i-1,1,k)-cc(m,ic-1,4, &
1435 : k))- (cc(m,i,3,k)+cc(m,ic,2,k))) - &
1436 : wa1(i-1)* ((cc(m,i,1,k)+cc(m,ic,4,k))+ &
1437 0 : (cc(m,i-1,3,k)-cc(m,ic-1,2,k)))
1438 : ch(m,i,k,2) = wa1(i-2)* ((cc(m,i,1,k)+cc(m,ic,4,k))+ &
1439 : (cc(m,i-1,3,k)-cc(m,ic-1,2,k))) + &
1440 : wa1(i-1)* ((cc(m,i-1,1,k)-cc(m,ic-1,4,k))- &
1441 0 : (cc(m,i,3,k)+cc(m,ic,2,k)))
1442 : ch(m,i-1,k,3) = wa2(i-2)* ((cc(m,i-1,1,k)+cc(m,ic-1,4, &
1443 : k))- (cc(m,i-1,3,k)+cc(m,ic-1,2,k))) - &
1444 : wa2(i-1)* ((cc(m,i,1,k)-cc(m,ic,4,k))- &
1445 0 : (cc(m,i,3,k)-cc(m,ic,2,k)))
1446 : ch(m,i,k,3) = wa2(i-2)* ((cc(m,i,1,k)-cc(m,ic,4,k))- &
1447 : (cc(m,i,3,k)-cc(m,ic,2,k))) + &
1448 : wa2(i-1)* ((cc(m,i-1,1,k)+cc(m,ic-1,4,k))- &
1449 0 : (cc(m,i-1,3,k)+cc(m,ic-1,2,k)))
1450 : ch(m,i-1,k,4) = wa3(i-2)* ((cc(m,i-1,1,k)-cc(m,ic-1,4, &
1451 : k))+ (cc(m,i,3,k)+cc(m,ic,2,k))) - &
1452 : wa3(i-1)* ((cc(m,i,1,k)+cc(m,ic,4,k))- &
1453 0 : (cc(m,i-1,3,k)-cc(m,ic-1,2,k)))
1454 : ch(m,i,k,4) = wa3(i-2)* ((cc(m,i,1,k)+cc(m,ic,4,k))- &
1455 : (cc(m,i-1,3,k)-cc(m,ic-1,2,k))) + &
1456 : wa3(i-1)* ((cc(m,i-1,1,k)-cc(m,ic-1,4,k))+ &
1457 0 : (cc(m,i,3,k)+cc(m,ic,2,k)))
1458 0 : 40 END DO
1459 0 : 50 END DO
1460 0 : 60 END DO
1461 0 : IF (mod(ido,2) == 1) RETURN
1462 : 70 CONTINUE
1463 0 : DO 90 k = 1,l1
1464 0 : DO 80 m = 1,mp
1465 : ch(m,ido,k,1) = (cc(m,ido,1,k)+cc(m,ido,3,k)) + &
1466 0 : (cc(m,ido,1,k)+cc(m,ido,3,k))
1467 : ch(m,ido,k,2) = sqrt2* ((cc(m,ido,1,k)-cc(m,ido,3,k))- &
1468 0 : (cc(m,1,2,k)+cc(m,1,4,k)))
1469 : ch(m,ido,k,3) = (cc(m,1,4,k)-cc(m,1,2,k)) + &
1470 0 : (cc(m,1,4,k)-cc(m,1,2,k))
1471 : ch(m,ido,k,4) = -sqrt2* ((cc(m,ido,1,k)-cc(m,ido,3,k))+ &
1472 0 : (cc(m,1,2,k)+cc(m,1,4,k)))
1473 0 : 80 END DO
1474 0 : 90 END DO
1475 : 100 RETURN
1476 : END SUBROUTINE vradb4
1477 :
1478 0 : SUBROUTINE vradb5(mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3,wa4)
1479 :
1480 : ! VRFFTPK, VERSION 1, AUGUST 1985
1481 :
1482 : USE m_constants, ONLY : pimach
1483 : IMPLICIT NONE
1484 : INTEGER, INTENT(IN) :: mp, ido, l1, mdimc
1485 : REAL, INTENT(IN) :: cc(mdimc,ido,5,l1), wa1(ido), wa2(ido), &
1486 : wa3(ido),wa4(ido)
1487 : REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,5)
1488 :
1489 : INTEGER :: i, k, m, ic, idp2
1490 : REAL :: arg, tr11, ti11, tr12, ti12
1491 :
1492 0 : arg = 2.*pimach()/5.
1493 0 : tr11 = cos(arg)
1494 0 : ti11 = sin(arg)
1495 0 : tr12 = cos(2.*arg)
1496 0 : ti12 = sin(2.*arg)
1497 0 : DO 20 k = 1,l1
1498 0 : DO 10 m = 1,mp
1499 : ch(m,1,k,1) = cc(m,1,1,k) + 2.*cc(m,ido,2,k) + &
1500 0 : & 2.*cc(m,ido,4,k)
1501 : ch(m,1,k,2) = (cc(m,1,1,k)+tr11*2.*cc(m,ido,2,k)+ &
1502 : tr12*2.*cc(m,ido,4,k)) - &
1503 0 : (ti11*2.*cc(m,1,3,k)+ti12*2.*cc(m,1,5,k))
1504 : ch(m,1,k,3) = (cc(m,1,1,k)+tr12*2.*cc(m,ido,2,k)+ &
1505 : tr11*2.*cc(m,ido,4,k)) - &
1506 0 : (ti12*2.*cc(m,1,3,k)-ti11*2.*cc(m,1,5,k))
1507 : ch(m,1,k,4) = (cc(m,1,1,k)+tr12*2.*cc(m,ido,2,k)+ &
1508 : tr11*2.*cc(m,ido,4,k)) + &
1509 0 : (ti12*2.*cc(m,1,3,k)-ti11*2.*cc(m,1,5,k))
1510 : ch(m,1,k,5) = (cc(m,1,1,k)+tr11*2.*cc(m,ido,2,k)+ &
1511 : tr12*2.*cc(m,ido,4,k)) + &
1512 0 : (ti11*2.*cc(m,1,3,k)+ti12*2.*cc(m,1,5,k))
1513 0 : 10 END DO
1514 0 : 20 END DO
1515 0 : IF (ido == 1) RETURN
1516 0 : idp2 = ido + 2
1517 0 : DO 50 k = 1,l1
1518 0 : DO 40 i = 3,ido,2
1519 0 : ic = idp2 - i
1520 0 : DO 30 m = 1,mp
1521 : ch(m,i-1,k,1) = cc(m,i-1,1,k) + &
1522 : (cc(m,i-1,3,k)+cc(m,ic-1,2,k)) + &
1523 0 : (cc(m,i-1,5,k)+cc(m,ic-1,4,k))
1524 : ch(m,i,k,1) = cc(m,i,1,k) + (cc(m,i,3,k)-cc(m,ic,2,k)) + &
1525 0 : (cc(m,i,5,k)-cc(m,ic,4,k))
1526 : ch(m,i-1,k,2) = wa1(i-2)* ((cc(m,i-1,1,k)+tr11* (cc(m, &
1527 : i-1,3,k)+cc(m,ic-1,2,k))+tr12* (cc(m,i-1, &
1528 : & 5,k)+cc(m,ic-1,4,k)))- &
1529 : (ti11* (cc(m,i,3,k)+cc(m,ic,2, &
1530 : k))+ti12* (cc(m,i,5,k)+cc(m,ic,4,k)))) - &
1531 : wa1(i-1)* ((cc(m,i,1,k)+tr11* (cc(m,i,3, &
1532 : k)-cc(m,ic,2,k))+tr12* (cc(m,i,5,k)-cc(m, &
1533 : ic,4,k)))+ (ti11* (cc(m,i-1,3,k)-cc(m, &
1534 : ic-1,2,k))+ti12* (cc(m,i-1,5,k)-cc(m, &
1535 0 : ic-1,4,k))))
1536 : ch(m,i,k,2) = wa1(i-2)* ((cc(m,i,1,k)+tr11* (cc(m,i,3, &
1537 : k)-cc(m,ic,2,k))+tr12* (cc(m,i,5,k)-cc(m, &
1538 : ic,4,k)))+ (ti11* (cc(m,i-1,3,k)-cc(m,ic-1, &
1539 : & 2,k))+ti12* (cc(m,i-1,5,k)-cc(m,ic-1,4, &
1540 : k)))) + wa1(i-1)* ((cc(m,i-1,1, &
1541 : k)+tr11* (cc(m,i-1,3,k)+cc(m,ic-1,2, &
1542 : k))+tr12* (cc(m,i-1,5,k)+cc(m,ic-1,4,k)))- &
1543 : (ti11* (cc(m,i,3,k)+cc(m,ic,2, &
1544 0 : k))+ti12* (cc(m,i,5,k)+cc(m,ic,4,k))))
1545 : ch(m,i-1,k,3) = wa2(i-2)* ((cc(m,i-1,1,k)+tr12* (cc(m, &
1546 : i-1,3,k)+cc(m,ic-1,2,k))+tr11* (cc(m,i-1, &
1547 : & 5,k)+cc(m,ic-1,4,k)))- &
1548 : (ti12* (cc(m,i,3,k)+cc(m,ic,2, &
1549 : k))-ti11* (cc(m,i,5,k)+cc(m,ic,4,k)))) - &
1550 : wa2(i-1)* ((cc(m,i,1,k)+tr12* (cc(m,i,3, &
1551 : k)-cc(m,ic,2,k))+tr11* (cc(m,i,5,k)-cc(m, &
1552 : ic,4,k)))+ (ti12* (cc(m,i-1,3,k)-cc(m, &
1553 : ic-1,2,k))-ti11* (cc(m,i-1,5,k)-cc(m, &
1554 0 : ic-1,4,k))))
1555 : ch(m,i,k,3) = wa2(i-2)* ((cc(m,i,1,k)+tr12* (cc(m,i,3, &
1556 : k)-cc(m,ic,2,k))+tr11* (cc(m,i,5,k)-cc(m, &
1557 : ic,4,k)))+ (ti12* (cc(m,i-1,3,k)-cc(m,ic-1, &
1558 : & 2,k))-ti11* (cc(m,i-1,5,k)-cc(m,ic-1,4, &
1559 : k)))) + wa2(i-1)* ((cc(m,i-1,1, &
1560 : k)+tr12* (cc(m,i-1,3,k)+cc(m,ic-1,2, &
1561 : k))+tr11* (cc(m,i-1,5,k)+cc(m,ic-1,4,k)))- &
1562 : (ti12* (cc(m,i,3,k)+cc(m,ic,2, &
1563 0 : k))-ti11* (cc(m,i,5,k)+cc(m,ic,4,k))))
1564 : ch(m,i-1,k,4) = wa3(i-2)* ((cc(m,i-1,1,k)+tr12* (cc(m, &
1565 : i-1,3,k)+cc(m,ic-1,2,k))+tr11* (cc(m,i-1, &
1566 : & 5,k)+cc(m,ic-1,4,k)))+ &
1567 : (ti12* (cc(m,i,3,k)+cc(m,ic,2, &
1568 : k))-ti11* (cc(m,i,5,k)+cc(m,ic,4,k)))) - &
1569 : wa3(i-1)* ((cc(m,i,1,k)+tr12* (cc(m,i,3, &
1570 : k)-cc(m,ic,2,k))+tr11* (cc(m,i,5,k)-cc(m, &
1571 : ic,4,k)))- (ti12* (cc(m,i-1,3,k)-cc(m, &
1572 : ic-1,2,k))-ti11* (cc(m,i-1,5,k)-cc(m, &
1573 0 : ic-1,4,k))))
1574 : ch(m,i,k,4) = wa3(i-2)* ((cc(m,i,1,k)+tr12* (cc(m,i,3, &
1575 : k)-cc(m,ic,2,k))+tr11* (cc(m,i,5,k)-cc(m, &
1576 : ic,4,k)))- (ti12* (cc(m,i-1,3,k)-cc(m,ic-1, &
1577 : & 2,k))-ti11* (cc(m,i-1,5,k)-cc(m,ic-1,4, &
1578 : k)))) + wa3(i-1)* ((cc(m,i-1,1, &
1579 : k)+tr12* (cc(m,i-1,3,k)+cc(m,ic-1,2, &
1580 : k))+tr11* (cc(m,i-1,5,k)+cc(m,ic-1,4,k)))+ &
1581 : (ti12* (cc(m,i,3,k)+cc(m,ic,2, &
1582 0 : k))-ti11* (cc(m,i,5,k)+cc(m,ic,4,k))))
1583 : ch(m,i-1,k,5) = wa4(i-2)* ((cc(m,i-1,1,k)+tr11* (cc(m, &
1584 : i-1,3,k)+cc(m,ic-1,2,k))+tr12* (cc(m,i-1, &
1585 : & 5,k)+cc(m,ic-1,4,k)))+ &
1586 : (ti11* (cc(m,i,3,k)+cc(m,ic,2, &
1587 : k))+ti12* (cc(m,i,5,k)+cc(m,ic,4,k)))) - &
1588 : wa4(i-1)* ((cc(m,i,1,k)+tr11* (cc(m,i,3, &
1589 : k)-cc(m,ic,2,k))+tr12* (cc(m,i,5,k)-cc(m, &
1590 : ic,4,k)))- (ti11* (cc(m,i-1,3,k)-cc(m, &
1591 : ic-1,2,k))+ti12* (cc(m,i-1,5,k)-cc(m, &
1592 0 : ic-1,4,k))))
1593 : ch(m,i,k,5) = wa4(i-2)* ((cc(m,i,1,k)+tr11* (cc(m,i,3, &
1594 : k)-cc(m,ic,2,k))+tr12* (cc(m,i,5,k)-cc(m, &
1595 : ic,4,k)))- (ti11* (cc(m,i-1,3,k)-cc(m,ic-1, &
1596 : & 2,k))+ti12* (cc(m,i-1,5,k)-cc(m,ic-1,4, &
1597 : k)))) + wa4(i-1)* ((cc(m,i-1,1, &
1598 : k)+tr11* (cc(m,i-1,3,k)+cc(m,ic-1,2, &
1599 : k))+tr12* (cc(m,i-1,5,k)+cc(m,ic-1,4,k)))+ &
1600 : (ti11* (cc(m,i,3,k)+cc(m,ic,2, &
1601 0 : k))+ti12* (cc(m,i,5,k)+cc(m,ic,4,k))))
1602 0 : 30 END DO
1603 0 : 40 END DO
1604 0 : 50 END DO
1605 : RETURN
1606 : END SUBROUTINE vradb5
1607 :
1608 0 : SUBROUTINE vradbg(mp,ido,ip,l1,idl1,cc,c1,c2,ch,ch2, &
1609 0 : mdimc,wa)
1610 :
1611 : ! VRFFTPK, VERSION 1, AUGUST 1985
1612 :
1613 : USE m_constants, ONLY : pimach
1614 : IMPLICIT NONE
1615 : INTEGER, INTENT(IN) :: mp, ido, ip, l1, idl1, mdimc
1616 : REAL, INTENT(IN) :: cc(mdimc,ido,ip,l1), wa(ido)
1617 : REAL, INTENT(INOUT) :: ch(mdimc,ido,l1,ip), ch2(mdimc,idl1,ip), &
1618 : c1(mdimc,ido,l1,ip), c2(mdimc,idl1,ip)
1619 :
1620 : INTEGER :: i, j, j2, k, l, m, ik, is, idp2, nbd, ipp2, ipph, &
1621 : ic, jc, lc, idij
1622 : REAL :: ai1, ar1, ar2, ai2, ar1h, ar2h, tpi, arg, dcp, dsp, &
1623 : dc2, ds2
1624 :
1625 0 : tpi = 2.*pimach()
1626 0 : arg = tpi/real(ip)
1627 0 : dcp = cos(arg)
1628 0 : dsp = sin(arg)
1629 0 : idp2 = ido + 2
1630 0 : nbd = (ido-1)/2
1631 0 : ipp2 = ip + 2
1632 0 : ipph = (ip+1)/2
1633 0 : IF (ido < l1) GO TO 40
1634 0 : DO 30 k = 1,l1
1635 0 : DO 20 i = 1,ido
1636 0 : DO 10 m = 1,mp
1637 0 : ch(m,i,k,1) = cc(m,i,1,k)
1638 0 : 10 END DO
1639 0 : 20 END DO
1640 0 : 30 END DO
1641 0 : GO TO 80
1642 0 : 40 DO 70 i = 1,ido
1643 0 : DO 60 k = 1,l1
1644 0 : DO 50 m = 1,mp
1645 0 : ch(m,i,k,1) = cc(m,i,1,k)
1646 0 : 50 END DO
1647 0 : 60 END DO
1648 0 : 70 END DO
1649 0 : 80 DO 110 j = 2,ipph
1650 0 : jc = ipp2 - j
1651 0 : j2 = j + j
1652 0 : DO 100 k = 1,l1
1653 0 : DO 90 m = 1,mp
1654 0 : ch(m,1,k,j) = cc(m,ido,j2-2,k) + cc(m,ido,j2-2,k)
1655 0 : ch(m,1,k,jc) = cc(m,1,j2-1,k) + cc(m,1,j2-1,k)
1656 0 : 90 END DO
1657 0 : 100 END DO
1658 0 : 110 END DO
1659 0 : IF (ido == 1) GO TO 210
1660 0 : IF (nbd < l1) GO TO 160
1661 0 : DO 150 j = 2,ipph
1662 0 : jc = ipp2 - j
1663 0 : DO 140 k = 1,l1
1664 0 : DO 130 i = 3,ido,2
1665 0 : ic = idp2 - i
1666 0 : DO 120 m = 1,mp
1667 0 : ch(m,i-1,k,j) = cc(m,i-1,2*j-1,k) + cc(m,ic-1,2*j-2,k)
1668 : ch(m,i-1,k,jc) = cc(m,i-1,2*j-1,k) - &
1669 0 : cc(m,ic-1,2*j-2,k)
1670 0 : ch(m,i,k,j) = cc(m,i,2*j-1,k) - cc(m,ic,2*j-2,k)
1671 0 : ch(m,i,k,jc) = cc(m,i,2*j-1,k) + cc(m,ic,2*j-2,k)
1672 0 : 120 END DO
1673 0 : 130 END DO
1674 0 : 140 END DO
1675 0 : 150 END DO
1676 0 : GO TO 210
1677 0 : 160 DO 200 j = 2,ipph
1678 0 : jc = ipp2 - j
1679 0 : DO 190 i = 3,ido,2
1680 0 : ic = idp2 - i
1681 0 : DO 180 k = 1,l1
1682 0 : DO 170 m = 1,mp
1683 0 : ch(m,i-1,k,j) = cc(m,i-1,2*j-1,k) + cc(m,ic-1,2*j-2,k)
1684 : ch(m,i-1,k,jc) = cc(m,i-1,2*j-1,k) - &
1685 0 : cc(m,ic-1,2*j-2,k)
1686 0 : ch(m,i,k,j) = cc(m,i,2*j-1,k) - cc(m,ic,2*j-2,k)
1687 0 : ch(m,i,k,jc) = cc(m,i,2*j-1,k) + cc(m,ic,2*j-2,k)
1688 0 : 170 END DO
1689 0 : 180 END DO
1690 0 : 190 END DO
1691 0 : 200 END DO
1692 : 210 ar1 = 1.
1693 : ai1 = 0.
1694 0 : DO 270 l = 2,ipph
1695 0 : lc = ipp2 - l
1696 0 : ar1h = dcp*ar1 - dsp*ai1
1697 0 : ai1 = dcp*ai1 + dsp*ar1
1698 0 : ar1 = ar1h
1699 0 : DO 230 ik = 1,idl1
1700 0 : DO 220 m = 1,mp
1701 0 : c2(m,ik,l) = ch2(m,ik,1) + ar1*ch2(m,ik,2)
1702 0 : c2(m,ik,lc) = ai1*ch2(m,ik,ip)
1703 0 : 220 END DO
1704 0 : 230 END DO
1705 : dc2 = ar1
1706 : ds2 = ai1
1707 : ar2 = ar1
1708 : ai2 = ai1
1709 0 : DO 260 j = 3,ipph
1710 0 : jc = ipp2 - j
1711 0 : ar2h = dc2*ar2 - ds2*ai2
1712 0 : ai2 = dc2*ai2 + ds2*ar2
1713 0 : ar2 = ar2h
1714 0 : DO 250 ik = 1,idl1
1715 0 : DO 240 m = 1,mp
1716 0 : c2(m,ik,l) = c2(m,ik,l) + ar2*ch2(m,ik,j)
1717 0 : c2(m,ik,lc) = c2(m,ik,lc) + ai2*ch2(m,ik,jc)
1718 0 : 240 END DO
1719 0 : 250 END DO
1720 0 : 260 END DO
1721 0 : 270 END DO
1722 0 : DO 300 j = 2,ipph
1723 0 : DO 290 ik = 1,idl1
1724 0 : DO 280 m = 1,mp
1725 0 : ch2(m,ik,1) = ch2(m,ik,1) + ch2(m,ik,j)
1726 0 : 280 END DO
1727 0 : 290 END DO
1728 0 : 300 END DO
1729 0 : DO 330 j = 2,ipph
1730 0 : jc = ipp2 - j
1731 0 : DO 320 k = 1,l1
1732 0 : DO 310 m = 1,mp
1733 0 : ch(m,1,k,j) = c1(m,1,k,j) - c1(m,1,k,jc)
1734 0 : ch(m,1,k,jc) = c1(m,1,k,j) + c1(m,1,k,jc)
1735 0 : 310 END DO
1736 0 : 320 END DO
1737 0 : 330 END DO
1738 0 : IF (ido == 1) GO TO 430
1739 0 : IF (nbd < l1) GO TO 380
1740 0 : DO 370 j = 2,ipph
1741 0 : jc = ipp2 - j
1742 0 : DO 360 k = 1,l1
1743 0 : DO 350 i = 3,ido,2
1744 0 : DO 340 m = 1,mp
1745 0 : ch(m,i-1,k,j) = c1(m,i-1,k,j) - c1(m,i,k,jc)
1746 0 : ch(m,i-1,k,jc) = c1(m,i-1,k,j) + c1(m,i,k,jc)
1747 0 : ch(m,i,k,j) = c1(m,i,k,j) + c1(m,i-1,k,jc)
1748 0 : ch(m,i,k,jc) = c1(m,i,k,j) - c1(m,i-1,k,jc)
1749 0 : 340 END DO
1750 0 : 350 END DO
1751 0 : 360 END DO
1752 0 : 370 END DO
1753 0 : GO TO 430
1754 0 : 380 DO 420 j = 2,ipph
1755 0 : jc = ipp2 - j
1756 0 : DO 410 i = 3,ido,2
1757 0 : DO 400 k = 1,l1
1758 0 : DO 390 m = 1,mp
1759 0 : ch(m,i-1,k,j) = c1(m,i-1,k,j) - c1(m,i,k,jc)
1760 0 : ch(m,i-1,k,jc) = c1(m,i-1,k,j) + c1(m,i,k,jc)
1761 0 : ch(m,i,k,j) = c1(m,i,k,j) + c1(m,i-1,k,jc)
1762 0 : ch(m,i,k,jc) = c1(m,i,k,j) - c1(m,i-1,k,jc)
1763 0 : 390 END DO
1764 0 : 400 END DO
1765 0 : 410 END DO
1766 0 : 420 END DO
1767 : 430 CONTINUE
1768 0 : IF (ido == 1) RETURN
1769 0 : DO 450 ik = 1,idl1
1770 0 : DO 440 m = 1,mp
1771 0 : c2(m,ik,1) = ch2(m,ik,1)
1772 0 : 440 END DO
1773 0 : 450 END DO
1774 0 : DO 480 j = 2,ip
1775 0 : DO 470 k = 1,l1
1776 0 : DO 460 m = 1,mp
1777 0 : c1(m,1,k,j) = ch(m,1,k,j)
1778 0 : 460 END DO
1779 0 : 470 END DO
1780 0 : 480 END DO
1781 0 : IF (nbd > l1) GO TO 530
1782 0 : is = -ido
1783 0 : DO 520 j = 2,ip
1784 0 : is = is + ido
1785 0 : idij = is
1786 0 : DO 510 i = 3,ido,2
1787 0 : idij = idij + 2
1788 0 : DO 500 k = 1,l1
1789 0 : DO 490 m = 1,mp
1790 : c1(m,i-1,k,j) = wa(idij-1)*ch(m,i-1,k,j) - &
1791 0 : wa(idij)*ch(m,i,k,j)
1792 : c1(m,i,k,j) = wa(idij-1)*ch(m,i,k,j) + &
1793 0 : wa(idij)*ch(m,i-1,k,j)
1794 0 : 490 END DO
1795 0 : 500 END DO
1796 0 : 510 END DO
1797 0 : 520 END DO
1798 0 : GO TO 580
1799 0 : 530 is = -ido
1800 0 : DO 570 j = 2,ip
1801 0 : is = is + ido
1802 0 : DO 560 k = 1,l1
1803 0 : idij = is
1804 0 : DO 550 i = 3,ido,2
1805 0 : idij = idij + 2
1806 0 : DO 540 m = 1,mp
1807 : c1(m,i-1,k,j) = wa(idij-1)*ch(m,i-1,k,j) - &
1808 0 : wa(idij)*ch(m,i,k,j)
1809 : c1(m,i,k,j) = wa(idij-1)*ch(m,i,k,j) + &
1810 0 : wa(idij)*ch(m,i-1,k,j)
1811 0 : 540 END DO
1812 0 : 550 END DO
1813 0 : 560 END DO
1814 0 : 570 END DO
1815 : 580 RETURN
1816 : END SUBROUTINE vradbg
1817 :
1818 0 : SUBROUTINE vrftb1(m,n,c,ch,mdimc,wa,fac)
1819 :
1820 : ! VRFFTPK, VERSION 1, AUGUST 1985
1821 :
1822 : IMPLICIT NONE
1823 : INTEGER, INTENT(IN) :: m, n, mdimc
1824 : REAL, INTENT(IN) :: wa(n), fac(15)
1825 : REAL, INTENT(INOUT) :: c(mdimc,n), ch(mdimc,n)
1826 :
1827 : INTEGER :: i, j, nf, na, l1, iw, k1, ip, ipo, idl1, ix2, ix3, ix4
1828 : INTEGER :: ido, l2
1829 : REAL :: scale
1830 :
1831 0 : nf = fac(2)
1832 0 : na = 0
1833 0 : l1 = 1
1834 0 : iw = 1
1835 0 : DO 160 k1 = 1,nf
1836 0 : ip = fac(k1+2)
1837 0 : l2 = ip*l1
1838 0 : ido = n/l2
1839 0 : idl1 = ido*l1
1840 0 : IF (ip /= 4) GO TO 30
1841 0 : ix2 = iw + ido
1842 0 : ix3 = ix2 + ido
1843 0 : IF (na /= 0) GO TO 10
1844 0 : CALL vradb4(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3))
1845 0 : GO TO 20
1846 0 : 10 CALL vradb4(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3))
1847 0 : 20 na = 1 - na
1848 0 : GO TO 150
1849 0 : 30 IF (ip /= 2) GO TO 60
1850 0 : IF (na /= 0) GO TO 40
1851 0 : CALL vradb2(m,ido,l1,c,ch,mdimc,wa(iw))
1852 0 : GO TO 50
1853 0 : 40 CALL vradb2(m,ido,l1,ch,c,mdimc,wa(iw))
1854 0 : 50 na = 1 - na
1855 0 : GO TO 150
1856 0 : 60 IF (ip /= 3) GO TO 90
1857 0 : ix2 = iw + ido
1858 0 : IF (na /= 0) GO TO 70
1859 0 : CALL vradb3(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2))
1860 0 : GO TO 80
1861 0 : 70 CALL vradb3(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2))
1862 0 : 80 na = 1 - na
1863 0 : GO TO 150
1864 0 : 90 IF (ip /= 5) GO TO 120
1865 0 : ix2 = iw + ido
1866 0 : ix3 = ix2 + ido
1867 0 : ix4 = ix3 + ido
1868 0 : IF (na /= 0) GO TO 100
1869 0 : CALL vradb5(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
1870 0 : GO TO 110
1871 0 : 100 CALL vradb5(m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
1872 0 : 110 na = 1 - na
1873 0 : GO TO 150
1874 0 : 120 IF (na /= 0) GO TO 130
1875 0 : CALL vradbg(m,ido,ip,l1,idl1,c,c,c,ch,ch,mdimc,wa(iw))
1876 0 : GO TO 140
1877 0 : 130 CALL vradbg(m,ido,ip,l1,idl1,ch,ch,ch,c,c,mdimc,wa(iw))
1878 0 : 140 IF (ido == 1) na = 1 - na
1879 0 : 150 l1 = l2
1880 0 : iw = iw + (ip-1)*ido
1881 0 : 160 END DO
1882 0 : scale = sqrt(1./n)
1883 0 : IF (na == 0) GO TO 190
1884 0 : DO 180 j = 1,n
1885 0 : DO 170 i = 1,m
1886 0 : c(i,j) = scale*ch(i,j)
1887 0 : 170 END DO
1888 0 : 180 END DO
1889 0 : RETURN
1890 0 : 190 DO 210 j = 1,n
1891 0 : DO 200 i = 1,m
1892 0 : c(i,j) = scale*c(i,j)
1893 0 : 200 END DO
1894 0 : 210 END DO
1895 : RETURN
1896 : END SUBROUTINE vrftb1
1897 :
1898 : END SUBROUTINE vrfftb
1899 :
1900 : END MODULE m_rfft
|
