NCEPLIBS-sp  2.3.3
•All Files Functions
fftpack.F
1 #if defined LINUX || defined APPLE
2  SUBROUTINE dcrft(init,x,ldx,y,ldy,n,m,isign,scale,
3  & table,n1,wrk,n2,z,nz)
4 
5  implicit none
6  integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
7  real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk,z
8 
9  IF (init.ne.0) THEN
10  CALL rffti(n,table)
11  ELSE
12 !OCL NOVREC
13  DO j=1,m
14  y(1,j)=x(1,j)
15  DO i=2,n
16  y(i,j)=x(i+1,j)
17  ENDDO
18  CALL rfftb(n,y(1,j),table)
19  DO i=1,n
20  y(i,j)=scale*y(i,j)
21  ENDDO
22  ENDDO
23  ENDIF
24 
25  RETURN
26  END
27 
28  SUBROUTINE scrft(init,x,ldx,y,ldy,n,m,isign,scale,
29  & table,n1,wrk,n2,z,nz)
30 
31  implicit none
32  integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
33  real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk,z
34 
35  IF (init.ne.0) THEN
36  CALL rffti(n,table)
37  ELSE
38 !OCL NOVREC
39  DO j=1,m
40  y(1,j)=x(1,j)
41  DO i=2,n
42  y(i,j)=x(i+1,j)
43  ENDDO
44  CALL rfftb(n,y(1,j),table)
45  DO i=1,n
46  y(i,j)=scale*y(i,j)
47  ENDDO
48  ENDDO
49  ENDIF
50 
51  RETURN
52  END
53 c
54 c***********************************************************************
55 c
56  SUBROUTINE csfft(isign,n,scale,x,y,table,work,isys)
57 
58  implicit none
59  integer isign,n,isys,i
60  real scale,x(*),y(*),table(*),work(*)
61 
62  IF (isign.eq.0) THEN
63  CALL rffti(n,table)
64  ENDIF
65  IF (isign.eq.1) THEN
66  y(1)=x(1)
67  DO i=2,n
68  y(i)=x(i+1)
69  ENDDO
70  CALL rfftb(n,y,table)
71  DO i=1,n
72  y(i)=scale*y(i)
73  ENDDO
74  ENDIF
75 
76  RETURN
77  END
78 c
79 c***********************************************************************
80 c
81  SUBROUTINE drcft(init,x,ldx,y,ldy,n,m,isign,scale,
82  & table,n1,wrk,n2,z,nz)
83 
84  implicit none
85  integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
86  real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk,z
87 
88  IF (init.ne.0) THEN
89  CALL rffti(n,table)
90  ELSE
91  DO j=1,m
92  DO i=1,n
93  y(i,j)=x(i,j)
94  ENDDO
95  CALL rfftf(n,y(1,j),table)
96  DO i=1,n
97  y(i,j)=scale*y(i,j)
98  ENDDO
99  DO i=n,2,-1
100  y(i+1,j)=y(i,j)
101  ENDDO
102  y(2,j)=0.
103 C 01/17/2013 vvvvvvvvvvvvv E.Mirvis added ver 2.0.1 by S.Moorthi request. No +|- demo.
104  y(n+2,j) = 0.
105  ENDDO
106  ENDIF
107 
108  RETURN
109  END
110 
111  SUBROUTINE srcft(init,x,ldx,y,ldy,n,m,isign,scale,
112  & table,n1,wrk,n2,z,nz)
113 
114  implicit none
115  integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
116  real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk,z
117 
118  IF (init.ne.0) THEN
119  CALL rffti(n,table)
120  ELSE
121  DO j=1,m
122  DO i=1,n
123  y(i,j)=x(i,j)
124  ENDDO
125  CALL rfftf(n,y(1,j),table)
126  DO i=1,n
127  y(i,j)=scale*y(i,j)
128  ENDDO
129  DO i=n,2,-1
130  y(i+1,j)=y(i,j)
131  ENDDO
132  y(2,j)=0.
133  y(n+2,j) = 0.
134 C 01/17/2013 ^^^^^^^^^^E.Mirvis added ver 2.0.1 by S.Moorthi request. No +|- demo.
135  ENDDO
136  ENDIF
137 
138  RETURN
139  END
140 c
141 c***********************************************************************
142 c
143  SUBROUTINE scfft(isign,n,scale,x,y,table,work,isys)
144 
145  implicit none
146  integer isign,n,isys,i
147  real scale,x(*),y(*),table(*),work(*)
148 
149  IF (isign.eq.0) THEN
150  CALL rffti(n,table)
151  ENDIF
152  IF (isign.eq.-1) THEN
153  DO i=1,n
154  y(i)=x(i)
155  ENDDO
156  CALL rfftf(n,y,table)
157  DO i=1,n
158  y(i)=scale*y(i)
159  ENDDO
160  DO i=n,2,-1
161  y(i+1)=y(i)
162  ENDDO
163  y(2)=0.
164  ENDIF
165 
166  RETURN
167  END
168 c
169 c ******************************************************************
170 c ******************************************************************
171 c ****** ******
172 c ****** FFTPACK ******
173 c ****** ******
174 c ******************************************************************
175 c ******************************************************************
176 c
177  SUBROUTINE rfftf (N,R,WSAVE)
178  dimension r(1) ,wsave(1)
179  IF (n .EQ. 1) RETURN
180  CALL rfftf1 (n,r,wsave,wsave(n+1),wsave(2*n+1))
181  RETURN
182  END
183  SUBROUTINE rfftb (N,R,WSAVE)
184  dimension r(1) ,wsave(1)
185  IF (n .EQ. 1) RETURN
186  CALL rfftb1 (n,r,wsave,wsave(n+1),wsave(2*n+1))
187  RETURN
188  END
189  SUBROUTINE rffti (N,WSAVE)
190  dimension wsave(1)
191  IF (n .EQ. 1) RETURN
192  CALL rffti1 (n,wsave(n+1),wsave(2*n+1))
193  RETURN
194  END
195  SUBROUTINE rfftb1 (N,C,CH,WA,IFAC)
196  dimension ch(1) ,c(1) ,wa(1) ,ifac(*)
197  nf = ifac(2)
198  na = 0
199  l1 = 1
200  iw = 1
201  DO 116 k1=1,nf
202  ip = ifac(k1+2)
203  l2 = ip*l1
204  ido = n/l2
205  idl1 = ido*l1
206  IF (ip .NE. 4) GO TO 103
207  ix2 = iw+ido
208  ix3 = ix2+ido
209  IF (na .NE. 0) GO TO 101
210  CALL radb4 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3))
211  GO TO 102
212  101 CALL radb4 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3))
213  102 na = 1-na
214  GO TO 115
215  103 IF (ip .NE. 2) GO TO 106
216  IF (na .NE. 0) GO TO 104
217  CALL radb2 (ido,l1,c,ch,wa(iw))
218  GO TO 105
219  104 CALL radb2 (ido,l1,ch,c,wa(iw))
220  105 na = 1-na
221  GO TO 115
222  106 IF (ip .NE. 3) GO TO 109
223  ix2 = iw+ido
224  IF (na .NE. 0) GO TO 107
225  CALL radb3 (ido,l1,c,ch,wa(iw),wa(ix2))
226  GO TO 108
227  107 CALL radb3 (ido,l1,ch,c,wa(iw),wa(ix2))
228  108 na = 1-na
229  GO TO 115
230  109 IF (ip .NE. 5) GO TO 112
231  ix2 = iw+ido
232  ix3 = ix2+ido
233  ix4 = ix3+ido
234  IF (na .NE. 0) GO TO 110
235  CALL radb5 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))
236  GO TO 111
237  110 CALL radb5 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))
238  111 na = 1-na
239  GO TO 115
240  112 IF (na .NE. 0) GO TO 113
241  CALL radbg (ido,ip,l1,idl1,c,c,c,ch,ch,wa(iw))
242  GO TO 114
243  113 CALL radbg (ido,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))
244  114 IF (ido .EQ. 1) na = 1-na
245  115 l1 = l2
246  iw = iw+(ip-1)*ido
247  116 CONTINUE
248  IF (na .EQ. 0) RETURN
249  DO 117 i=1,n
250  c(i) = ch(i)
251  117 CONTINUE
252  RETURN
253  END
254 
255 
256  SUBROUTINE rfftf1 (N,C,CH,WA,IFAC)
257  dimension ch(1) ,c(1) ,wa(1) ,ifac(*)
258  nf = ifac(2)
259  na = 1
260  l2 = n
261  iw = n
262  DO 111 k1=1,nf
263  kh = nf-k1
264  ip = ifac(kh+3)
265  l1 = l2/ip
266  ido = n/l2
267  idl1 = ido*l1
268  iw = iw-(ip-1)*ido
269  na = 1-na
270  IF (ip .NE. 4) GO TO 102
271  ix2 = iw+ido
272  ix3 = ix2+ido
273  IF (na .NE. 0) GO TO 101
274  CALL radf4 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3))
275  GO TO 110
276  101 CALL radf4 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3))
277  GO TO 110
278  102 IF (ip .NE. 2) GO TO 104
279  IF (na .NE. 0) GO TO 103
280  CALL radf2 (ido,l1,c,ch,wa(iw))
281  GO TO 110
282  103 CALL radf2 (ido,l1,ch,c,wa(iw))
283  GO TO 110
284  104 IF (ip .NE. 3) GO TO 106
285  ix2 = iw+ido
286  IF (na .NE. 0) GO TO 105
287  CALL radf3 (ido,l1,c,ch,wa(iw),wa(ix2))
288  GO TO 110
289  105 CALL radf3 (ido,l1,ch,c,wa(iw),wa(ix2))
290  GO TO 110
291  106 IF (ip .NE. 5) GO TO 108
292  ix2 = iw+ido
293  ix3 = ix2+ido
294  ix4 = ix3+ido
295  IF (na .NE. 0) GO TO 107
296  CALL radf5 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))
297  GO TO 110
298  107 CALL radf5 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))
299  GO TO 110
300  108 IF (ido .EQ. 1) na = 1-na
301  IF (na .NE. 0) GO TO 109
302  CALL radfg (ido,ip,l1,idl1,c,c,c,ch,ch,wa(iw))
303  na = 1
304  GO TO 110
305  109 CALL radfg (ido,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))
306  na = 0
307  110 l2 = l1
308  111 CONTINUE
309  IF (na .EQ. 1) RETURN
310  DO 112 i=1,n
311  c(i) = ch(i)
312  112 CONTINUE
313  RETURN
314  END
315 
316 
317  SUBROUTINE rffti1 (N,WA,IFAC)
318  dimension wa(1) ,ifac(*) ,ntryh(4)
319  DATA ntryh(1),ntryh(2),ntryh(3),ntryh(4)/4,2,3,5/
320  nl = n
321  nf = 0
322  j = 0
323  101 j = j+1
324  IF (j-4) 102,102,103
325  102 ntry = ntryh(j)
326  GO TO 104
327  103 ntry = ntry+2
328  104 nq = nl/ntry
329  nr = nl-ntry*nq
330  IF (nr) 101,105,101
331  105 nf = nf+1
332  ifac(nf+2) = ntry
333  nl = nq
334  IF (ntry .NE. 2) GO TO 107
335  IF (nf .EQ. 1) GO TO 107
336  DO 106 i=2,nf
337  ib = nf-i+2
338  ifac(ib+2) = ifac(ib+1)
339  106 CONTINUE
340  ifac(3) = 2
341  107 IF (nl .NE. 1) GO TO 104
342  ifac(1) = n
343  ifac(2) = nf
344  tpi = 6.28318530717959
345  argh = tpi/float(n)
346  is = 0
347  nfm1 = nf-1
348  l1 = 1
349  IF (nfm1 .EQ. 0) RETURN
350 !OCL NOVREC
351  DO 110 k1=1,nfm1
352  ip = ifac(k1+2)
353  ld = 0
354  l2 = l1*ip
355  ido = n/l2
356  ipm = ip-1
357  DO 109 j=1,ipm
358  ld = ld+l1
359  i = is
360  argld = float(ld)*argh
361  fi = 0
362 !OCL SCALAR
363  DO 108 ii=3,ido,2
364  i = i+2
365  fi = fi+1
366  arg = fi*argld
367  wa(i-1) = cos(arg)
368  wa(i) = sin(arg)
369  108 CONTINUE
370  is = is+ido
371  109 CONTINUE
372  l1 = l2
373  110 CONTINUE
374  RETURN
375  END
376 
377 
378  SUBROUTINE radb2 (IDO,L1,CC,CH,WA1)
379  dimension cc(ido,2,l1) ,ch(ido,l1,2) ,
380  1 wa1(1)
381  DO 101 k=1,l1
382  ch(1,k,1) = cc(1,1,k)+cc(ido,2,k)
383  ch(1,k,2) = cc(1,1,k)-cc(ido,2,k)
384  101 CONTINUE
385  IF (ido-2) 107,105,102
386  102 idp2 = ido+2
387 !OCL NOVREC
388  DO 104 k=1,l1
389  DO 103 i=3,ido,2
390  ic = idp2-i
391  ch(i-1,k,1) = cc(i-1,1,k)+cc(ic-1,2,k)
392  tr2 = cc(i-1,1,k)-cc(ic-1,2,k)
393  ch(i,k,1) = cc(i,1,k)-cc(ic,2,k)
394  ti2 = cc(i,1,k)+cc(ic,2,k)
395  ch(i-1,k,2) = wa1(i-2)*tr2-wa1(i-1)*ti2
396  ch(i,k,2) = wa1(i-2)*ti2+wa1(i-1)*tr2
397  103 CONTINUE
398  104 CONTINUE
399  IF (mod(ido,2) .EQ. 1) RETURN
400  105 DO 106 k=1,l1
401  ch(ido,k,1) = cc(ido,1,k)+cc(ido,1,k)
402  ch(ido,k,2) = -(cc(1,2,k)+cc(1,2,k))
403  106 CONTINUE
404  107 RETURN
405  END
406 
407 
408  SUBROUTINE radb3 (IDO,L1,CC,CH,WA1,WA2)
409  dimension cc(ido,3,l1) ,ch(ido,l1,3) ,
410  1 wa1(1) ,wa2(1)
411  DATA taur,taui /-.5,.866025403784439/
412  DO 101 k=1,l1
413  tr2 = cc(ido,2,k)+cc(ido,2,k)
414  cr2 = cc(1,1,k)+taur*tr2
415  ch(1,k,1) = cc(1,1,k)+tr2
416  ci3 = taui*(cc(1,3,k)+cc(1,3,k))
417  ch(1,k,2) = cr2-ci3
418  ch(1,k,3) = cr2+ci3
419  101 CONTINUE
420  IF (ido .EQ. 1) RETURN
421  idp2 = ido+2
422 !OCL NOVREC
423  DO 103 k=1,l1
424  DO 102 i=3,ido,2
425  ic = idp2-i
426  tr2 = cc(i-1,3,k)+cc(ic-1,2,k)
427  cr2 = cc(i-1,1,k)+taur*tr2
428  ch(i-1,k,1) = cc(i-1,1,k)+tr2
429  ti2 = cc(i,3,k)-cc(ic,2,k)
430  ci2 = cc(i,1,k)+taur*ti2
431  ch(i,k,1) = cc(i,1,k)+ti2
432  cr3 = taui*(cc(i-1,3,k)-cc(ic-1,2,k))
433  ci3 = taui*(cc(i,3,k)+cc(ic,2,k))
434  dr2 = cr2-ci3
435  dr3 = cr2+ci3
436  di2 = ci2+cr3
437  di3 = ci2-cr3
438  ch(i-1,k,2) = wa1(i-2)*dr2-wa1(i-1)*di2
439  ch(i,k,2) = wa1(i-2)*di2+wa1(i-1)*dr2
440  ch(i-1,k,3) = wa2(i-2)*dr3-wa2(i-1)*di3
441  ch(i,k,3) = wa2(i-2)*di3+wa2(i-1)*dr3
442  102 CONTINUE
443  103 CONTINUE
444  RETURN
445  END
446 
447 
448  SUBROUTINE radb4 (IDO,L1,CC,CH,WA1,WA2,WA3)
449  dimension cc(ido,4,l1) ,ch(ido,l1,4) ,
450  1 wa1(1) ,wa2(1) ,wa3(1)
451  DATA sqrt2 /1.414213562373095/
452  DO 101 k=1,l1
453  tr1 = cc(1,1,k)-cc(ido,4,k)
454  tr2 = cc(1,1,k)+cc(ido,4,k)
455  tr3 = cc(ido,2,k)+cc(ido,2,k)
456  tr4 = cc(1,3,k)+cc(1,3,k)
457  ch(1,k,1) = tr2+tr3
458  ch(1,k,2) = tr1-tr4
459  ch(1,k,3) = tr2-tr3
460  ch(1,k,4) = tr1+tr4
461  101 CONTINUE
462  IF (ido-2) 107,105,102
463  102 idp2 = ido+2
464 !OCL NOVREC
465  DO 104 k=1,l1
466  DO 103 i=3,ido,2
467  ic = idp2-i
468  ti1 = cc(i,1,k)+cc(ic,4,k)
469  ti2 = cc(i,1,k)-cc(ic,4,k)
470  ti3 = cc(i,3,k)-cc(ic,2,k)
471  tr4 = cc(i,3,k)+cc(ic,2,k)
472  tr1 = cc(i-1,1,k)-cc(ic-1,4,k)
473  tr2 = cc(i-1,1,k)+cc(ic-1,4,k)
474  ti4 = cc(i-1,3,k)-cc(ic-1,2,k)
475  tr3 = cc(i-1,3,k)+cc(ic-1,2,k)
476  ch(i-1,k,1) = tr2+tr3
477  cr3 = tr2-tr3
478  ch(i,k,1) = ti2+ti3
479  ci3 = ti2-ti3
480  cr2 = tr1-tr4
481  cr4 = tr1+tr4
482  ci2 = ti1+ti4
483  ci4 = ti1-ti4
484  ch(i-1,k,2) = wa1(i-2)*cr2-wa1(i-1)*ci2
485  ch(i,k,2) = wa1(i-2)*ci2+wa1(i-1)*cr2
486  ch(i-1,k,3) = wa2(i-2)*cr3-wa2(i-1)*ci3
487  ch(i,k,3) = wa2(i-2)*ci3+wa2(i-1)*cr3
488  ch(i-1,k,4) = wa3(i-2)*cr4-wa3(i-1)*ci4
489  ch(i,k,4) = wa3(i-2)*ci4+wa3(i-1)*cr4
490  103 CONTINUE
491  104 CONTINUE
492  IF (mod(ido,2) .EQ. 1) RETURN
493  105 CONTINUE
494  DO 106 k=1,l1
495  ti1 = cc(1,2,k)+cc(1,4,k)
496  ti2 = cc(1,4,k)-cc(1,2,k)
497  tr1 = cc(ido,1,k)-cc(ido,3,k)
498  tr2 = cc(ido,1,k)+cc(ido,3,k)
499  ch(ido,k,1) = tr2+tr2
500  ch(ido,k,2) = sqrt2*(tr1-ti1)
501  ch(ido,k,3) = ti2+ti2
502  ch(ido,k,4) = -sqrt2*(tr1+ti1)
503  106 CONTINUE
504  107 RETURN
505  END
506 
507 
508  SUBROUTINE radb5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
509  dimension cc(ido,5,l1) ,ch(ido,l1,5) ,
510  1 wa1(1) ,wa2(1) ,wa3(1) ,wa4(1)
511  DATA tr11,ti11,tr12,ti12 /.309016994374947,.951056516295154,
512  1-.809016994374947,.587785252292473/
513  DO 101 k=1,l1
514  ti5 = cc(1,3,k)+cc(1,3,k)
515  ti4 = cc(1,5,k)+cc(1,5,k)
516  tr2 = cc(ido,2,k)+cc(ido,2,k)
517  tr3 = cc(ido,4,k)+cc(ido,4,k)
518  ch(1,k,1) = cc(1,1,k)+tr2+tr3
519  cr2 = cc(1,1,k)+tr11*tr2+tr12*tr3
520  cr3 = cc(1,1,k)+tr12*tr2+tr11*tr3
521  ci5 = ti11*ti5+ti12*ti4
522  ci4 = ti12*ti5-ti11*ti4
523  ch(1,k,2) = cr2-ci5
524  ch(1,k,3) = cr3-ci4
525  ch(1,k,4) = cr3+ci4
526  ch(1,k,5) = cr2+ci5
527  101 CONTINUE
528  IF (ido .EQ. 1) RETURN
529  idp2 = ido+2
530  DO 103 k=1,l1
531  DO 102 i=3,ido,2
532  ic = idp2-i
533  ti5 = cc(i,3,k)+cc(ic,2,k)
534  ti2 = cc(i,3,k)-cc(ic,2,k)
535  ti4 = cc(i,5,k)+cc(ic,4,k)
536  ti3 = cc(i,5,k)-cc(ic,4,k)
537  tr5 = cc(i-1,3,k)-cc(ic-1,2,k)
538  tr2 = cc(i-1,3,k)+cc(ic-1,2,k)
539  tr4 = cc(i-1,5,k)-cc(ic-1,4,k)
540  tr3 = cc(i-1,5,k)+cc(ic-1,4,k)
541  ch(i-1,k,1) = cc(i-1,1,k)+tr2+tr3
542  ch(i,k,1) = cc(i,1,k)+ti2+ti3
543  cr2 = cc(i-1,1,k)+tr11*tr2+tr12*tr3
544  ci2 = cc(i,1,k)+tr11*ti2+tr12*ti3
545  cr3 = cc(i-1,1,k)+tr12*tr2+tr11*tr3
546  ci3 = cc(i,1,k)+tr12*ti2+tr11*ti3
547  cr5 = ti11*tr5+ti12*tr4
548  ci5 = ti11*ti5+ti12*ti4
549  cr4 = ti12*tr5-ti11*tr4
550  ci4 = ti12*ti5-ti11*ti4
551  dr3 = cr3-ci4
552  dr4 = cr3+ci4
553  di3 = ci3+cr4
554  di4 = ci3-cr4
555  dr5 = cr2+ci5
556  dr2 = cr2-ci5
557  di5 = ci2-cr5
558  di2 = ci2+cr5
559  ch(i-1,k,2) = wa1(i-2)*dr2-wa1(i-1)*di2
560  ch(i,k,2) = wa1(i-2)*di2+wa1(i-1)*dr2
561  ch(i-1,k,3) = wa2(i-2)*dr3-wa2(i-1)*di3
562  ch(i,k,3) = wa2(i-2)*di3+wa2(i-1)*dr3
563  ch(i-1,k,4) = wa3(i-2)*dr4-wa3(i-1)*di4
564  ch(i,k,4) = wa3(i-2)*di4+wa3(i-1)*dr4
565  ch(i-1,k,5) = wa4(i-2)*dr5-wa4(i-1)*di5
566  ch(i,k,5) = wa4(i-2)*di5+wa4(i-1)*dr5
567  102 CONTINUE
568  103 CONTINUE
569  RETURN
570  END
571 
572 
573  SUBROUTINE radbg (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
574  dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
575  1 c1(ido,l1,ip) ,c2(idl1,ip),
576  2 ch2(idl1,ip) ,wa(1)
577  DATA tpi/6.28318530717959/
578  arg = tpi/float(ip)
579  dcp = cos(arg)
580  dsp = sin(arg)
581  idp2 = ido+2
582  nbd = (ido-1)/2
583  ipp2 = ip+2
584  ipph = (ip+1)/2
585  IF (ido .LT. l1) GO TO 103
586  DO 102 k=1,l1
587  DO 101 i=1,ido
588  ch(i,k,1) = cc(i,1,k)
589  101 CONTINUE
590  102 CONTINUE
591  GO TO 106
592  103 DO 105 i=1,ido
593  DO 104 k=1,l1
594  ch(i,k,1) = cc(i,1,k)
595  104 CONTINUE
596  105 CONTINUE
597 !OCL NOVREC
598  106 DO 108 j=2,ipph
599  jc = ipp2-j
600  j2 = j+j
601  DO 107 k=1,l1
602  ch(1,k,j) = cc(ido,j2-2,k)+cc(ido,j2-2,k)
603  ch(1,k,jc) = cc(1,j2-1,k)+cc(1,j2-1,k)
604  107 CONTINUE
605  108 CONTINUE
606  IF (ido .EQ. 1) GO TO 116
607  IF (nbd .LT. l1) GO TO 112
608 !OCL NOVREC
609  DO 111 j=2,ipph
610  jc = ipp2-j
611  DO 110 k=1,l1
612  DO 109 i=3,ido,2
613  ic = idp2-i
614  ch(i-1,k,j) = cc(i-1,2*j-1,k)+cc(ic-1,2*j-2,k)
615  ch(i-1,k,jc) = cc(i-1,2*j-1,k)-cc(ic-1,2*j-2,k)
616  ch(i,k,j) = cc(i,2*j-1,k)-cc(ic,2*j-2,k)
617  ch(i,k,jc) = cc(i,2*j-1,k)+cc(ic,2*j-2,k)
618  109 CONTINUE
619  110 CONTINUE
620  111 CONTINUE
621  GO TO 116
622  112 DO 115 j=2,ipph
623  jc = ipp2-j
624  DO 114 i=3,ido,2
625  ic = idp2-i
626  DO 113 k=1,l1
627  ch(i-1,k,j) = cc(i-1,2*j-1,k)+cc(ic-1,2*j-2,k)
628  ch(i-1,k,jc) = cc(i-1,2*j-1,k)-cc(ic-1,2*j-2,k)
629  ch(i,k,j) = cc(i,2*j-1,k)-cc(ic,2*j-2,k)
630  ch(i,k,jc) = cc(i,2*j-1,k)+cc(ic,2*j-2,k)
631  113 CONTINUE
632  114 CONTINUE
633  115 CONTINUE
634  116 ar1 = 1.
635  ai1 = 0.
636 !OCL NOVREC
637  DO 120 l=2,ipph
638  lc = ipp2-l
639  ar1h = dcp*ar1-dsp*ai1
640  ai1 = dcp*ai1+dsp*ar1
641  ar1 = ar1h
642  DO 117 ik=1,idl1
643  c2(ik,l) = ch2(ik,1)+ar1*ch2(ik,2)
644  c2(ik,lc) = ai1*ch2(ik,ip)
645  117 CONTINUE
646  dc2 = ar1
647  ds2 = ai1
648  ar2 = ar1
649  ai2 = ai1
650 !OCL NOVREC
651  DO 119 j=3,ipph
652  jc = ipp2-j
653  ar2h = dc2*ar2-ds2*ai2
654  ai2 = dc2*ai2+ds2*ar2
655  ar2 = ar2h
656  DO 118 ik=1,idl1
657  c2(ik,l) = c2(ik,l)+ar2*ch2(ik,j)
658  c2(ik,lc) = c2(ik,lc)+ai2*ch2(ik,jc)
659  118 CONTINUE
660  119 CONTINUE
661  120 CONTINUE
662 !OCL NOVREC
663  DO 122 j=2,ipph
664  DO 121 ik=1,idl1
665  ch2(ik,1) = ch2(ik,1)+ch2(ik,j)
666  121 CONTINUE
667  122 CONTINUE
668 !OCL NOVREC
669  DO 124 j=2,ipph
670  jc = ipp2-j
671  DO 123 k=1,l1
672  ch(1,k,j) = c1(1,k,j)-c1(1,k,jc)
673  ch(1,k,jc) = c1(1,k,j)+c1(1,k,jc)
674  123 CONTINUE
675  124 CONTINUE
676  IF (ido .EQ. 1) GO TO 132
677  IF (nbd .LT. l1) GO TO 128
678 !OCL NOVREC
679  DO 127 j=2,ipph
680  jc = ipp2-j
681  DO 126 k=1,l1
682  DO 125 i=3,ido,2
683  ch(i-1,k,j) = c1(i-1,k,j)-c1(i,k,jc)
684  ch(i-1,k,jc) = c1(i-1,k,j)+c1(i,k,jc)
685  ch(i,k,j) = c1(i,k,j)+c1(i-1,k,jc)
686  ch(i,k,jc) = c1(i,k,j)-c1(i-1,k,jc)
687  125 CONTINUE
688  126 CONTINUE
689  127 CONTINUE
690  GO TO 132
691  128 DO 131 j=2,ipph
692  jc = ipp2-j
693  DO 130 i=3,ido,2
694  DO 129 k=1,l1
695  ch(i-1,k,j) = c1(i-1,k,j)-c1(i,k,jc)
696  ch(i-1,k,jc) = c1(i-1,k,j)+c1(i,k,jc)
697  ch(i,k,j) = c1(i,k,j)+c1(i-1,k,jc)
698  ch(i,k,jc) = c1(i,k,j)-c1(i-1,k,jc)
699  129 CONTINUE
700  130 CONTINUE
701  131 CONTINUE
702  132 CONTINUE
703  IF (ido .EQ. 1) RETURN
704  DO 133 ik=1,idl1
705  c2(ik,1) = ch2(ik,1)
706  133 CONTINUE
707  DO 135 j=2,ip
708  DO 134 k=1,l1
709  c1(1,k,j) = ch(1,k,j)
710  134 CONTINUE
711  135 CONTINUE
712  IF (nbd .GT. l1) GO TO 139
713  is = -ido
714  DO 138 j=2,ip
715  is = is+ido
716  idij = is
717  DO 137 i=3,ido,2
718  idij = idij+2
719  DO 136 k=1,l1
720  c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
721  c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
722  136 CONTINUE
723  137 CONTINUE
724  138 CONTINUE
725  GO TO 143
726  139 is = -ido
727 !OCL NOVREC
728  DO 142 j=2,ip
729  is = is+ido
730  DO 141 k=1,l1
731  idij = is
732  DO 140 i=3,ido,2
733  idij = idij+2
734  c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
735  c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
736  140 CONTINUE
737  141 CONTINUE
738  142 CONTINUE
739  143 RETURN
740  END
741 
742 
743  SUBROUTINE radf2 (IDO,L1,CC,CH,WA1)
744  dimension ch(ido,2,l1) ,cc(ido,l1,2) ,
745  1 wa1(1)
746  DO 101 k=1,l1
747  ch(1,1,k) = cc(1,k,1)+cc(1,k,2)
748  ch(ido,2,k) = cc(1,k,1)-cc(1,k,2)
749  101 CONTINUE
750  IF (ido-2) 107,105,102
751  102 idp2 = ido+2
752  DO 104 k=1,l1
753  DO 103 i=3,ido,2
754  ic = idp2-i
755  tr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
756  ti2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
757  ch(i,1,k) = cc(i,k,1)+ti2
758  ch(ic,2,k) = ti2-cc(i,k,1)
759  ch(i-1,1,k) = cc(i-1,k,1)+tr2
760  ch(ic-1,2,k) = cc(i-1,k,1)-tr2
761  103 CONTINUE
762  104 CONTINUE
763  IF (mod(ido,2) .EQ. 1) RETURN
764  105 DO 106 k=1,l1
765  ch(1,2,k) = -cc(ido,k,2)
766  ch(ido,1,k) = cc(ido,k,1)
767  106 CONTINUE
768  107 RETURN
769  END
770 
771 
772  SUBROUTINE radf3 (IDO,L1,CC,CH,WA1,WA2)
773  dimension ch(ido,3,l1) ,cc(ido,l1,3) ,
774  1 wa1(1) ,wa2(1)
775  DATA taur,taui /-.5,.866025403784439/
776  DO 101 k=1,l1
777  cr2 = cc(1,k,2)+cc(1,k,3)
778  ch(1,1,k) = cc(1,k,1)+cr2
779  ch(1,3,k) = taui*(cc(1,k,3)-cc(1,k,2))
780  ch(ido,2,k) = cc(1,k,1)+taur*cr2
781  101 CONTINUE
782  IF (ido .EQ. 1) RETURN
783  idp2 = ido+2
784  DO 103 k=1,l1
785  DO 102 i=3,ido,2
786  ic = idp2-i
787  dr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
788  di2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
789  dr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
790  di3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
791  cr2 = dr2+dr3
792  ci2 = di2+di3
793  ch(i-1,1,k) = cc(i-1,k,1)+cr2
794  ch(i,1,k) = cc(i,k,1)+ci2
795  tr2 = cc(i-1,k,1)+taur*cr2
796  ti2 = cc(i,k,1)+taur*ci2
797  tr3 = taui*(di2-di3)
798  ti3 = taui*(dr3-dr2)
799  ch(i-1,3,k) = tr2+tr3
800  ch(ic-1,2,k) = tr2-tr3
801  ch(i,3,k) = ti2+ti3
802  ch(ic,2,k) = ti3-ti2
803  102 CONTINUE
804  103 CONTINUE
805  RETURN
806  END
807  SUBROUTINE radf4 (IDO,L1,CC,CH,WA1,WA2,WA3)
808  dimension cc(ido,l1,4) ,ch(ido,4,l1) ,
809  1 wa1(1) ,wa2(1) ,wa3(1)
810  DATA hsqt2 /.7071067811865475/
811  DO 101 k=1,l1
812  tr1 = cc(1,k,2)+cc(1,k,4)
813  tr2 = cc(1,k,1)+cc(1,k,3)
814  ch(1,1,k) = tr1+tr2
815  ch(ido,4,k) = tr2-tr1
816  ch(ido,2,k) = cc(1,k,1)-cc(1,k,3)
817  ch(1,3,k) = cc(1,k,4)-cc(1,k,2)
818  101 CONTINUE
819  IF (ido-2) 107,105,102
820  102 idp2 = ido+2
821 !OCL NOVREC
822  DO 104 k=1,l1
823  DO 103 i=3,ido,2
824  ic = idp2-i
825  cr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
826  ci2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
827  cr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
828  ci3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
829  cr4 = wa3(i-2)*cc(i-1,k,4)+wa3(i-1)*cc(i,k,4)
830  ci4 = wa3(i-2)*cc(i,k,4)-wa3(i-1)*cc(i-1,k,4)
831  tr1 = cr2+cr4
832  tr4 = cr4-cr2
833  ti1 = ci2+ci4
834  ti4 = ci2-ci4
835  ti2 = cc(i,k,1)+ci3
836  ti3 = cc(i,k,1)-ci3
837  tr2 = cc(i-1,k,1)+cr3
838  tr3 = cc(i-1,k,1)-cr3
839  ch(i-1,1,k) = tr1+tr2
840  ch(ic-1,4,k) = tr2-tr1
841  ch(i,1,k) = ti1+ti2
842  ch(ic,4,k) = ti1-ti2
843  ch(i-1,3,k) = ti4+tr3
844  ch(ic-1,2,k) = tr3-ti4
845  ch(i,3,k) = tr4+ti3
846  ch(ic,2,k) = tr4-ti3
847  103 CONTINUE
848  104 CONTINUE
849  IF (mod(ido,2) .EQ. 1) RETURN
850  105 CONTINUE
851  DO 106 k=1,l1
852  ti1 = -hsqt2*(cc(ido,k,2)+cc(ido,k,4))
853  tr1 = hsqt2*(cc(ido,k,2)-cc(ido,k,4))
854  ch(ido,1,k) = tr1+cc(ido,k,1)
855  ch(ido,3,k) = cc(ido,k,1)-tr1
856  ch(1,2,k) = ti1-cc(ido,k,3)
857  ch(1,4,k) = ti1+cc(ido,k,3)
858  106 CONTINUE
859  107 RETURN
860  END
861 
862 
863  SUBROUTINE radf5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
864  dimension cc(ido,l1,5) ,ch(ido,5,l1) ,
865  1 wa1(1) ,wa2(1) ,wa3(1) ,wa4(1)
866  DATA tr11,ti11,tr12,ti12 /.309016994374947,.951056516295154,
867  1-.809016994374947,.587785252292473/
868  DO 101 k=1,l1
869  cr2 = cc(1,k,5)+cc(1,k,2)
870  ci5 = cc(1,k,5)-cc(1,k,2)
871  cr3 = cc(1,k,4)+cc(1,k,3)
872  ci4 = cc(1,k,4)-cc(1,k,3)
873  ch(1,1,k) = cc(1,k,1)+cr2+cr3
874  ch(ido,2,k) = cc(1,k,1)+tr11*cr2+tr12*cr3
875  ch(1,3,k) = ti11*ci5+ti12*ci4
876  ch(ido,4,k) = cc(1,k,1)+tr12*cr2+tr11*cr3
877  ch(1,5,k) = ti12*ci5-ti11*ci4
878  101 CONTINUE
879  IF (ido .EQ. 1) RETURN
880  idp2 = ido+2
881  DO 103 k=1,l1
882  DO 102 i=3,ido,2
883  ic = idp2-i
884  dr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
885  di2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
886  dr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
887  di3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
888  dr4 = wa3(i-2)*cc(i-1,k,4)+wa3(i-1)*cc(i,k,4)
889  di4 = wa3(i-2)*cc(i,k,4)-wa3(i-1)*cc(i-1,k,4)
890  dr5 = wa4(i-2)*cc(i-1,k,5)+wa4(i-1)*cc(i,k,5)
891  di5 = wa4(i-2)*cc(i,k,5)-wa4(i-1)*cc(i-1,k,5)
892  cr2 = dr2+dr5
893  ci5 = dr5-dr2
894  cr5 = di2-di5
895  ci2 = di2+di5
896  cr3 = dr3+dr4
897  ci4 = dr4-dr3
898  cr4 = di3-di4
899  ci3 = di3+di4
900  ch(i-1,1,k) = cc(i-1,k,1)+cr2+cr3
901  ch(i,1,k) = cc(i,k,1)+ci2+ci3
902  tr2 = cc(i-1,k,1)+tr11*cr2+tr12*cr3
903  ti2 = cc(i,k,1)+tr11*ci2+tr12*ci3
904  tr3 = cc(i-1,k,1)+tr12*cr2+tr11*cr3
905  ti3 = cc(i,k,1)+tr12*ci2+tr11*ci3
906  tr5 = ti11*cr5+ti12*cr4
907  ti5 = ti11*ci5+ti12*ci4
908  tr4 = ti12*cr5-ti11*cr4
909  ti4 = ti12*ci5-ti11*ci4
910  ch(i-1,3,k) = tr2+tr5
911  ch(ic-1,2,k) = tr2-tr5
912  ch(i,3,k) = ti2+ti5
913  ch(ic,2,k) = ti5-ti2
914  ch(i-1,5,k) = tr3+tr4
915  ch(ic-1,4,k) = tr3-tr4
916  ch(i,5,k) = ti3+ti4
917  ch(ic,4,k) = ti4-ti3
918  102 CONTINUE
919  103 CONTINUE
920  RETURN
921  END
922 
923 
924  SUBROUTINE radfg (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
925  dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
926  1 c1(ido,l1,ip) ,c2(idl1,ip),
927  2 ch2(idl1,ip) ,wa(1)
928  DATA tpi/6.28318530717959/
929  arg = tpi/float(ip)
930  dcp = cos(arg)
931  dsp = sin(arg)
932  ipph = (ip+1)/2
933  ipp2 = ip+2
934  idp2 = ido+2
935  nbd = (ido-1)/2
936  IF (ido .EQ. 1) GO TO 119
937  DO 101 ik=1,idl1
938  ch2(ik,1) = c2(ik,1)
939  101 CONTINUE
940  DO 103 j=2,ip
941  DO 102 k=1,l1
942  ch(1,k,j) = c1(1,k,j)
943  102 CONTINUE
944  103 CONTINUE
945  IF (nbd .GT. l1) GO TO 107
946  is = -ido
947  DO 106 j=2,ip
948  is = is+ido
949  idij = is
950  DO 105 i=3,ido,2
951  idij = idij+2
952  DO 104 k=1,l1
953  ch(i-1,k,j) = wa(idij-1)*c1(i-1,k,j)+wa(idij)*c1(i,k,j)
954  ch(i,k,j) = wa(idij-1)*c1(i,k,j)-wa(idij)*c1(i-1,k,j)
955  104 CONTINUE
956  105 CONTINUE
957  106 CONTINUE
958  GO TO 111
959  107 is = -ido
960  DO 110 j=2,ip
961  is = is+ido
962  DO 109 k=1,l1
963  idij = is
964  DO 108 i=3,ido,2
965  idij = idij+2
966  ch(i-1,k,j) = wa(idij-1)*c1(i-1,k,j)+wa(idij)*c1(i,k,j)
967  ch(i,k,j) = wa(idij-1)*c1(i,k,j)-wa(idij)*c1(i-1,k,j)
968  108 CONTINUE
969  109 CONTINUE
970  110 CONTINUE
971  111 IF (nbd .LT. l1) GO TO 115
972  DO 114 j=2,ipph
973  jc = ipp2-j
974  DO 113 k=1,l1
975  DO 112 i=3,ido,2
976  c1(i-1,k,j) = ch(i-1,k,j)+ch(i-1,k,jc)
977  c1(i-1,k,jc) = ch(i,k,j)-ch(i,k,jc)
978  c1(i,k,j) = ch(i,k,j)+ch(i,k,jc)
979  c1(i,k,jc) = ch(i-1,k,jc)-ch(i-1,k,j)
980  112 CONTINUE
981  113 CONTINUE
982  114 CONTINUE
983  GO TO 121
984  115 DO 118 j=2,ipph
985  jc = ipp2-j
986  DO 117 i=3,ido,2
987  DO 116 k=1,l1
988  c1(i-1,k,j) = ch(i-1,k,j)+ch(i-1,k,jc)
989  c1(i-1,k,jc) = ch(i,k,j)-ch(i,k,jc)
990  c1(i,k,j) = ch(i,k,j)+ch(i,k,jc)
991  c1(i,k,jc) = ch(i-1,k,jc)-ch(i-1,k,j)
992  116 CONTINUE
993  117 CONTINUE
994  118 CONTINUE
995  GO TO 121
996  119 DO 120 ik=1,idl1
997  c2(ik,1) = ch2(ik,1)
998  120 CONTINUE
999  121 DO 123 j=2,ipph
1000  jc = ipp2-j
1001  DO 122 k=1,l1
1002  c1(1,k,j) = ch(1,k,j)+ch(1,k,jc)
1003  c1(1,k,jc) = ch(1,k,jc)-ch(1,k,j)
1004  122 CONTINUE
1005  123 CONTINUE
1006 C
1007  ar1 = 1.
1008  ai1 = 0.
1009  DO 127 l=2,ipph
1010  lc = ipp2-l
1011  ar1h = dcp*ar1-dsp*ai1
1012  ai1 = dcp*ai1+dsp*ar1
1013  ar1 = ar1h
1014  DO 124 ik=1,idl1
1015  ch2(ik,l) = c2(ik,1)+ar1*c2(ik,2)
1016  ch2(ik,lc) = ai1*c2(ik,ip)
1017  124 CONTINUE
1018  dc2 = ar1
1019  ds2 = ai1
1020  ar2 = ar1
1021  ai2 = ai1
1022  DO 126 j=3,ipph
1023  jc = ipp2-j
1024  ar2h = dc2*ar2-ds2*ai2
1025  ai2 = dc2*ai2+ds2*ar2
1026  ar2 = ar2h
1027  DO 125 ik=1,idl1
1028  ch2(ik,l) = ch2(ik,l)+ar2*c2(ik,j)
1029  ch2(ik,lc) = ch2(ik,lc)+ai2*c2(ik,jc)
1030  125 CONTINUE
1031  126 CONTINUE
1032  127 CONTINUE
1033  DO 129 j=2,ipph
1034  DO 128 ik=1,idl1
1035  ch2(ik,1) = ch2(ik,1)+c2(ik,j)
1036  128 CONTINUE
1037  129 CONTINUE
1038 C
1039  IF (ido .LT. l1) GO TO 132
1040  DO 131 k=1,l1
1041  DO 130 i=1,ido
1042  cc(i,1,k) = ch(i,k,1)
1043  130 CONTINUE
1044  131 CONTINUE
1045  GO TO 135
1046  132 DO 134 i=1,ido
1047  DO 133 k=1,l1
1048  cc(i,1,k) = ch(i,k,1)
1049  133 CONTINUE
1050  134 CONTINUE
1051  135 DO 137 j=2,ipph
1052  jc = ipp2-j
1053  j2 = j+j
1054  DO 136 k=1,l1
1055  cc(ido,j2-2,k) = ch(1,k,j)
1056  cc(1,j2-1,k) = ch(1,k,jc)
1057  136 CONTINUE
1058  137 CONTINUE
1059  IF (ido .EQ. 1) RETURN
1060  IF (nbd .LT. l1) GO TO 141
1061  DO 140 j=2,ipph
1062  jc = ipp2-j
1063  j2 = j+j
1064  DO 139 k=1,l1
1065  DO 138 i=3,ido,2
1066  ic = idp2-i
1067  cc(i-1,j2-1,k) = ch(i-1,k,j)+ch(i-1,k,jc)
1068  cc(ic-1,j2-2,k) = ch(i-1,k,j)-ch(i-1,k,jc)
1069  cc(i,j2-1,k) = ch(i,k,j)+ch(i,k,jc)
1070  cc(ic,j2-2,k) = ch(i,k,jc)-ch(i,k,j)
1071  138 CONTINUE
1072  139 CONTINUE
1073  140 CONTINUE
1074  RETURN
1075  141 DO 144 j=2,ipph
1076  jc = ipp2-j
1077  j2 = j+j
1078  DO 143 i=3,ido,2
1079  ic = idp2-i
1080  DO 142 k=1,l1
1081  cc(i-1,j2-1,k) = ch(i-1,k,j)+ch(i-1,k,jc)
1082  cc(ic-1,j2-2,k) = ch(i-1,k,j)-ch(i-1,k,jc)
1083  cc(i,j2-1,k) = ch(i,k,j)+ch(i,k,jc)
1084  cc(ic,j2-2,k) = ch(i,k,jc)-ch(i,k,j)
1085  142 CONTINUE
1086  143 CONTINUE
1087  144 CONTINUE
1088  RETURN
1089  END
1090 
1091 #endif