NCEPLIBS-sp 2.4.0
fftpack.F
Go to the documentation of this file.
1C> @file
2C> @brief A concatination of the (FFTPACK)[https://netlib.org/fftpack/] library code.
3C>
4C> FFTPACK is a package of Fortran subprograms for the fast Fourier
5C> transform of periodic and other symmetric sequences. It includes
6C> complex, real, sine, cosine, and quarter-wave transforms.
7C>
8C>Reference:
9C>- P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations
10C>(G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83.
11C>
12C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
13
14C> dcrft
15C>
16C> @param init
17C> @param x
18C> @param ldx
19C> @param y
20C> @param ldy
21C> @param n
22C> @param m
23C> @param isign
24C> @param scale
25C> @param table
26C> @param n1
27C> @param wrk
28C> @param n2
29C> @param z
30C> @param nz
31C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
32 SUBROUTINE dcrft(init,x,ldx,y,ldy,n,m,isign,scale,
33 & table,n1,wrk,n2,z,nz)
34
35 implicit none
36 integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
37 real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk,z
38
39 IF (init.ne.0) THEN
40 CALL rffti(n,table)
41 ELSE
42!OCL NOVREC
43 DO j=1,m
44 y(1,j)=x(1,j)
45 DO i=2,n
46 y(i,j)=x(i+1,j)
47 ENDDO
48 CALL rfftb(n,y(1,j),table)
49 DO i=1,n
50 y(i,j)=scale*y(i,j)
51 ENDDO
52 ENDDO
53 ENDIF
54
55 RETURN
56 END
57
58C> scrft
59C>
60C> @param init
61C> @param x
62C> @param ldx
63C> @param y
64C> @param ldy
65C> @param n
66C> @param m
67C> @param isign
68C> @param scale
69C> @param table
70C> @param n1
71C> @param wrk
72C> @param n2
73C> @param z
74C> @param nz
75C>
76C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
77
78 SUBROUTINE scrft(init,x,ldx,y,ldy,n,m,isign,scale,
79 & table,n1,wrk,n2,z,nz)
80
81 implicit none
82 integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
83 real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk,z
84
85 IF (init.ne.0) THEN
86 CALL rffti(n,table)
87 ELSE
88!OCL NOVREC
89 DO j=1,m
90 y(1,j)=x(1,j)
91 DO i=2,n
92 y(i,j)=x(i+1,j)
93 ENDDO
94 CALL rfftb(n,y(1,j),table)
95 DO i=1,n
96 y(i,j)=scale*y(i,j)
97 ENDDO
98 ENDDO
99 ENDIF
100
101 RETURN
102 END
103
104C> csfft
105C>
106C> @param isign
107C> @param n
108C> @param scale
109C> @param x
110C> @param y
111C> @param table
112C> @param work
113C> @param isys
114C>
115C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
116 SUBROUTINE csfft(isign,n,scale,x,y,table,work,isys)
117
118 implicit none
119 integer isign,n,isys,i
120 real scale,x(*),y(*),table(*),work(*)
121
122 IF (isign.eq.0) THEN
123 CALL rffti(n,table)
124 ENDIF
125 IF (isign.eq.1) THEN
126 y(1)=x(1)
127 DO i=2,n
128 y(i)=x(i+1)
129 ENDDO
130 CALL rfftb(n,y,table)
131 DO i=1,n
132 y(i)=scale*y(i)
133 ENDDO
134 ENDIF
135
136 RETURN
137 END
138
139C> drcft
140C>
141C> @param init
142C> @param x
143C> @param ldx
144C> @param y
145C> @param ldy
146C> @param n
147C> @param m
148C> @param isign
149C> @param scale
150C> @param table
151C> @param n1
152C> @param wrk
153C> @param n2
154C> @param z
155C> @param nz
156C>
157C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
158 SUBROUTINE drcft(init,x,ldx,y,ldy,n,m,isign,scale,
159 & table,n1,wrk,n2,z,nz)
160
161 implicit none
162 integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
163 real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk,z
164
165 IF (init.ne.0) THEN
166 CALL rffti(n,table)
167 ELSE
168 DO j=1,m
169 DO i=1,n
170 y(i,j)=x(i,j)
171 ENDDO
172 CALL rfftf(n,y(1,j),table)
173 DO i=1,n
174 y(i,j)=scale*y(i,j)
175 ENDDO
176 DO i=n,2,-1
177 y(i+1,j)=y(i,j)
178 ENDDO
179 y(2,j)=0.
180C 01/17/2013 vvvvvvvvvvvvv E.Mirvis added ver 2.0.1 by S.Moorthi request. No +|- demo.
181 y(n+2,j) = 0.
182 ENDDO
183 ENDIF
184
185 RETURN
186 END
187
188C> srcft
189C>
190C> @param init
191C> @param x
192C> @param ldx
193C> @param y
194C> @param ldy
195C> @param n
196C> @param m
197C> @param isign
198C> @param scale
199C> @param table
200C> @param n1
201C> @param wrk
202C> @param n2
203C> @param z
204C> @param nz
205C>
206C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
207 SUBROUTINE srcft(init,x,ldx,y,ldy,n,m,isign,scale,
208 & table,n1,wrk,n2,z,nz)
209
210 implicit none
211 integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
212 real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk,z
213
214 IF (init.ne.0) THEN
215 CALL rffti(n,table)
216 ELSE
217 DO j=1,m
218 DO i=1,n
219 y(i,j)=x(i,j)
220 ENDDO
221 CALL rfftf(n,y(1,j),table)
222 DO i=1,n
223 y(i,j)=scale*y(i,j)
224 ENDDO
225 DO i=n,2,-1
226 y(i+1,j)=y(i,j)
227 ENDDO
228 y(2,j)=0.
229 y(n+2,j) = 0.
230C 01/17/2013 ^^^^^^^^^^E.Mirvis added ver 2.0.1 by S.Moorthi request. No +|- demo.
231 ENDDO
232 ENDIF
233
234 RETURN
235 END
236
237C> scfft
238C>
239C> @param isign
240C> @param n
241C> @param scale
242C> @param x
243C> @param y
244C> @param table
245C> @param work
246C> @param isys
247C>
248C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
249 SUBROUTINE scfft(isign,n,scale,x,y,table,work,isys)
250
251 implicit none
252 integer isign,n,isys,i
253 real scale,x(*),y(*),table(*),work(*)
254
255 IF (isign.eq.0) THEN
256 CALL rffti(n,table)
257 ENDIF
258 IF (isign.eq.-1) THEN
259 DO i=1,n
260 y(i)=x(i)
261 ENDDO
262 CALL rfftf(n,y,table)
263 DO i=1,n
264 y(i)=scale*y(i)
265 ENDDO
266 DO i=n,2,-1
267 y(i+1)=y(i)
268 ENDDO
269 y(2)=0.
270 ENDIF
271
272 RETURN
273 END
274
275C> RFFTF
276C>
277C> @param N
278C> @param R
279C> @param WSAVE
280C>
281C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
282 SUBROUTINE rfftf (N,R,WSAVE)
283 dimension r(1) ,wsave(1)
284 IF (n .EQ. 1) RETURN
285 CALL rfftf1 (n,r,wsave,wsave(n+1),wsave(2*n+1))
286 RETURN
287 END
288
289C> RFFTB
290C>
291C> @param N
292C> @param R
293C> @param WSAVE
294C>
295C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
296 SUBROUTINE rfftb (N,R,WSAVE)
297 dimension r(1) ,wsave(1)
298 IF (n .EQ. 1) RETURN
299 CALL rfftb1 (n,r,wsave,wsave(n+1),wsave(2*n+1))
300 RETURN
301 END
302
303C> RFFTI
304C>
305C> @param N
306C> @param WSAVE
307C>
308C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
309 SUBROUTINE rffti (N,WSAVE)
310 dimension wsave(1)
311 IF (n .EQ. 1) RETURN
312 CALL rffti1 (n,wsave(n+1),wsave(2*n+1))
313 RETURN
314 END
315
316C> RFFTB1
317C>
318C> @param N
319C> @param C
320C> @param CH
321C> @param WA
322C> @param IFAC
323C>
324C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
325 SUBROUTINE rfftb1 (N,C,CH,WA,IFAC)
326 dimension ch(1) ,c(1) ,wa(1) ,ifac(*)
327 nf = ifac(2)
328 na = 0
329 l1 = 1
330 iw = 1
331 DO 116 k1=1,nf
332 ip = ifac(k1+2)
333 l2 = ip*l1
334 ido = n/l2
335 idl1 = ido*l1
336 IF (ip .NE. 4) GO TO 103
337 ix2 = iw+ido
338 ix3 = ix2+ido
339 IF (na .NE. 0) GO TO 101
340 CALL radb4 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3))
341 GO TO 102
342 101 CALL radb4 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3))
343 102 na = 1-na
344 GO TO 115
345 103 IF (ip .NE. 2) GO TO 106
346 IF (na .NE. 0) GO TO 104
347 CALL radb2 (ido,l1,c,ch,wa(iw))
348 GO TO 105
349 104 CALL radb2 (ido,l1,ch,c,wa(iw))
350 105 na = 1-na
351 GO TO 115
352 106 IF (ip .NE. 3) GO TO 109
353 ix2 = iw+ido
354 IF (na .NE. 0) GO TO 107
355 CALL radb3 (ido,l1,c,ch,wa(iw),wa(ix2))
356 GO TO 108
357 107 CALL radb3 (ido,l1,ch,c,wa(iw),wa(ix2))
358 108 na = 1-na
359 GO TO 115
360 109 IF (ip .NE. 5) GO TO 112
361 ix2 = iw+ido
362 ix3 = ix2+ido
363 ix4 = ix3+ido
364 IF (na .NE. 0) GO TO 110
365 CALL radb5 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))
366 GO TO 111
367 110 CALL radb5 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))
368 111 na = 1-na
369 GO TO 115
370 112 IF (na .NE. 0) GO TO 113
371 CALL radbg (ido,ip,l1,idl1,c,c,c,ch,ch,wa(iw))
372 GO TO 114
373 113 CALL radbg (ido,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))
374 114 IF (ido .EQ. 1) na = 1-na
375 115 l1 = l2
376 iw = iw+(ip-1)*ido
377 116 CONTINUE
378 IF (na .EQ. 0) RETURN
379 DO 117 i=1,n
380 c(i) = ch(i)
381 117 CONTINUE
382 RETURN
383 END
384
385C> RFFTF1
386C>
387C> @param N
388C> @param C
389C> @param CH
390C> @param WA
391C> @param IFAC
392C>
393C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
394 SUBROUTINE rfftf1 (N,C,CH,WA,IFAC)
395 dimension ch(1) ,c(1) ,wa(1) ,ifac(*)
396 nf = ifac(2)
397 na = 1
398 l2 = n
399 iw = n
400 DO 111 k1=1,nf
401 kh = nf-k1
402 ip = ifac(kh+3)
403 l1 = l2/ip
404 ido = n/l2
405 idl1 = ido*l1
406 iw = iw-(ip-1)*ido
407 na = 1-na
408 IF (ip .NE. 4) GO TO 102
409 ix2 = iw+ido
410 ix3 = ix2+ido
411 IF (na .NE. 0) GO TO 101
412 CALL radf4 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3))
413 GO TO 110
414 101 CALL radf4 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3))
415 GO TO 110
416 102 IF (ip .NE. 2) GO TO 104
417 IF (na .NE. 0) GO TO 103
418 CALL radf2 (ido,l1,c,ch,wa(iw))
419 GO TO 110
420 103 CALL radf2 (ido,l1,ch,c,wa(iw))
421 GO TO 110
422 104 IF (ip .NE. 3) GO TO 106
423 ix2 = iw+ido
424 IF (na .NE. 0) GO TO 105
425 CALL radf3 (ido,l1,c,ch,wa(iw),wa(ix2))
426 GO TO 110
427 105 CALL radf3 (ido,l1,ch,c,wa(iw),wa(ix2))
428 GO TO 110
429 106 IF (ip .NE. 5) GO TO 108
430 ix2 = iw+ido
431 ix3 = ix2+ido
432 ix4 = ix3+ido
433 IF (na .NE. 0) GO TO 107
434 CALL radf5 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))
435 GO TO 110
436 107 CALL radf5 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))
437 GO TO 110
438 108 IF (ido .EQ. 1) na = 1-na
439 IF (na .NE. 0) GO TO 109
440 CALL radfg (ido,ip,l1,idl1,c,c,c,ch,ch,wa(iw))
441 na = 1
442 GO TO 110
443 109 CALL radfg (ido,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))
444 na = 0
445 110 l2 = l1
446 111 CONTINUE
447 IF (na .EQ. 1) RETURN
448 DO 112 i=1,n
449 c(i) = ch(i)
450 112 CONTINUE
451 RETURN
452 END
453
454C> RFFTI1
455C>
456C> @param N
457C> @param WA
458C> @param IFAC
459C>
460C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
461 SUBROUTINE rffti1 (N,WA,IFAC)
462 dimension wa(1) ,ifac(*) ,ntryh(4)
463 DATA ntryh(1),ntryh(2),ntryh(3),ntryh(4)/4,2,3,5/
464 nl = n
465 nf = 0
466 j = 0
467 101 j = j+1
468 IF (j-4) 102,102,103
469 102 ntry = ntryh(j)
470 GO TO 104
471 103 ntry = ntry+2
472 104 nq = nl/ntry
473 nr = nl-ntry*nq
474 IF (nr) 101,105,101
475 105 nf = nf+1
476 ifac(nf+2) = ntry
477 nl = nq
478 IF (ntry .NE. 2) GO TO 107
479 IF (nf .EQ. 1) GO TO 107
480 DO 106 i=2,nf
481 ib = nf-i+2
482 ifac(ib+2) = ifac(ib+1)
483 106 CONTINUE
484 ifac(3) = 2
485 107 IF (nl .NE. 1) GO TO 104
486 ifac(1) = n
487 ifac(2) = nf
488 tpi = 6.28318530717959
489 argh = tpi/float(n)
490 is = 0
491 nfm1 = nf-1
492 l1 = 1
493 IF (nfm1 .EQ. 0) RETURN
494!OCL NOVREC
495 DO 110 k1=1,nfm1
496 ip = ifac(k1+2)
497 ld = 0
498 l2 = l1*ip
499 ido = n/l2
500 ipm = ip-1
501 DO 109 j=1,ipm
502 ld = ld+l1
503 i = is
504 argld = float(ld)*argh
505 fi = 0
506!OCL SCALAR
507 DO 108 ii=3,ido,2
508 i = i+2
509 fi = fi+1
510 arg = fi*argld
511 wa(i-1) = cos(arg)
512 wa(i) = sin(arg)
513 108 CONTINUE
514 is = is+ido
515 109 CONTINUE
516 l1 = l2
517 110 CONTINUE
518 RETURN
519 END
520
521C> RADB2
522C>
523C> @param IDO
524C> @param L1
525C> @param CC
526C> @param CH
527C> @param WA1
528C>
529C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
530 SUBROUTINE radb2 (IDO,L1,CC,CH,WA1)
531 dimension cc(ido,2,l1) ,ch(ido,l1,2) ,
532 1 wa1(1)
533 DO 101 k=1,l1
534 ch(1,k,1) = cc(1,1,k)+cc(ido,2,k)
535 ch(1,k,2) = cc(1,1,k)-cc(ido,2,k)
536 101 CONTINUE
537 IF (ido-2) 107,105,102
538 102 idp2 = ido+2
539!OCL NOVREC
540 DO 104 k=1,l1
541 DO 103 i=3,ido,2
542 ic = idp2-i
543 ch(i-1,k,1) = cc(i-1,1,k)+cc(ic-1,2,k)
544 tr2 = cc(i-1,1,k)-cc(ic-1,2,k)
545 ch(i,k,1) = cc(i,1,k)-cc(ic,2,k)
546 ti2 = cc(i,1,k)+cc(ic,2,k)
547 ch(i-1,k,2) = wa1(i-2)*tr2-wa1(i-1)*ti2
548 ch(i,k,2) = wa1(i-2)*ti2+wa1(i-1)*tr2
549 103 CONTINUE
550 104 CONTINUE
551 IF (mod(ido,2) .EQ. 1) RETURN
552 105 DO 106 k=1,l1
553 ch(ido,k,1) = cc(ido,1,k)+cc(ido,1,k)
554 ch(ido,k,2) = -(cc(1,2,k)+cc(1,2,k))
555 106 CONTINUE
556 107 RETURN
557 END
558
559C> RADB3
560C>
561C> @param IDO
562C> @param L1
563C> @param CC
564C> @param CH
565C> @param WA1
566C> @param WA2
567C>
568C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
569 SUBROUTINE radb3 (IDO,L1,CC,CH,WA1,WA2)
570 dimension cc(ido,3,l1) ,ch(ido,l1,3) ,
571 1 wa1(1) ,wa2(1)
572 DATA taur,taui /-.5,.866025403784439/
573 DO 101 k=1,l1
574 tr2 = cc(ido,2,k)+cc(ido,2,k)
575 cr2 = cc(1,1,k)+taur*tr2
576 ch(1,k,1) = cc(1,1,k)+tr2
577 ci3 = taui*(cc(1,3,k)+cc(1,3,k))
578 ch(1,k,2) = cr2-ci3
579 ch(1,k,3) = cr2+ci3
580 101 CONTINUE
581 IF (ido .EQ. 1) RETURN
582 idp2 = ido+2
583!OCL NOVREC
584 DO 103 k=1,l1
585 DO 102 i=3,ido,2
586 ic = idp2-i
587 tr2 = cc(i-1,3,k)+cc(ic-1,2,k)
588 cr2 = cc(i-1,1,k)+taur*tr2
589 ch(i-1,k,1) = cc(i-1,1,k)+tr2
590 ti2 = cc(i,3,k)-cc(ic,2,k)
591 ci2 = cc(i,1,k)+taur*ti2
592 ch(i,k,1) = cc(i,1,k)+ti2
593 cr3 = taui*(cc(i-1,3,k)-cc(ic-1,2,k))
594 ci3 = taui*(cc(i,3,k)+cc(ic,2,k))
595 dr2 = cr2-ci3
596 dr3 = cr2+ci3
597 di2 = ci2+cr3
598 di3 = ci2-cr3
599 ch(i-1,k,2) = wa1(i-2)*dr2-wa1(i-1)*di2
600 ch(i,k,2) = wa1(i-2)*di2+wa1(i-1)*dr2
601 ch(i-1,k,3) = wa2(i-2)*dr3-wa2(i-1)*di3
602 ch(i,k,3) = wa2(i-2)*di3+wa2(i-1)*dr3
603 102 CONTINUE
604 103 CONTINUE
605 RETURN
606 END
607
608C> RADB4
609C>
610C> @param IDO
611C> @param L1
612C> @param CC
613C> @param CH
614C> @param WA1
615C> @param WA2
616C> @param WA3
617C>
618C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
619 SUBROUTINE radb4 (IDO,L1,CC,CH,WA1,WA2,WA3)
620 dimension cc(ido,4,l1) ,ch(ido,l1,4) ,
621 1 wa1(1) ,wa2(1) ,wa3(1)
622 DATA sqrt2 /1.414213562373095/
623 DO 101 k=1,l1
624 tr1 = cc(1,1,k)-cc(ido,4,k)
625 tr2 = cc(1,1,k)+cc(ido,4,k)
626 tr3 = cc(ido,2,k)+cc(ido,2,k)
627 tr4 = cc(1,3,k)+cc(1,3,k)
628 ch(1,k,1) = tr2+tr3
629 ch(1,k,2) = tr1-tr4
630 ch(1,k,3) = tr2-tr3
631 ch(1,k,4) = tr1+tr4
632 101 CONTINUE
633 IF (ido-2) 107,105,102
634 102 idp2 = ido+2
635!OCL NOVREC
636 DO 104 k=1,l1
637 DO 103 i=3,ido,2
638 ic = idp2-i
639 ti1 = cc(i,1,k)+cc(ic,4,k)
640 ti2 = cc(i,1,k)-cc(ic,4,k)
641 ti3 = cc(i,3,k)-cc(ic,2,k)
642 tr4 = cc(i,3,k)+cc(ic,2,k)
643 tr1 = cc(i-1,1,k)-cc(ic-1,4,k)
644 tr2 = cc(i-1,1,k)+cc(ic-1,4,k)
645 ti4 = cc(i-1,3,k)-cc(ic-1,2,k)
646 tr3 = cc(i-1,3,k)+cc(ic-1,2,k)
647 ch(i-1,k,1) = tr2+tr3
648 cr3 = tr2-tr3
649 ch(i,k,1) = ti2+ti3
650 ci3 = ti2-ti3
651 cr2 = tr1-tr4
652 cr4 = tr1+tr4
653 ci2 = ti1+ti4
654 ci4 = ti1-ti4
655 ch(i-1,k,2) = wa1(i-2)*cr2-wa1(i-1)*ci2
656 ch(i,k,2) = wa1(i-2)*ci2+wa1(i-1)*cr2
657 ch(i-1,k,3) = wa2(i-2)*cr3-wa2(i-1)*ci3
658 ch(i,k,3) = wa2(i-2)*ci3+wa2(i-1)*cr3
659 ch(i-1,k,4) = wa3(i-2)*cr4-wa3(i-1)*ci4
660 ch(i,k,4) = wa3(i-2)*ci4+wa3(i-1)*cr4
661 103 CONTINUE
662 104 CONTINUE
663 IF (mod(ido,2) .EQ. 1) RETURN
664 105 CONTINUE
665 DO 106 k=1,l1
666 ti1 = cc(1,2,k)+cc(1,4,k)
667 ti2 = cc(1,4,k)-cc(1,2,k)
668 tr1 = cc(ido,1,k)-cc(ido,3,k)
669 tr2 = cc(ido,1,k)+cc(ido,3,k)
670 ch(ido,k,1) = tr2+tr2
671 ch(ido,k,2) = sqrt2*(tr1-ti1)
672 ch(ido,k,3) = ti2+ti2
673 ch(ido,k,4) = -sqrt2*(tr1+ti1)
674 106 CONTINUE
675 107 RETURN
676 END
677
678C> RADB5
679C>
680C> @param IDO
681C> @param L1
682C> @param CC
683C> @param CH
684C> @param WA1
685C> @param WA2
686C> @param WA3
687C> @param WA4
688C>
689C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
690 SUBROUTINE radb5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
691 dimension cc(ido,5,l1) ,ch(ido,l1,5) ,
692 1 wa1(1) ,wa2(1) ,wa3(1) ,wa4(1)
693 DATA tr11,ti11,tr12,ti12 /.309016994374947,.951056516295154,
694 1-.809016994374947,.587785252292473/
695 DO 101 k=1,l1
696 ti5 = cc(1,3,k)+cc(1,3,k)
697 ti4 = cc(1,5,k)+cc(1,5,k)
698 tr2 = cc(ido,2,k)+cc(ido,2,k)
699 tr3 = cc(ido,4,k)+cc(ido,4,k)
700 ch(1,k,1) = cc(1,1,k)+tr2+tr3
701 cr2 = cc(1,1,k)+tr11*tr2+tr12*tr3
702 cr3 = cc(1,1,k)+tr12*tr2+tr11*tr3
703 ci5 = ti11*ti5+ti12*ti4
704 ci4 = ti12*ti5-ti11*ti4
705 ch(1,k,2) = cr2-ci5
706 ch(1,k,3) = cr3-ci4
707 ch(1,k,4) = cr3+ci4
708 ch(1,k,5) = cr2+ci5
709 101 CONTINUE
710 IF (ido .EQ. 1) RETURN
711 idp2 = ido+2
712 DO 103 k=1,l1
713 DO 102 i=3,ido,2
714 ic = idp2-i
715 ti5 = cc(i,3,k)+cc(ic,2,k)
716 ti2 = cc(i,3,k)-cc(ic,2,k)
717 ti4 = cc(i,5,k)+cc(ic,4,k)
718 ti3 = cc(i,5,k)-cc(ic,4,k)
719 tr5 = cc(i-1,3,k)-cc(ic-1,2,k)
720 tr2 = cc(i-1,3,k)+cc(ic-1,2,k)
721 tr4 = cc(i-1,5,k)-cc(ic-1,4,k)
722 tr3 = cc(i-1,5,k)+cc(ic-1,4,k)
723 ch(i-1,k,1) = cc(i-1,1,k)+tr2+tr3
724 ch(i,k,1) = cc(i,1,k)+ti2+ti3
725 cr2 = cc(i-1,1,k)+tr11*tr2+tr12*tr3
726 ci2 = cc(i,1,k)+tr11*ti2+tr12*ti3
727 cr3 = cc(i-1,1,k)+tr12*tr2+tr11*tr3
728 ci3 = cc(i,1,k)+tr12*ti2+tr11*ti3
729 cr5 = ti11*tr5+ti12*tr4
730 ci5 = ti11*ti5+ti12*ti4
731 cr4 = ti12*tr5-ti11*tr4
732 ci4 = ti12*ti5-ti11*ti4
733 dr3 = cr3-ci4
734 dr4 = cr3+ci4
735 di3 = ci3+cr4
736 di4 = ci3-cr4
737 dr5 = cr2+ci5
738 dr2 = cr2-ci5
739 di5 = ci2-cr5
740 di2 = ci2+cr5
741 ch(i-1,k,2) = wa1(i-2)*dr2-wa1(i-1)*di2
742 ch(i,k,2) = wa1(i-2)*di2+wa1(i-1)*dr2
743 ch(i-1,k,3) = wa2(i-2)*dr3-wa2(i-1)*di3
744 ch(i,k,3) = wa2(i-2)*di3+wa2(i-1)*dr3
745 ch(i-1,k,4) = wa3(i-2)*dr4-wa3(i-1)*di4
746 ch(i,k,4) = wa3(i-2)*di4+wa3(i-1)*dr4
747 ch(i-1,k,5) = wa4(i-2)*dr5-wa4(i-1)*di5
748 ch(i,k,5) = wa4(i-2)*di5+wa4(i-1)*dr5
749 102 CONTINUE
750 103 CONTINUE
751 RETURN
752 END
753
754C> RADBG
755C>
756C> @param IDO
757C> @param IP
758C> @param L1
759C> @param IDL1
760C> @param CC
761C> @param C1
762C> @param C2
763C> @param CH
764C> @param CH2
765C> @param WA
766C>
767C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
768 SUBROUTINE radbg (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
769 dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
770 1 c1(ido,l1,ip) ,c2(idl1,ip),
771 2 ch2(idl1,ip) ,wa(1)
772 DATA tpi/6.28318530717959/
773 arg = tpi/float(ip)
774 dcp = cos(arg)
775 dsp = sin(arg)
776 idp2 = ido+2
777 nbd = (ido-1)/2
778 ipp2 = ip+2
779 ipph = (ip+1)/2
780 IF (ido .LT. l1) GO TO 103
781 DO 102 k=1,l1
782 DO 101 i=1,ido
783 ch(i,k,1) = cc(i,1,k)
784 101 CONTINUE
785 102 CONTINUE
786 GO TO 106
787 103 DO 105 i=1,ido
788 DO 104 k=1,l1
789 ch(i,k,1) = cc(i,1,k)
790 104 CONTINUE
791 105 CONTINUE
792!OCL NOVREC
793 106 DO 108 j=2,ipph
794 jc = ipp2-j
795 j2 = j+j
796 DO 107 k=1,l1
797 ch(1,k,j) = cc(ido,j2-2,k)+cc(ido,j2-2,k)
798 ch(1,k,jc) = cc(1,j2-1,k)+cc(1,j2-1,k)
799 107 CONTINUE
800 108 CONTINUE
801 IF (ido .EQ. 1) GO TO 116
802 IF (nbd .LT. l1) GO TO 112
803!OCL NOVREC
804 DO 111 j=2,ipph
805 jc = ipp2-j
806 DO 110 k=1,l1
807 DO 109 i=3,ido,2
808 ic = idp2-i
809 ch(i-1,k,j) = cc(i-1,2*j-1,k)+cc(ic-1,2*j-2,k)
810 ch(i-1,k,jc) = cc(i-1,2*j-1,k)-cc(ic-1,2*j-2,k)
811 ch(i,k,j) = cc(i,2*j-1,k)-cc(ic,2*j-2,k)
812 ch(i,k,jc) = cc(i,2*j-1,k)+cc(ic,2*j-2,k)
813 109 CONTINUE
814 110 CONTINUE
815 111 CONTINUE
816 GO TO 116
817 112 DO 115 j=2,ipph
818 jc = ipp2-j
819 DO 114 i=3,ido,2
820 ic = idp2-i
821 DO 113 k=1,l1
822 ch(i-1,k,j) = cc(i-1,2*j-1,k)+cc(ic-1,2*j-2,k)
823 ch(i-1,k,jc) = cc(i-1,2*j-1,k)-cc(ic-1,2*j-2,k)
824 ch(i,k,j) = cc(i,2*j-1,k)-cc(ic,2*j-2,k)
825 ch(i,k,jc) = cc(i,2*j-1,k)+cc(ic,2*j-2,k)
826 113 CONTINUE
827 114 CONTINUE
828 115 CONTINUE
829 116 ar1 = 1.
830 ai1 = 0.
831!OCL NOVREC
832 DO 120 l=2,ipph
833 lc = ipp2-l
834 ar1h = dcp*ar1-dsp*ai1
835 ai1 = dcp*ai1+dsp*ar1
836 ar1 = ar1h
837 DO 117 ik=1,idl1
838 c2(ik,l) = ch2(ik,1)+ar1*ch2(ik,2)
839 c2(ik,lc) = ai1*ch2(ik,ip)
840 117 CONTINUE
841 dc2 = ar1
842 ds2 = ai1
843 ar2 = ar1
844 ai2 = ai1
845!OCL NOVREC
846 DO 119 j=3,ipph
847 jc = ipp2-j
848 ar2h = dc2*ar2-ds2*ai2
849 ai2 = dc2*ai2+ds2*ar2
850 ar2 = ar2h
851 DO 118 ik=1,idl1
852 c2(ik,l) = c2(ik,l)+ar2*ch2(ik,j)
853 c2(ik,lc) = c2(ik,lc)+ai2*ch2(ik,jc)
854 118 CONTINUE
855 119 CONTINUE
856 120 CONTINUE
857!OCL NOVREC
858 DO 122 j=2,ipph
859 DO 121 ik=1,idl1
860 ch2(ik,1) = ch2(ik,1)+ch2(ik,j)
861 121 CONTINUE
862 122 CONTINUE
863!OCL NOVREC
864 DO 124 j=2,ipph
865 jc = ipp2-j
866 DO 123 k=1,l1
867 ch(1,k,j) = c1(1,k,j)-c1(1,k,jc)
868 ch(1,k,jc) = c1(1,k,j)+c1(1,k,jc)
869 123 CONTINUE
870 124 CONTINUE
871 IF (ido .EQ. 1) GO TO 132
872 IF (nbd .LT. l1) GO TO 128
873!OCL NOVREC
874 DO 127 j=2,ipph
875 jc = ipp2-j
876 DO 126 k=1,l1
877 DO 125 i=3,ido,2
878 ch(i-1,k,j) = c1(i-1,k,j)-c1(i,k,jc)
879 ch(i-1,k,jc) = c1(i-1,k,j)+c1(i,k,jc)
880 ch(i,k,j) = c1(i,k,j)+c1(i-1,k,jc)
881 ch(i,k,jc) = c1(i,k,j)-c1(i-1,k,jc)
882 125 CONTINUE
883 126 CONTINUE
884 127 CONTINUE
885 GO TO 132
886 128 DO 131 j=2,ipph
887 jc = ipp2-j
888 DO 130 i=3,ido,2
889 DO 129 k=1,l1
890 ch(i-1,k,j) = c1(i-1,k,j)-c1(i,k,jc)
891 ch(i-1,k,jc) = c1(i-1,k,j)+c1(i,k,jc)
892 ch(i,k,j) = c1(i,k,j)+c1(i-1,k,jc)
893 ch(i,k,jc) = c1(i,k,j)-c1(i-1,k,jc)
894 129 CONTINUE
895 130 CONTINUE
896 131 CONTINUE
897 132 CONTINUE
898 IF (ido .EQ. 1) RETURN
899 DO 133 ik=1,idl1
900 c2(ik,1) = ch2(ik,1)
901 133 CONTINUE
902 DO 135 j=2,ip
903 DO 134 k=1,l1
904 c1(1,k,j) = ch(1,k,j)
905 134 CONTINUE
906 135 CONTINUE
907 IF (nbd .GT. l1) GO TO 139
908 is = -ido
909 DO 138 j=2,ip
910 is = is+ido
911 idij = is
912 DO 137 i=3,ido,2
913 idij = idij+2
914 DO 136 k=1,l1
915 c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
916 c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
917 136 CONTINUE
918 137 CONTINUE
919 138 CONTINUE
920 GO TO 143
921 139 is = -ido
922!OCL NOVREC
923 DO 142 j=2,ip
924 is = is+ido
925 DO 141 k=1,l1
926 idij = is
927 DO 140 i=3,ido,2
928 idij = idij+2
929 c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
930 c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
931 140 CONTINUE
932 141 CONTINUE
933 142 CONTINUE
934 143 RETURN
935 END
936
937C> RADBG
938C>
939C> @param IDO
940C> @param L1
941C> @param CC
942C> @param CH
943C> @param WA1
944C>
945C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
946 SUBROUTINE radf2 (IDO,L1,CC,CH,WA1)
947 dimension ch(ido,2,l1) ,cc(ido,l1,2) ,
948 1 wa1(1)
949 DO 101 k=1,l1
950 ch(1,1,k) = cc(1,k,1)+cc(1,k,2)
951 ch(ido,2,k) = cc(1,k,1)-cc(1,k,2)
952 101 CONTINUE
953 IF (ido-2) 107,105,102
954 102 idp2 = ido+2
955 DO 104 k=1,l1
956 DO 103 i=3,ido,2
957 ic = idp2-i
958 tr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
959 ti2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
960 ch(i,1,k) = cc(i,k,1)+ti2
961 ch(ic,2,k) = ti2-cc(i,k,1)
962 ch(i-1,1,k) = cc(i-1,k,1)+tr2
963 ch(ic-1,2,k) = cc(i-1,k,1)-tr2
964 103 CONTINUE
965 104 CONTINUE
966 IF (mod(ido,2) .EQ. 1) RETURN
967 105 DO 106 k=1,l1
968 ch(1,2,k) = -cc(ido,k,2)
969 ch(ido,1,k) = cc(ido,k,1)
970 106 CONTINUE
971 107 RETURN
972 END
973
974C> RADF3
975C>
976C> @param IDO
977C> @param L1
978C> @param CC
979C> @param CH
980C> @param WA1
981C> @param WA2
982C>
983C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
984 SUBROUTINE radf3 (IDO,L1,CC,CH,WA1,WA2)
985 dimension ch(ido,3,l1) ,cc(ido,l1,3) ,
986 1 wa1(1) ,wa2(1)
987 DATA taur,taui /-.5,.866025403784439/
988 DO 101 k=1,l1
989 cr2 = cc(1,k,2)+cc(1,k,3)
990 ch(1,1,k) = cc(1,k,1)+cr2
991 ch(1,3,k) = taui*(cc(1,k,3)-cc(1,k,2))
992 ch(ido,2,k) = cc(1,k,1)+taur*cr2
993 101 CONTINUE
994 IF (ido .EQ. 1) RETURN
995 idp2 = ido+2
996 DO 103 k=1,l1
997 DO 102 i=3,ido,2
998 ic = idp2-i
999 dr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1000 di2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1001 dr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1002 di3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1003 cr2 = dr2+dr3
1004 ci2 = di2+di3
1005 ch(i-1,1,k) = cc(i-1,k,1)+cr2
1006 ch(i,1,k) = cc(i,k,1)+ci2
1007 tr2 = cc(i-1,k,1)+taur*cr2
1008 ti2 = cc(i,k,1)+taur*ci2
1009 tr3 = taui*(di2-di3)
1010 ti3 = taui*(dr3-dr2)
1011 ch(i-1,3,k) = tr2+tr3
1012 ch(ic-1,2,k) = tr2-tr3
1013 ch(i,3,k) = ti2+ti3
1014 ch(ic,2,k) = ti3-ti2
1015 102 CONTINUE
1016 103 CONTINUE
1017 RETURN
1018 END
1019
1020C> RADF4
1021C>
1022C> @param IDO
1023C> @param L1
1024C> @param CC
1025C> @param CH
1026C> @param WA1
1027C> @param WA2
1028C> @param WA3
1029C>
1030C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1031 SUBROUTINE radf4 (IDO,L1,CC,CH,WA1,WA2,WA3)
1032 dimension cc(ido,l1,4) ,ch(ido,4,l1) ,
1033 1 wa1(1) ,wa2(1) ,wa3(1)
1034 DATA hsqt2 /.7071067811865475/
1035 DO 101 k=1,l1
1036 tr1 = cc(1,k,2)+cc(1,k,4)
1037 tr2 = cc(1,k,1)+cc(1,k,3)
1038 ch(1,1,k) = tr1+tr2
1039 ch(ido,4,k) = tr2-tr1
1040 ch(ido,2,k) = cc(1,k,1)-cc(1,k,3)
1041 ch(1,3,k) = cc(1,k,4)-cc(1,k,2)
1042 101 CONTINUE
1043 IF (ido-2) 107,105,102
1044 102 idp2 = ido+2
1045!OCL NOVREC
1046 DO 104 k=1,l1
1047 DO 103 i=3,ido,2
1048 ic = idp2-i
1049 cr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1050 ci2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1051 cr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1052 ci3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1053 cr4 = wa3(i-2)*cc(i-1,k,4)+wa3(i-1)*cc(i,k,4)
1054 ci4 = wa3(i-2)*cc(i,k,4)-wa3(i-1)*cc(i-1,k,4)
1055 tr1 = cr2+cr4
1056 tr4 = cr4-cr2
1057 ti1 = ci2+ci4
1058 ti4 = ci2-ci4
1059 ti2 = cc(i,k,1)+ci3
1060 ti3 = cc(i,k,1)-ci3
1061 tr2 = cc(i-1,k,1)+cr3
1062 tr3 = cc(i-1,k,1)-cr3
1063 ch(i-1,1,k) = tr1+tr2
1064 ch(ic-1,4,k) = tr2-tr1
1065 ch(i,1,k) = ti1+ti2
1066 ch(ic,4,k) = ti1-ti2
1067 ch(i-1,3,k) = ti4+tr3
1068 ch(ic-1,2,k) = tr3-ti4
1069 ch(i,3,k) = tr4+ti3
1070 ch(ic,2,k) = tr4-ti3
1071 103 CONTINUE
1072 104 CONTINUE
1073 IF (mod(ido,2) .EQ. 1) RETURN
1074 105 CONTINUE
1075 DO 106 k=1,l1
1076 ti1 = -hsqt2*(cc(ido,k,2)+cc(ido,k,4))
1077 tr1 = hsqt2*(cc(ido,k,2)-cc(ido,k,4))
1078 ch(ido,1,k) = tr1+cc(ido,k,1)
1079 ch(ido,3,k) = cc(ido,k,1)-tr1
1080 ch(1,2,k) = ti1-cc(ido,k,3)
1081 ch(1,4,k) = ti1+cc(ido,k,3)
1082 106 CONTINUE
1083 107 RETURN
1084 END
1085
1086C> RADF5
1087C>
1088C> @param IDO
1089C> @param L1
1090C> @param CC
1091C> @param CH
1092C> @param WA1
1093C> @param WA2
1094C> @param WA3
1095C> @param WA4
1096C>
1097C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1098 SUBROUTINE radf5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1099 dimension cc(ido,l1,5) ,ch(ido,5,l1) ,
1100 1 wa1(1) ,wa2(1) ,wa3(1) ,wa4(1)
1101 DATA tr11,ti11,tr12,ti12 /.309016994374947,.951056516295154,
1102 1-.809016994374947,.587785252292473/
1103 DO 101 k=1,l1
1104 cr2 = cc(1,k,5)+cc(1,k,2)
1105 ci5 = cc(1,k,5)-cc(1,k,2)
1106 cr3 = cc(1,k,4)+cc(1,k,3)
1107 ci4 = cc(1,k,4)-cc(1,k,3)
1108 ch(1,1,k) = cc(1,k,1)+cr2+cr3
1109 ch(ido,2,k) = cc(1,k,1)+tr11*cr2+tr12*cr3
1110 ch(1,3,k) = ti11*ci5+ti12*ci4
1111 ch(ido,4,k) = cc(1,k,1)+tr12*cr2+tr11*cr3
1112 ch(1,5,k) = ti12*ci5-ti11*ci4
1113 101 CONTINUE
1114 IF (ido .EQ. 1) RETURN
1115 idp2 = ido+2
1116 DO 103 k=1,l1
1117 DO 102 i=3,ido,2
1118 ic = idp2-i
1119 dr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1120 di2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1121 dr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1122 di3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1123 dr4 = wa3(i-2)*cc(i-1,k,4)+wa3(i-1)*cc(i,k,4)
1124 di4 = wa3(i-2)*cc(i,k,4)-wa3(i-1)*cc(i-1,k,4)
1125 dr5 = wa4(i-2)*cc(i-1,k,5)+wa4(i-1)*cc(i,k,5)
1126 di5 = wa4(i-2)*cc(i,k,5)-wa4(i-1)*cc(i-1,k,5)
1127 cr2 = dr2+dr5
1128 ci5 = dr5-dr2
1129 cr5 = di2-di5
1130 ci2 = di2+di5
1131 cr3 = dr3+dr4
1132 ci4 = dr4-dr3
1133 cr4 = di3-di4
1134 ci3 = di3+di4
1135 ch(i-1,1,k) = cc(i-1,k,1)+cr2+cr3
1136 ch(i,1,k) = cc(i,k,1)+ci2+ci3
1137 tr2 = cc(i-1,k,1)+tr11*cr2+tr12*cr3
1138 ti2 = cc(i,k,1)+tr11*ci2+tr12*ci3
1139 tr3 = cc(i-1,k,1)+tr12*cr2+tr11*cr3
1140 ti3 = cc(i,k,1)+tr12*ci2+tr11*ci3
1141 tr5 = ti11*cr5+ti12*cr4
1142 ti5 = ti11*ci5+ti12*ci4
1143 tr4 = ti12*cr5-ti11*cr4
1144 ti4 = ti12*ci5-ti11*ci4
1145 ch(i-1,3,k) = tr2+tr5
1146 ch(ic-1,2,k) = tr2-tr5
1147 ch(i,3,k) = ti2+ti5
1148 ch(ic,2,k) = ti5-ti2
1149 ch(i-1,5,k) = tr3+tr4
1150 ch(ic-1,4,k) = tr3-tr4
1151 ch(i,5,k) = ti3+ti4
1152 ch(ic,4,k) = ti4-ti3
1153 102 CONTINUE
1154 103 CONTINUE
1155 RETURN
1156 END
1157
1158C> RADFG
1159C>
1160C> @param IDO
1161C> @param IP
1162C> @param L1
1163C> @param IDL1
1164C> @param CC
1165C> @param C1
1166C> @param C2
1167C> @param CH
1168C> @param CH2
1169C> @param WA
1170C>
1171C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1172 SUBROUTINE radfg (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1173 dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
1174 1 c1(ido,l1,ip) ,c2(idl1,ip),
1175 2 ch2(idl1,ip) ,wa(1)
1176 DATA tpi/6.28318530717959/
1177 arg = tpi/float(ip)
1178 dcp = cos(arg)
1179 dsp = sin(arg)
1180 ipph = (ip+1)/2
1181 ipp2 = ip+2
1182 idp2 = ido+2
1183 nbd = (ido-1)/2
1184 IF (ido .EQ. 1) GO TO 119
1185 DO 101 ik=1,idl1
1186 ch2(ik,1) = c2(ik,1)
1187 101 CONTINUE
1188 DO 103 j=2,ip
1189 DO 102 k=1,l1
1190 ch(1,k,j) = c1(1,k,j)
1191 102 CONTINUE
1192 103 CONTINUE
1193 IF (nbd .GT. l1) GO TO 107
1194 is = -ido
1195 DO 106 j=2,ip
1196 is = is+ido
1197 idij = is
1198 DO 105 i=3,ido,2
1199 idij = idij+2
1200 DO 104 k=1,l1
1201 ch(i-1,k,j) = wa(idij-1)*c1(i-1,k,j)+wa(idij)*c1(i,k,j)
1202 ch(i,k,j) = wa(idij-1)*c1(i,k,j)-wa(idij)*c1(i-1,k,j)
1203 104 CONTINUE
1204 105 CONTINUE
1205 106 CONTINUE
1206 GO TO 111
1207 107 is = -ido
1208 DO 110 j=2,ip
1209 is = is+ido
1210 DO 109 k=1,l1
1211 idij = is
1212 DO 108 i=3,ido,2
1213 idij = idij+2
1214 ch(i-1,k,j) = wa(idij-1)*c1(i-1,k,j)+wa(idij)*c1(i,k,j)
1215 ch(i,k,j) = wa(idij-1)*c1(i,k,j)-wa(idij)*c1(i-1,k,j)
1216 108 CONTINUE
1217 109 CONTINUE
1218 110 CONTINUE
1219 111 IF (nbd .LT. l1) GO TO 115
1220 DO 114 j=2,ipph
1221 jc = ipp2-j
1222 DO 113 k=1,l1
1223 DO 112 i=3,ido,2
1224 c1(i-1,k,j) = ch(i-1,k,j)+ch(i-1,k,jc)
1225 c1(i-1,k,jc) = ch(i,k,j)-ch(i,k,jc)
1226 c1(i,k,j) = ch(i,k,j)+ch(i,k,jc)
1227 c1(i,k,jc) = ch(i-1,k,jc)-ch(i-1,k,j)
1228 112 CONTINUE
1229 113 CONTINUE
1230 114 CONTINUE
1231 GO TO 121
1232 115 DO 118 j=2,ipph
1233 jc = ipp2-j
1234 DO 117 i=3,ido,2
1235 DO 116 k=1,l1
1236 c1(i-1,k,j) = ch(i-1,k,j)+ch(i-1,k,jc)
1237 c1(i-1,k,jc) = ch(i,k,j)-ch(i,k,jc)
1238 c1(i,k,j) = ch(i,k,j)+ch(i,k,jc)
1239 c1(i,k,jc) = ch(i-1,k,jc)-ch(i-1,k,j)
1240 116 CONTINUE
1241 117 CONTINUE
1242 118 CONTINUE
1243 GO TO 121
1244 119 DO 120 ik=1,idl1
1245 c2(ik,1) = ch2(ik,1)
1246 120 CONTINUE
1247 121 DO 123 j=2,ipph
1248 jc = ipp2-j
1249 DO 122 k=1,l1
1250 c1(1,k,j) = ch(1,k,j)+ch(1,k,jc)
1251 c1(1,k,jc) = ch(1,k,jc)-ch(1,k,j)
1252 122 CONTINUE
1253 123 CONTINUE
1254C
1255 ar1 = 1.
1256 ai1 = 0.
1257 DO 127 l=2,ipph
1258 lc = ipp2-l
1259 ar1h = dcp*ar1-dsp*ai1
1260 ai1 = dcp*ai1+dsp*ar1
1261 ar1 = ar1h
1262 DO 124 ik=1,idl1
1263 ch2(ik,l) = c2(ik,1)+ar1*c2(ik,2)
1264 ch2(ik,lc) = ai1*c2(ik,ip)
1265 124 CONTINUE
1266 dc2 = ar1
1267 ds2 = ai1
1268 ar2 = ar1
1269 ai2 = ai1
1270 DO 126 j=3,ipph
1271 jc = ipp2-j
1272 ar2h = dc2*ar2-ds2*ai2
1273 ai2 = dc2*ai2+ds2*ar2
1274 ar2 = ar2h
1275 DO 125 ik=1,idl1
1276 ch2(ik,l) = ch2(ik,l)+ar2*c2(ik,j)
1277 ch2(ik,lc) = ch2(ik,lc)+ai2*c2(ik,jc)
1278 125 CONTINUE
1279 126 CONTINUE
1280 127 CONTINUE
1281 DO 129 j=2,ipph
1282 DO 128 ik=1,idl1
1283 ch2(ik,1) = ch2(ik,1)+c2(ik,j)
1284 128 CONTINUE
1285 129 CONTINUE
1286C
1287 IF (ido .LT. l1) GO TO 132
1288 DO 131 k=1,l1
1289 DO 130 i=1,ido
1290 cc(i,1,k) = ch(i,k,1)
1291 130 CONTINUE
1292 131 CONTINUE
1293 GO TO 135
1294 132 DO 134 i=1,ido
1295 DO 133 k=1,l1
1296 cc(i,1,k) = ch(i,k,1)
1297 133 CONTINUE
1298 134 CONTINUE
1299 135 DO 137 j=2,ipph
1300 jc = ipp2-j
1301 j2 = j+j
1302 DO 136 k=1,l1
1303 cc(ido,j2-2,k) = ch(1,k,j)
1304 cc(1,j2-1,k) = ch(1,k,jc)
1305 136 CONTINUE
1306 137 CONTINUE
1307 IF (ido .EQ. 1) RETURN
1308 IF (nbd .LT. l1) GO TO 141
1309 DO 140 j=2,ipph
1310 jc = ipp2-j
1311 j2 = j+j
1312 DO 139 k=1,l1
1313 DO 138 i=3,ido,2
1314 ic = idp2-i
1315 cc(i-1,j2-1,k) = ch(i-1,k,j)+ch(i-1,k,jc)
1316 cc(ic-1,j2-2,k) = ch(i-1,k,j)-ch(i-1,k,jc)
1317 cc(i,j2-1,k) = ch(i,k,j)+ch(i,k,jc)
1318 cc(ic,j2-2,k) = ch(i,k,jc)-ch(i,k,j)
1319 138 CONTINUE
1320 139 CONTINUE
1321 140 CONTINUE
1322 RETURN
1323 141 DO 144 j=2,ipph
1324 jc = ipp2-j
1325 j2 = j+j
1326 DO 143 i=3,ido,2
1327 ic = idp2-i
1328 DO 142 k=1,l1
1329 cc(i-1,j2-1,k) = ch(i-1,k,j)+ch(i-1,k,jc)
1330 cc(ic-1,j2-2,k) = ch(i-1,k,j)-ch(i-1,k,jc)
1331 cc(i,j2-1,k) = ch(i,k,j)+ch(i,k,jc)
1332 cc(ic,j2-2,k) = ch(i,k,jc)-ch(i,k,j)
1333 142 CONTINUE
1334 143 CONTINUE
1335 144 CONTINUE
1336 RETURN
1337 END
radb5
subroutine radb5(IDO, L1, CC, CH, WA1, WA2, WA3, WA4)
RADB5.
Definition: fftpack.F:691
scfft
subroutine scfft(isign, n, scale, x, y, table, work, isys)
scfft
Definition: fftpack.F:250
radb4
subroutine radb4(IDO, L1, CC, CH, WA1, WA2, WA3)
RADB4.
Definition: fftpack.F:620
rffti
subroutine rffti(N, WSAVE)
RFFTI.
Definition: fftpack.F:310
radf3
subroutine radf3(IDO, L1, CC, CH, WA1, WA2)
RADF3.
Definition: fftpack.F:985
drcft
subroutine drcft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
drcft
Definition: fftpack.F:160
radb3
subroutine radb3(IDO, L1, CC, CH, WA1, WA2)
RADB3
Definition: fftpack.F:570
radf5
subroutine radf5(IDO, L1, CC, CH, WA1, WA2, WA3, WA4)
RADF5.
Definition: fftpack.F:1099
radf2
subroutine radf2(IDO, L1, CC, CH, WA1)
RADBG.
Definition: fftpack.F:947
dcrft
subroutine dcrft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
dcrft
Definition: fftpack.F:34
rfftb
subroutine rfftb(N, R, WSAVE)
RFFTB.
Definition: fftpack.F:297
scrft
subroutine scrft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
scrft
Definition: fftpack.F:80
rfftf1
subroutine rfftf1(N, C, CH, WA, IFAC)
RFFTF1.
Definition: fftpack.F:395
srcft
subroutine srcft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
srcft
Definition: fftpack.F:209
csfft
subroutine csfft(isign, n, scale, x, y, table, work, isys)
csfft
Definition: fftpack.F:117
radfg
subroutine radfg(IDO, IP, L1, IDL1, CC, C1, C2, CH, CH2, WA)
RADFG.
Definition: fftpack.F:1173
rffti1
subroutine rffti1(N, WA, IFAC)
RFFTI1.
Definition: fftpack.F:462
radb2
subroutine radb2(IDO, L1, CC, CH, WA1)
RADB2.
Definition: fftpack.F:531
radf4
subroutine radf4(IDO, L1, CC, CH, WA1, WA2, WA3)
RADF4.
Definition: fftpack.F:1032
radbg
subroutine radbg(IDO, IP, L1, IDL1, CC, C1, C2, CH, CH2, WA)
RADBG.
Definition: fftpack.F:769
rfftb1
subroutine rfftb1(N, C, CH, WA, IFAC)
RFFTB1.
Definition: fftpack.F:326
rfftf
subroutine rfftf(N, R, WSAVE)
RFFTF.
Definition: fftpack.F:283