Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/SOURCE/temperton_fft.f90 @ 99

Last change on this file since 99 was 4, checked in by raasch, 18 years ago

Id keyword set as property for all *.f90 files
  • Property svn:keywords set to Id
File size: 51.2 KB
Line 
1 MODULE temperton_fft
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: temperton_fft.f90 4 2007-02-13 11:33:16Z raasch $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.2  2003/04/16 12:49:25  raasch
14! Abort in case of illegal factors
15!
16! Revision 1.1  2003/03/12 16:41:59  raasch
17! Initial revision
18!
19!
20! Description:
21! ------------
22! Fast Fourier transformation developed by Clive Temperton, ECMWF.
23!------------------------------------------------------------------------------!
24
25    IMPLICIT NONE
26
27    PRIVATE
28
29    PUBLIC set99, fft991cy
30
31
32    INTEGER           ::  nfax(10)   !  array used by *fft991*.
33    REAL, ALLOCATABLE ::  trig(:)    !  array used by *fft991*.
34
35!
36!-- nfft: maximum length of calls to *fft.
37#if defined( __nec )
38    INTEGER, PARAMETER ::  nfft = 256
39#else
40    INTEGER, PARAMETER ::  nfft =  32
41#endif
42
43    INTEGER, PARAMETER ::  nout =   6  ! standard output stream
44
45CONTAINS
46
47  SUBROUTINE fft991cy(a,work,trigs,ifax,inc,jump,n,lot,isign)
48
49    ! Description:
50    !
51    ! Calls fortran-versions of fft's.
52    !
53    ! Method:
54    !
55    ! Subroutine 'fft991cy' - multiple fast real periodic transform
56    ! supersedes previous routine 'fft991cy'.
57    !
58    ! Real transform of length n performed by removing redundant
59    ! operations from complex transform of length n.
60    !
61    ! a       is the array containing input & output data.
62    ! work    is an area of size (n+1)*min(lot,nfft).
63    ! trigs   is a previously prepared list of trig function values.
64    ! ifax    is a previously prepared list of factors of n.
65    ! inc     is the increment within each data 'vector'
66    !         (e.g. inc=1 for consecutively stored data).
67    ! jump    is the increment between the start of each data vector.
68    ! n       is the length of the data vectors.
69    ! lot     is the number of data vectors.
70    ! isign = +1 for transform from spectral to gridpoint
71    !       = -1 for transform from gridpoint to spectral.
72    !
73    ! ordering of coefficients:
74    ! a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2)
75    ! where b(0)=b(n/2)=0; (n+2) locations required.
76    !
77    ! ordering of data:
78    ! x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required.
79    !
80    ! Vectorization is achieved on cray by doing the transforms
81    ! in parallel.
82    !
83    ! n must be composed of factors 2,3 & 5 but does not have to be even.
84    !
85    ! definition of transforms:
86    !
87    ! isign=+1: x(j)=sum(k=0,.,n-1)(c(k)*exp(2*i*j*k*pi/n))
88    ! where c(k)=a(k)+i*b(k) and c(n-k)=a(k)-i*b(k)
89
90    ! isign=-1: a(k)=(1/n)*sum(j=0,.,n-1)(x(j)*cos(2*j*k*pi/n))
91    ! b(k)=-(1/n)*sum(j=0,.,n-1)(x(j)*sin(2*j*k*pi/n))
92
93    ! calls fortran-versions of fft's  !!!
94    ! dimension a(n),work(n),trigs(n),ifax(1)
95
96
97    IMPLICIT NONE
98
99    !  Scalar arguments
100    INTEGER :: inc, isign, jump, lot, n
101
102    !  Array arguments
103    REAL :: a(*), trigs(*), work(*)
104    INTEGER :: ifax(*)
105
106    !  Local scalars:
107    INTEGER :: i, ia, ibase, ierr, ifac, igo, ii, istart, ix, iz, j, jbase, jj, &
108         &      k, la, nb, nblox, nfax, nvex, nx
109
110    !  Intrinsic functions
111    INTRINSIC MOD
112
113
114    !  Executable statements
115
116    IF (ifax(10)/=n) CALL set99(trigs,ifax,n)
117    nfax = ifax(1)
118    nx = n + 1
119    IF (MOD(n,2)==1) nx = n
120    nblox = 1 + (lot-1)/nfft
121    nvex = lot - (nblox-1)*nfft
122    IF (isign==-1) GO TO 50
123
124    ! isign=+1, spectral to gridpoint transform
125
126    istart = 1
127    DO nb = 1, nblox
128       ia = istart
129       i = istart
130!DIR$ IVDEP
131!CDIR NODEP
132!OCL NOVREC
133       DO j = 1, nvex
134          a(i+inc) = 0.5*a(i)
135          i = i + jump
136       END DO
137       IF (MOD(n,2)==1) GO TO 10
138       i = istart + n*inc
139       DO j = 1, nvex
140          a(i) = 0.5*a(i)
141          i = i + jump
142       END DO
14310     CONTINUE
144       ia = istart + inc
145       la = 1
146       igo = + 1
147
148       DO k = 1, nfax
149          ifac = ifax(k+1)
150          ierr = -1
151          IF (igo==-1) GO TO 20
152          CALL rpassm(a(ia),a(ia+la*inc),work(1),work(ifac*la+1),trigs,inc,1, &
153               &          jump,nx,nvex,n,ifac,la,ierr)
154          GO TO 30
15520        CONTINUE
156          CALL rpassm(work(1),work(la+1),a(ia),a(ia+ifac*la*inc),trigs,1,inc,nx, &
157               &          jump,nvex,n,ifac,la,ierr)
15830        CONTINUE
159          IF (ierr/=0) GO TO 100
160          la = ifac*la
161          igo = -igo
162          ia = istart
163       END DO
164
165       ! If necessary, copy results back to a
166
167       IF (MOD(nfax,2)==0) GO TO 40
168       ibase = 1
169       jbase = ia
170       DO jj = 1, nvex
171          i = ibase
172          j = jbase
173          DO ii = 1, n
174             a(j) = work(i)
175             i = i + 1
176             j = j + inc
177          END DO
178          ibase = ibase + nx
179          jbase = jbase + jump
180       END DO
18140     CONTINUE
182
183       ! Fill in zeros at end
184
185       ix = istart + n*inc
186!DIR$ IVDEP
187!CDIR NODEP
188!OCL NOVREC
189       DO j = 1, nvex
190          a(ix) = 0.0
191          a(ix+inc) = 0.0
192          ix = ix + jump
193       END DO
194
195       istart = istart + nvex*jump
196       nvex = nfft
197    END DO
198    RETURN
199
200    ! isign=-1, gridpoint to spectral transform
201
20250  CONTINUE
203    istart = 1
204    DO nb = 1, nblox
205       ia = istart
206       la = n
207       igo = + 1
208
209       DO k = 1, nfax
210          ifac = ifax(nfax+2-k)
211          la = la/ifac
212          ierr = -1
213          IF (igo==-1) GO TO 60
214          CALL qpassm(a(ia),a(ia+ifac*la*inc),work(1),work(la+1),trigs,inc,1, &
215               &          jump,nx,nvex,n,ifac,la,ierr)
216          GO TO 70
21760        CONTINUE
218          CALL qpassm(work(1),work(ifac*la+1),a(ia),a(ia+la*inc),trigs,1,inc,nx, &
219               &          jump,nvex,n,ifac,la,ierr)
22070        CONTINUE
221          IF (ierr/=0) GO TO 100
222          igo = -igo
223          ia = istart + inc
224       END DO
225
226       ! If necessary, copy results back to a
227
228       IF (MOD(nfax,2)==0) GO TO 80
229       ibase = 1
230       jbase = ia
231       DO jj = 1, nvex
232          i = ibase
233          j = jbase
234          DO ii = 1, n
235             a(j) = work(i)
236             i = i + 1
237             j = j + inc
238          END DO
239          ibase = ibase + nx
240          jbase = jbase + jump
241       END DO
24280     CONTINUE
243
244       ! Shift a(0) & fill in zero imag parts
245
246       ix = istart
247!DIR$ IVDEP
248!CDIR NODEP
249!OCL NOVREC
250       DO j = 1, nvex
251          a(ix) = a(ix+inc)
252          a(ix+inc) = 0.0
253          ix = ix + jump
254       END DO
255       IF (MOD(n,2)==1) GO TO 90
256       iz = istart + (n+1)*inc
257       DO j = 1, nvex
258          a(iz) = 0.0
259          iz = iz + jump
260       END DO
26190     CONTINUE
262
263       istart = istart + nvex*jump
264       nvex = nfft
265    END DO
266    RETURN
267
268    ! Error messages
269
270100 CONTINUE
271
272    SELECT CASE (ierr)
273    CASE (:-1)
274       WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
275    CASE (0)
276       WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
277    CASE (1:)
278       WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
279    END SELECT
280
281    RETURN
282  END SUBROUTINE fft991cy
283
284  SUBROUTINE qpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
285
286    ! Description:
287    !
288    ! Performs one pass through data as part of
289    ! multiple real fft (fourier analysis) routine.
290    !
291    ! Method:
292    !
293    ! a       is first real input vector
294    !         equivalence b(1) with a(ifac*la*inc1+1)
295    ! c       is first real output vector
296    !         equivalence d(1) with c(la*inc2+1)
297    ! trigs   is a precalculated list of sines & cosines
298    ! inc1    is the addressing increment for a
299    ! inc2    is the addressing increment for c
300    ! inc3    is the increment between input vectors a
301    ! inc4    is the increment between output vectors c
302    ! lot     is the number of vectors
303    ! n       is the length of the vectors
304    ! ifac    is the current factor of n
305    !         la = n/(product of factors used so far)
306    ! ierr    is an error indicator:
307    !         0 - pass completed without error
308    !         1 - lot greater than nfft
309    !         2 - ifac not catered for
310    !         3 - ifac only catered for if la=n/ifac
311    !
312
313    IMPLICIT NONE
314
315    !  Scalar arguments
316    INTEGER :: ierr, ifac, inc1, inc2, inc3, inc4, la, lot, n
317
318    !  Array arguments
319    ! REAL :: a(n),b(n),c(n),d(n),trigs(n)
320    REAL :: a(*), b(*), c(*), d(*), trigs(*)
321
322    !  Local scalars:
323    REAL :: a0, a1, a10, a11, a2, a20, a21, a3, a4, a5, a6, b0, b1, b10, b11, &
324         &      b2, b20, b21, b3, b4, b5, b6, c1, c2, c3, c4, c5, qrt5, s1, s2, s3, s4, &
325         &      s5, sin36, sin45, sin60, sin72, z, zqrt5, zsin36, zsin45, zsin60, &
326         &      zsin72
327    INTEGER :: i, ia, ib, ibad, ibase, ic, id, ie, if, ig, igo, ih, iink, ijk, &
328         &      ijump, j, ja, jb, jbase, jc, jd, je, jf, jink, k, kb, kc, kd, ke, kf, &
329         &      kstop, l, m
330
331    !  Intrinsic functions
332    INTRINSIC REAL, SQRT
333
334    !  Data statements
335    DATA sin36/0.587785252292473/, sin72/0.951056516295154/, &
336         &      qrt5/0.559016994374947/, sin60/0.866025403784437/
337
338
339    !  Executable statements
340
341    m = n/ifac
342    iink = la*inc1
343    jink = la*inc2
344    ijump = (ifac-1)*iink
345    kstop = (n-ifac)/(2*ifac)
346
347    ibad = 1
348    IF (lot>nfft) GO TO 180
349    ibase = 0
350    jbase = 0
351    igo = ifac - 1
352    IF (igo==7) igo = 6
353    ibad = 2
354    IF (igo<1 .OR. igo>6) GO TO 180
355    GO TO (10,40,70,100,130,160) igo
356
357    ! Coding for factor 2
358
35910  CONTINUE
360    ia = 1
361    ib = ia + iink
362    ja = 1
363    jb = ja + (2*m-la)*inc2
364
365    IF (la==m) GO TO 30
366
367    DO l = 1, la
368       i = ibase
369       j = jbase
370!DIR$ IVDEP
371!CDIR NODEP
372!OCL NOVREC
373       DO ijk = 1, lot
374          c(ja+j) = a(ia+i) + a(ib+i)
375          c(jb+j) = a(ia+i) - a(ib+i)
376          i = i + inc3
377          j = j + inc4
378       END DO
379       ibase = ibase + inc1
380       jbase = jbase + inc2
381    END DO
382    ja = ja + jink
383    jink = 2*jink
384    jb = jb - jink
385    ibase = ibase + ijump
386    ijump = 2*ijump + iink
387    IF (ja==jb) GO TO 20
388    DO k = la, kstop, la
389       kb = k + k
390       c1 = trigs(kb+1)
391       s1 = trigs(kb+2)
392       jbase = 0
393       DO l = 1, la
394          i = ibase
395          j = jbase
396!DIR$ IVDEP
397!CDIR NODEP
398!OCL NOVREC
399          DO ijk = 1, lot
400             c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
401             c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
402             d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
403             d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
404             i = i + inc3
405             j = j + inc4
406          END DO
407          ibase = ibase + inc1
408          jbase = jbase + inc2
409       END DO
410       ibase = ibase + ijump
411       ja = ja + jink
412       jb = jb - jink
413    END DO
414    IF (ja>jb) GO TO 170
41520  CONTINUE
416    jbase = 0
417    DO l = 1, la
418       i = ibase
419       j = jbase
420!DIR$ IVDEP
421!CDIR NODEP
422!OCL NOVREC
423       DO ijk = 1, lot
424          c(ja+j) = a(ia+i)
425          d(ja+j) = -a(ib+i)
426          i = i + inc3
427          j = j + inc4
428       END DO
429       ibase = ibase + inc1
430       jbase = jbase + inc2
431    END DO
432
433    GO TO 170
43430  CONTINUE
435    z = 1.0/REAL(n)
436    DO l = 1, la
437       i = ibase
438       j = jbase
439!DIR$ IVDEP
440!CDIR NODEP
441!OCL NOVREC
442       DO ijk = 1, lot
443          c(ja+j) = z*(a(ia+i)+a(ib+i))
444          c(jb+j) = z*(a(ia+i)-a(ib+i))
445          i = i + inc3
446          j = j + inc4
447       END DO
448       ibase = ibase + inc1
449       jbase = jbase + inc2
450    END DO
451    GO TO 170
452
453    ! Coding for factor 3
454
45540  CONTINUE
456    ia = 1
457    ib = ia + iink
458    ic = ib + iink
459    ja = 1
460    jb = ja + (2*m-la)*inc2
461    jc = jb
462
463    IF (la==m) GO TO 60
464
465    DO l = 1, la
466       i = ibase
467       j = jbase
468!DIR$ IVDEP
469!CDIR NODEP
470!OCL NOVREC
471       DO ijk = 1, lot
472          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
473          c(jb+j) = a(ia+i) - 0.5*(a(ib+i)+a(ic+i))
474          d(jb+j) = sin60*(a(ic+i)-a(ib+i))
475          i = i + inc3
476          j = j + inc4
477       END DO
478       ibase = ibase + inc1
479       jbase = jbase + inc2
480    END DO
481    ja = ja + jink
482    jink = 2*jink
483    jb = jb + jink
484    jc = jc - jink
485    ibase = ibase + ijump
486    ijump = 2*ijump + iink
487    IF (ja==jc) GO TO 50
488    DO k = la, kstop, la
489       kb = k + k
490       kc = kb + kb
491       c1 = trigs(kb+1)
492       s1 = trigs(kb+2)
493       c2 = trigs(kc+1)
494       s2 = trigs(kc+2)
495       jbase = 0
496       DO l = 1, la
497          i = ibase
498          j = jbase
499!DIR$ IVDEP
500!CDIR NODEP
501!OCL NOVREC
502          DO ijk = 1, lot
503             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
504             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
505             a2 = a(ia+i) - 0.5*a1
506             b2 = b(ia+i) - 0.5*b1
507             a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
508             b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
509             c(ja+j) = a(ia+i) + a1
510             d(ja+j) = b(ia+i) + b1
511             c(jb+j) = a2 + b3
512             d(jb+j) = b2 - a3
513             c(jc+j) = a2 - b3
514             d(jc+j) = -(b2+a3)
515             i = i + inc3
516             j = j + inc4
517          END DO
518          ibase = ibase + inc1
519          jbase = jbase + inc2
520       END DO
521       ibase = ibase + ijump
522       ja = ja + jink
523       jb = jb + jink
524       jc = jc - jink
525    END DO
526    IF (ja>jc) GO TO 170
52750  CONTINUE
528    jbase = 0
529    DO l = 1, la
530       i = ibase
531       j = jbase
532!DIR$ IVDEP
533!CDIR NODEP
534!OCL NOVREC
535       DO ijk = 1, lot
536          c(ja+j) = a(ia+i) + 0.5*(a(ib+i)-a(ic+i))
537          d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
538          c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
539          i = i + inc3
540          j = j + inc4
541       END DO
542       ibase = ibase + inc1
543       jbase = jbase + inc2
544    END DO
545
546    GO TO 170
54760  CONTINUE
548    z = 1.0/REAL(n)
549    zsin60 = z*sin60
550    DO l = 1, la
551       i = ibase
552       j = jbase
553!DIR$ IVDEP
554!CDIR NODEP
555!OCL NOVREC
556       DO ijk = 1, lot
557          c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
558          c(jb+j) = z*(a(ia+i)-0.5*(a(ib+i)+a(ic+i)))
559          d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
560          i = i + inc3
561          j = j + inc4
562       END DO
563       ibase = ibase + inc1
564       jbase = jbase + inc2
565    END DO
566    GO TO 170
567
568    ! Coding for factor 4
569
57070  CONTINUE
571    ia = 1
572    ib = ia + iink
573    ic = ib + iink
574    id = ic + iink
575    ja = 1
576    jb = ja + (2*m-la)*inc2
577    jc = jb + 2*m*inc2
578    jd = jb
579
580    IF (la==m) GO TO 90
581
582    DO l = 1, la
583       i = ibase
584       j = jbase
585!DIR$ IVDEP
586!CDIR NODEP
587!OCL NOVREC
588       DO ijk = 1, lot
589          c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
590          c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
591          c(jb+j) = a(ia+i) - a(ic+i)
592          d(jb+j) = a(id+i) - a(ib+i)
593          i = i + inc3
594          j = j + inc4
595       END DO
596       ibase = ibase + inc1
597       jbase = jbase + inc2
598    END DO
599    ja = ja + jink
600    jink = 2*jink
601    jb = jb + jink
602    jc = jc - jink
603    jd = jd - jink
604    ibase = ibase + ijump
605    ijump = 2*ijump + iink
606    IF (jb==jc) GO TO 80
607    DO k = la, kstop, la
608       kb = k + k
609       kc = kb + kb
610       kd = kc + kb
611       c1 = trigs(kb+1)
612       s1 = trigs(kb+2)
613       c2 = trigs(kc+1)
614       s2 = trigs(kc+2)
615       c3 = trigs(kd+1)
616       s3 = trigs(kd+2)
617       jbase = 0
618       DO l = 1, la
619          i = ibase
620          j = jbase
621!DIR$ IVDEP
622!CDIR NODEP
623!OCL NOVREC
624          DO ijk = 1, lot
625             a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
626             a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
627             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
628             a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
629             b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
630             b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
631             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
632             b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
633             c(ja+j) = a0 + a1
634             c(jc+j) = a0 - a1
635             d(ja+j) = b0 + b1
636             d(jc+j) = b1 - b0
637             c(jb+j) = a2 + b3
638             c(jd+j) = a2 - b3
639             d(jb+j) = b2 - a3
640             d(jd+j) = -(b2+a3)
641             i = i + inc3
642             j = j + inc4
643          END DO
644          ibase = ibase + inc1
645          jbase = jbase + inc2
646       END DO
647       ibase = ibase + ijump
648       ja = ja + jink
649       jb = jb + jink
650       jc = jc - jink
651       jd = jd - jink
652    END DO
653    IF (jb>jc) GO TO 170
65480  CONTINUE
655    sin45 = SQRT(0.5)
656    jbase = 0
657    DO l = 1, la
658       i = ibase
659       j = jbase
660!DIR$ IVDEP
661!CDIR NODEP
662!OCL NOVREC
663       DO ijk = 1, lot
664          c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
665          c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
666          d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
667          d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
668          i = i + inc3
669          j = j + inc4
670       END DO
671       ibase = ibase + inc1
672       jbase = jbase + inc2
673    END DO
674
675    GO TO 170
67690  CONTINUE
677    z = 1.0/REAL(n)
678    DO l = 1, la
679       i = ibase
680       j = jbase
681!DIR$ IVDEP
682!CDIR NODEP
683!OCL NOVREC
684       DO ijk = 1, lot
685          c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
686          c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
687          c(jb+j) = z*(a(ia+i)-a(ic+i))
688          d(jb+j) = z*(a(id+i)-a(ib+i))
689          i = i + inc3
690          j = j + inc4
691       END DO
692       ibase = ibase + inc1
693       jbase = jbase + inc2
694    END DO
695    GO TO 170
696
697    ! Coding for factor 5
698
699100 CONTINUE
700    ia = 1
701    ib = ia + iink
702    ic = ib + iink
703    id = ic + iink
704    ie = id + iink
705    ja = 1
706    jb = ja + (2*m-la)*inc2
707    jc = jb + 2*m*inc2
708    jd = jc
709    je = jb
710
711    IF (la==m) GO TO 120
712
713    DO l = 1, la
714       i = ibase
715       j = jbase
716!DIR$ IVDEP
717!CDIR NODEP
718!OCL NOVREC
719       DO ijk = 1, lot
720          a1 = a(ib+i) + a(ie+i)
721          a3 = a(ib+i) - a(ie+i)
722          a2 = a(ic+i) + a(id+i)
723          a4 = a(ic+i) - a(id+i)
724          a5 = a(ia+i) - 0.25*(a1+a2)
725          a6 = qrt5*(a1-a2)
726          c(ja+j) = a(ia+i) + (a1+a2)
727          c(jb+j) = a5 + a6
728          c(jc+j) = a5 - a6
729          d(jb+j) = -sin72*a3 - sin36*a4
730          d(jc+j) = -sin36*a3 + sin72*a4
731          i = i + inc3
732          j = j + inc4
733       END DO
734       ibase = ibase + inc1
735       jbase = jbase + inc2
736    END DO
737    ja = ja + jink
738    jink = 2*jink
739    jb = jb + jink
740    jc = jc + jink
741    jd = jd - jink
742    je = je - jink
743    ibase = ibase + ijump
744    ijump = 2*ijump + iink
745    IF (jb==jd) GO TO 110
746    DO k = la, kstop, la
747       kb = k + k
748       kc = kb + kb
749       kd = kc + kb
750       ke = kd + kb
751       c1 = trigs(kb+1)
752       s1 = trigs(kb+2)
753       c2 = trigs(kc+1)
754       s2 = trigs(kc+2)
755       c3 = trigs(kd+1)
756       s3 = trigs(kd+2)
757       c4 = trigs(ke+1)
758       s4 = trigs(ke+2)
759       jbase = 0
760       DO l = 1, la
761          i = ibase
762          j = jbase
763!DIR$ IVDEP
764!CDIR NODEP
765!OCL NOVREC
766          DO ijk = 1, lot
767             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
768             a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
769             a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
770             a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
771             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
772             b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
773             b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
774             b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
775             a5 = a(ia+i) - 0.25*(a1+a2)
776             a6 = qrt5*(a1-a2)
777             b5 = b(ia+i) - 0.25*(b1+b2)
778             b6 = qrt5*(b1-b2)
779             a10 = a5 + a6
780             a20 = a5 - a6
781             b10 = b5 + b6
782             b20 = b5 - b6
783             a11 = sin72*b3 + sin36*b4
784             a21 = sin36*b3 - sin72*b4
785             b11 = sin72*a3 + sin36*a4
786             b21 = sin36*a3 - sin72*a4
787             c(ja+j) = a(ia+i) + (a1+a2)
788             c(jb+j) = a10 + a11
789             c(je+j) = a10 - a11
790             c(jc+j) = a20 + a21
791             c(jd+j) = a20 - a21
792             d(ja+j) = b(ia+i) + (b1+b2)
793             d(jb+j) = b10 - b11
794             d(je+j) = -(b10+b11)
795             d(jc+j) = b20 - b21
796             d(jd+j) = -(b20+b21)
797             i = i + inc3
798             j = j + inc4
799          END DO
800          ibase = ibase + inc1
801          jbase = jbase + inc2
802       END DO
803       ibase = ibase + ijump
804       ja = ja + jink
805       jb = jb + jink
806       jc = jc + jink
807       jd = jd - jink
808       je = je - jink
809    END DO
810    IF (jb>jd) GO TO 170
811110 CONTINUE
812    jbase = 0
813    DO l = 1, la
814       i = ibase
815       j = jbase
816!DIR$ IVDEP
817!CDIR NODEP
818!OCL NOVREC
819       DO ijk = 1, lot
820          a1 = a(ib+i) + a(ie+i)
821          a3 = a(ib+i) - a(ie+i)
822          a2 = a(ic+i) + a(id+i)
823          a4 = a(ic+i) - a(id+i)
824          a5 = a(ia+i) + 0.25*(a3-a4)
825          a6 = qrt5*(a3+a4)
826          c(ja+j) = a5 + a6
827          c(jb+j) = a5 - a6
828          c(jc+j) = a(ia+i) - (a3-a4)
829          d(ja+j) = -sin36*a1 - sin72*a2
830          d(jb+j) = -sin72*a1 + sin36*a2
831          i = i + inc3
832          j = j + inc4
833       END DO
834       ibase = ibase + inc1
835       jbase = jbase + inc2
836    END DO
837
838    GO TO 170
839120 CONTINUE
840    z = 1.0/REAL(n)
841    zqrt5 = z*qrt5
842    zsin36 = z*sin36
843    zsin72 = z*sin72
844    DO l = 1, la
845       i = ibase
846       j = jbase
847!DIR$ IVDEP
848!CDIR NODEP
849!OCL NOVREC
850       DO ijk = 1, lot
851          a1 = a(ib+i) + a(ie+i)
852          a3 = a(ib+i) - a(ie+i)
853          a2 = a(ic+i) + a(id+i)
854          a4 = a(ic+i) - a(id+i)
855          a5 = z*(a(ia+i)-0.25*(a1+a2))
856          a6 = zqrt5*(a1-a2)
857          c(ja+j) = z*(a(ia+i)+(a1+a2))
858          c(jb+j) = a5 + a6
859          c(jc+j) = a5 - a6
860          d(jb+j) = -zsin72*a3 - zsin36*a4
861          d(jc+j) = -zsin36*a3 + zsin72*a4
862          i = i + inc3
863          j = j + inc4
864       END DO
865       ibase = ibase + inc1
866       jbase = jbase + inc2
867    END DO
868    GO TO 170
869
870    ! Coding for factor 6
871
872130 CONTINUE
873    ia = 1
874    ib = ia + iink
875    ic = ib + iink
876    id = ic + iink
877    ie = id + iink
878    if = ie + iink
879    ja = 1
880    jb = ja + (2*m-la)*inc2
881    jc = jb + 2*m*inc2
882    jd = jc + 2*m*inc2
883    je = jc
884    jf = jb
885
886    IF (la==m) GO TO 150
887
888    DO l = 1, la
889       i = ibase
890       j = jbase
891!DIR$ IVDEP
892!CDIR NODEP
893!OCL NOVREC
894       DO ijk = 1, lot
895          a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
896          c(ja+j) = (a(ia+i)+a(id+i)) + a11
897          c(jc+j) = (a(ia+i)+a(id+i)-0.5*a11)
898          d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
899          a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
900          c(jb+j) = (a(ia+i)-a(id+i)) - 0.5*a11
901          d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
902          c(jd+j) = (a(ia+i)-a(id+i)) + a11
903          i = i + inc3
904          j = j + inc4
905       END DO
906       ibase = ibase + inc1
907       jbase = jbase + inc2
908    END DO
909    ja = ja + jink
910    jink = 2*jink
911    jb = jb + jink
912    jc = jc + jink
913    jd = jd - jink
914    je = je - jink
915    jf = jf - jink
916    ibase = ibase + ijump
917    ijump = 2*ijump + iink
918    IF (jc==jd) GO TO 140
919    DO k = la, kstop, la
920       kb = k + k
921       kc = kb + kb
922       kd = kc + kb
923       ke = kd + kb
924       kf = ke + kb
925       c1 = trigs(kb+1)
926       s1 = trigs(kb+2)
927       c2 = trigs(kc+1)
928       s2 = trigs(kc+2)
929       c3 = trigs(kd+1)
930       s3 = trigs(kd+2)
931       c4 = trigs(ke+1)
932       s4 = trigs(ke+2)
933       c5 = trigs(kf+1)
934       s5 = trigs(kf+2)
935       jbase = 0
936       DO l = 1, la
937          i = ibase
938          j = jbase
939!DIR$ IVDEP
940!CDIR NODEP
941!OCL NOVREC
942          DO ijk = 1, lot
943             a1 = c1*a(ib+i) + s1*b(ib+i)
944             b1 = c1*b(ib+i) - s1*a(ib+i)
945             a2 = c2*a(ic+i) + s2*b(ic+i)
946             b2 = c2*b(ic+i) - s2*a(ic+i)
947             a3 = c3*a(id+i) + s3*b(id+i)
948             b3 = c3*b(id+i) - s3*a(id+i)
949             a4 = c4*a(ie+i) + s4*b(ie+i)
950             b4 = c4*b(ie+i) - s4*a(ie+i)
951             a5 = c5*a(if+i) + s5*b(if+i)
952             b5 = c5*b(if+i) - s5*a(if+i)
953             a11 = (a2+a5) + (a1+a4)
954             a20 = (a(ia+i)+a3) - 0.5*a11
955             a21 = sin60*((a2+a5)-(a1+a4))
956             b11 = (b2+b5) + (b1+b4)
957             b20 = (b(ia+i)+b3) - 0.5*b11
958             b21 = sin60*((b2+b5)-(b1+b4))
959             c(ja+j) = (a(ia+i)+a3) + a11
960             d(ja+j) = (b(ia+i)+b3) + b11
961             c(jc+j) = a20 - b21
962             d(jc+j) = a21 + b20
963             c(je+j) = a20 + b21
964             d(je+j) = a21 - b20
965             a11 = (a2-a5) + (a4-a1)
966             a20 = (a(ia+i)-a3) - 0.5*a11
967             a21 = sin60*((a4-a1)-(a2-a5))
968             b11 = (b5-b2) - (b4-b1)
969             b20 = (b3-b(ia+i)) - 0.5*b11
970             b21 = sin60*((b5-b2)+(b4-b1))
971             c(jb+j) = a20 - b21
972             d(jb+j) = a21 - b20
973             c(jd+j) = a11 + (a(ia+i)-a3)
974             d(jd+j) = b11 + (b3-b(ia+i))
975             c(jf+j) = a20 + b21
976             d(jf+j) = a21 + b20
977             i = i + inc3
978             j = j + inc4
979          END DO
980          ibase = ibase + inc1
981          jbase = jbase + inc2
982       END DO
983       ibase = ibase + ijump
984       ja = ja + jink
985       jb = jb + jink
986       jc = jc + jink
987       jd = jd - jink
988       je = je - jink
989       jf = jf - jink
990    END DO
991    IF (jc>jd) GO TO 170
992140 CONTINUE
993    jbase = 0
994    DO l = 1, la
995       i = ibase
996       j = jbase
997!DIR$ IVDEP
998!CDIR NODEP
999!OCL NOVREC
1000       DO ijk = 1, lot
1001          c(ja+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
1002          d(ja+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
1003          c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
1004          d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
1005          c(jc+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
1006          d(jc+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
1007          i = i + inc3
1008          j = j + inc4
1009       END DO
1010       ibase = ibase + inc1
1011       jbase = jbase + inc2
1012    END DO
1013
1014    GO TO 170
1015150 CONTINUE
1016    z = 1.0/REAL(n)
1017    zsin60 = z*sin60
1018    DO l = 1, la
1019       i = ibase
1020       j = jbase
1021!DIR$ IVDEP
1022!CDIR NODEP
1023!OCL NOVREC
1024       DO ijk = 1, lot
1025          a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
1026          c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
1027          c(jc+j) = z*((a(ia+i)+a(id+i))-0.5*a11)
1028          d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
1029          a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
1030          c(jb+j) = z*((a(ia+i)-a(id+i))-0.5*a11)
1031          d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
1032          c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
1033          i = i + inc3
1034          j = j + inc4
1035       END DO
1036       ibase = ibase + inc1
1037       jbase = jbase + inc2
1038    END DO
1039    GO TO 170
1040
1041    ! Coding for factor 8
1042
1043160 CONTINUE
1044    ibad = 3
1045    IF (la/=m) GO TO 180
1046    ia = 1
1047    ib = ia + iink
1048    ic = ib + iink
1049    id = ic + iink
1050    ie = id + iink
1051    if = ie + iink
1052    ig = if + iink
1053    ih = ig + iink
1054    ja = 1
1055    jb = ja + la*inc2
1056    jc = jb + 2*m*inc2
1057    jd = jc + 2*m*inc2
1058    je = jd + 2*m*inc2
1059    z = 1.0/REAL(n)
1060    zsin45 = z*SQRT(0.5)
1061
1062    DO l = 1, la
1063       i = ibase
1064       j = jbase
1065!DIR$ IVDEP
1066!CDIR NODEP
1067!OCL NOVREC
1068       DO ijk = 1, lot
1069          c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ &
1070               &          a(ih+i))+(a(ib+i)+a(if+i))))
1071          c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ &
1072               &          a(ih+i))+(a(ib+i)+a(if+i))))
1073          c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
1074          d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
1075          c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+ &
1076               &          i)-a(ib+i)))
1077          c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+ &
1078               &          i)-a(ib+i)))
1079          d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + &
1080               &          z*(a(ig+i)-a(ic+i))
1081          d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - &
1082               &          z*(a(ig+i)-a(ic+i))
1083          i = i + inc3
1084          j = j + inc4
1085       END DO
1086       ibase = ibase + inc1
1087       jbase = jbase + inc2
1088    END DO
1089
1090    ! Return
1091
1092170 CONTINUE
1093    ibad = 0
1094180 CONTINUE
1095    ierr = ibad
1096    RETURN
1097  END SUBROUTINE qpassm
1098
1099  SUBROUTINE rpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
1100    ! Dimension a(n),b(n),c(n),d(n),trigs(n)
1101
1102    IMPLICIT NONE
1103
1104    !  Scalar arguments
1105    INTEGER :: ierr, ifac, inc1, inc2, inc3, inc4, la, lot, n
1106
1107    !  Array arguments
1108    REAL :: a(*), b(*), c(*), d(*), trigs(*)
1109
1110    !  Local scalars:
1111    REAL :: c1, c2, c3, c4, c5, qqrt5, qrt5, s1, s2, s3, s4, s5, sin36, sin45, &
1112         &      sin60, sin72, ssin36, ssin45, ssin60, ssin72
1113    INTEGER :: i, ia, ib, ibad, ibase, ic, id, ie, if, igo, iink, ijk, j, ja, &
1114         &      jb, jbase, jc, jd, je, jf, jg, jh, jink, jump, k, kb, kc, kd, ke, kf, &
1115         &      kstop, l, m
1116
1117    !  Local arrays:
1118    REAL :: a10(nfft), a11(nfft), a20(nfft), a21(nfft), b10(nfft), b11(nfft), b20(nfft), &
1119         &      b21(nfft)
1120
1121    !  Intrinsic functions
1122    INTRINSIC SQRT
1123
1124    !  Data statements
1125    DATA sin36/0.587785252292473/, sin72/0.951056516295154/, &
1126         &      qrt5/0.559016994374947/, sin60/0.866025403784437/
1127
1128
1129    !  Executable statements
1130
1131    m = n/ifac
1132    iink = la*inc1
1133    jink = la*inc2
1134    jump = (ifac-1)*jink
1135    kstop = (n-ifac)/(2*ifac)
1136
1137    ibad = 1
1138    IF (lot>nfft) GO TO 180
1139    ibase = 0
1140    jbase = 0
1141    igo = ifac - 1
1142    IF (igo==7) igo = 6
1143    ibad = 2
1144    IF (igo<1 .OR. igo>6) GO TO 180
1145    GO TO (10,40,70,100,130,160) igo
1146
1147    ! Coding for factor 2
1148
114910  CONTINUE
1150    ia = 1
1151    ib = ia + (2*m-la)*inc1
1152    ja = 1
1153    jb = ja + jink
1154
1155    IF (la==m) GO TO 30
1156
1157    DO l = 1, la
1158       i = ibase
1159       j = jbase
1160!DIR$ IVDEP
1161!CDIR NODEP
1162!OCL NOVREC
1163       DO ijk = 1, lot
1164          c(ja+j) = a(ia+i) + a(ib+i)
1165          c(jb+j) = a(ia+i) - a(ib+i)
1166          i = i + inc3
1167          j = j + inc4
1168       END DO
1169       ibase = ibase + inc1
1170       jbase = jbase + inc2
1171    END DO
1172    ia = ia + iink
1173    iink = 2*iink
1174    ib = ib - iink
1175    ibase = 0
1176    jbase = jbase + jump
1177    jump = 2*jump + jink
1178    IF (ia==ib) GO TO 20
1179    DO k = la, kstop, la
1180       kb = k + k
1181       c1 = trigs(kb+1)
1182       s1 = trigs(kb+2)
1183       ibase = 0
1184       DO l = 1, la
1185          i = ibase
1186          j = jbase
1187!DIR$ IVDEP
1188!CDIR NODEP
1189!OCL NOVREC
1190          DO ijk = 1, lot
1191             c(ja+j) = a(ia+i) + a(ib+i)
1192             d(ja+j) = b(ia+i) - b(ib+i)
1193             c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i))
1194             d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i))
1195             i = i + inc3
1196             j = j + inc4
1197          END DO
1198          ibase = ibase + inc1
1199          jbase = jbase + inc2
1200       END DO
1201       ia = ia + iink
1202       ib = ib - iink
1203       jbase = jbase + jump
1204    END DO
1205    IF (ia>ib) GO TO 170
120620  CONTINUE
1207    ibase = 0
1208    DO l = 1, la
1209       i = ibase
1210       j = jbase
1211!DIR$ IVDEP
1212!CDIR NODEP
1213!OCL NOVREC
1214       DO ijk = 1, lot
1215          c(ja+j) = a(ia+i)
1216          c(jb+j) = -b(ia+i)
1217          i = i + inc3
1218          j = j + inc4
1219       END DO
1220       ibase = ibase + inc1
1221       jbase = jbase + inc2
1222    END DO
1223
1224    GO TO 170
122530  CONTINUE
1226    DO l = 1, la
1227       i = ibase
1228       j = jbase
1229!DIR$ IVDEP
1230!CDIR NODEP
1231!OCL NOVREC
1232       DO ijk = 1, lot
1233          c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
1234          c(jb+j) = 2.0*(a(ia+i)-a(ib+i))
1235          i = i + inc3
1236          j = j + inc4
1237       END DO
1238       ibase = ibase + inc1
1239       jbase = jbase + inc2
1240    END DO
1241    GO TO 170
1242
1243    ! Coding for factor 3
1244
124540  CONTINUE
1246    ia = 1
1247    ib = ia + (2*m-la)*inc1
1248    ic = ib
1249    ja = 1
1250    jb = ja + jink
1251    jc = jb + jink
1252
1253    IF (la==m) GO TO 60
1254
1255    DO l = 1, la
1256       i = ibase
1257       j = jbase
1258!DIR$ IVDEP
1259!CDIR NODEP
1260!OCL NOVREC
1261       DO ijk = 1, lot
1262          c(ja+j) = a(ia+i) + a(ib+i)
1263          c(jb+j) = (a(ia+i)-0.5*a(ib+i)) - (sin60*(b(ib+i)))
1264          c(jc+j) = (a(ia+i)-0.5*a(ib+i)) + (sin60*(b(ib+i)))
1265          i = i + inc3
1266          j = j + inc4
1267       END DO
1268       ibase = ibase + inc1
1269       jbase = jbase + inc2
1270    END DO
1271    ia = ia + iink
1272    iink = 2*iink
1273    ib = ib + iink
1274    ic = ic - iink
1275    jbase = jbase + jump
1276    jump = 2*jump + jink
1277    IF (ia==ic) GO TO 50
1278    DO k = la, kstop, la
1279       kb = k + k
1280       kc = kb + kb
1281       c1 = trigs(kb+1)
1282       s1 = trigs(kb+2)
1283       c2 = trigs(kc+1)
1284       s2 = trigs(kc+2)
1285       ibase = 0
1286       DO l = 1, la
1287          i = ibase
1288          j = jbase
1289!DIR$ IVDEP
1290!CDIR NODEP
1291!OCL NOVREC
1292          DO ijk = 1, lot
1293             c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
1294             d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
1295             c(jb+j) = c1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
1296                  &            b(ic+i)))) - s1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
1297                  &            a(ic+i))))
1298             d(jb+j) = s1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
1299                  &            b(ic+i)))) + c1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
1300                  &            a(ic+i))))
1301             c(jc+j) = c2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
1302                  &            b(ic+i)))) - s2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
1303                  &            a(ic+i))))
1304             d(jc+j) = s2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
1305                  &            b(ic+i)))) + c2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
1306                  &            a(ic+i))))
1307             i = i + inc3
1308             j = j + inc4
1309          END DO
1310          ibase = ibase + inc1
1311          jbase = jbase + inc2
1312       END DO
1313       ia = ia + iink
1314       ib = ib + iink
1315       ic = ic - iink
1316       jbase = jbase + jump
1317    END DO
1318    IF (ia>ic) GO TO 170
131950  CONTINUE
1320    ibase = 0
1321    DO l = 1, la
1322       i = ibase
1323       j = jbase
1324!DIR$ IVDEP
1325!CDIR NODEP
1326!OCL NOVREC
1327       DO ijk = 1, lot
1328          c(ja+j) = a(ia+i) + a(ib+i)
1329          c(jb+j) = (0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
1330          c(jc+j) = -(0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
1331          i = i + inc3
1332          j = j + inc4
1333       END DO
1334       ibase = ibase + inc1
1335       jbase = jbase + inc2
1336    END DO
1337
1338    GO TO 170
133960  CONTINUE
1340    ssin60 = 2.0*sin60
1341    DO l = 1, la
1342       i = ibase
1343       j = jbase
1344!DIR$ IVDEP
1345!CDIR NODEP
1346!OCL NOVREC
1347       DO ijk = 1, lot
1348          c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
1349          c(jb+j) = (2.0*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
1350          c(jc+j) = (2.0*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
1351          i = i + inc3
1352          j = j + inc4
1353       END DO
1354       ibase = ibase + inc1
1355       jbase = jbase + inc2
1356    END DO
1357    GO TO 170
1358
1359    ! Coding for factor 4
1360
136170  CONTINUE
1362    ia = 1
1363    ib = ia + (2*m-la)*inc1
1364    ic = ib + 2*m*inc1
1365    id = ib
1366    ja = 1
1367    jb = ja + jink
1368    jc = jb + jink
1369    jd = jc + jink
1370
1371    IF (la==m) GO TO 90
1372
1373    DO l = 1, la
1374       i = ibase
1375       j = jbase
1376!DIR$ IVDEP
1377!CDIR NODEP
1378!OCL NOVREC
1379       DO ijk = 1, lot
1380          c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
1381          c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
1382          c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
1383          c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
1384          i = i + inc3
1385          j = j + inc4
1386       END DO
1387       ibase = ibase + inc1
1388       jbase = jbase + inc2
1389    END DO
1390    ia = ia + iink
1391    iink = 2*iink
1392    ib = ib + iink
1393    ic = ic - iink
1394    id = id - iink
1395    jbase = jbase + jump
1396    jump = 2*jump + jink
1397    IF (ib==ic) GO TO 80
1398    DO k = la, kstop, la
1399       kb = k + k
1400       kc = kb + kb
1401       kd = kc + kb
1402       c1 = trigs(kb+1)
1403       s1 = trigs(kb+2)
1404       c2 = trigs(kc+1)
1405       s2 = trigs(kc+2)
1406       c3 = trigs(kd+1)
1407       s3 = trigs(kd+2)
1408       ibase = 0
1409       DO l = 1, la
1410          i = ibase
1411          j = jbase
1412!DIR$ IVDEP
1413!CDIR NODEP
1414!OCL NOVREC
1415          DO ijk = 1, lot
1416             c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
1417             d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
1418             c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+ &
1419                  &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
1420             d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+ &
1421                  &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
1422             c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+ &
1423                  &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
1424             d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+ &
1425                  &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
1426             c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+ &
1427                  &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
1428             d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+ &
1429                  &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
1430             i = i + inc3
1431             j = j + inc4
1432          END DO
1433          ibase = ibase + inc1
1434          jbase = jbase + inc2
1435       END DO
1436       ia = ia + iink
1437       ib = ib + iink
1438       ic = ic - iink
1439       id = id - iink
1440       jbase = jbase + jump
1441    END DO
1442    IF (ib>ic) GO TO 170
144380  CONTINUE
1444    ibase = 0
1445    sin45 = SQRT(0.5)
1446    DO l = 1, la
1447       i = ibase
1448       j = jbase
1449!DIR$ IVDEP
1450!CDIR NODEP
1451!OCL NOVREC
1452       DO ijk = 1, lot
1453          c(ja+j) = a(ia+i) + a(ib+i)
1454          c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
1455          c(jc+j) = b(ib+i) - b(ia+i)
1456          c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
1457          i = i + inc3
1458          j = j + inc4
1459       END DO
1460       ibase = ibase + inc1
1461       jbase = jbase + inc2
1462    END DO
1463
1464    GO TO 170
146590  CONTINUE
1466    DO l = 1, la
1467       i = ibase
1468       j = jbase
1469!DIR$ IVDEP
1470!CDIR NODEP
1471!OCL NOVREC
1472       DO ijk = 1, lot
1473          c(ja+j) = 2.0*((a(ia+i)+a(ic+i))+a(ib+i))
1474          c(jb+j) = 2.0*((a(ia+i)-a(ic+i))-b(ib+i))
1475          c(jc+j) = 2.0*((a(ia+i)+a(ic+i))-a(ib+i))
1476          c(jd+j) = 2.0*((a(ia+i)-a(ic+i))+b(ib+i))
1477          i = i + inc3
1478          j = j + inc4
1479       END DO
1480       ibase = ibase + inc1
1481       jbase = jbase + inc2
1482    END DO
1483
1484    ! Coding for factor 5
1485
1486    GO TO 170
1487100 CONTINUE
1488    ia = 1
1489    ib = ia + (2*m-la)*inc1
1490    ic = ib + 2*m*inc1
1491    id = ic
1492    ie = ib
1493    ja = 1
1494    jb = ja + jink
1495    jc = jb + jink
1496    jd = jc + jink
1497    je = jd + jink
1498
1499    IF (la==m) GO TO 120
1500
1501    DO l = 1, la
1502       i = ibase
1503       j = jbase
1504!DIR$ IVDEP
1505!CDIR NODEP
1506!OCL NOVREC
1507       DO ijk = 1, lot
1508          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
1509          c(jb+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
1510               &          (sin72*b(ib+i)+sin36*b(ic+i))
1511          c(jc+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
1512               &          (sin36*b(ib+i)-sin72*b(ic+i))
1513          c(jd+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
1514               &          (sin36*b(ib+i)-sin72*b(ic+i))
1515          c(je+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
1516               &          (sin72*b(ib+i)+sin36*b(ic+i))
1517          i = i + inc3
1518          j = j + inc4
1519       END DO
1520       ibase = ibase + inc1
1521       jbase = jbase + inc2
1522    END DO
1523    ia = ia + iink
1524    iink = 2*iink
1525    ib = ib + iink
1526    ic = ic + iink
1527    id = id - iink
1528    ie = ie - iink
1529    jbase = jbase + jump
1530    jump = 2*jump + jink
1531    IF (ib==id) GO TO 110
1532    DO k = la, kstop, la
1533       kb = k + k
1534       kc = kb + kb
1535       kd = kc + kb
1536       ke = kd + kb
1537       c1 = trigs(kb+1)
1538       s1 = trigs(kb+2)
1539       c2 = trigs(kc+1)
1540       s2 = trigs(kc+2)
1541       c3 = trigs(kd+1)
1542       s3 = trigs(kd+2)
1543       c4 = trigs(ke+1)
1544       s4 = trigs(ke+2)
1545       ibase = 0
1546       DO l = 1, la
1547          i = ibase
1548          j = jbase
1549!DIR$ IVDEP
1550!CDIR NODEP
1551!OCL NOVREC
1552          DO ijk = 1, lot
1553
1554             a10(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
1555                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1556             a20(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
1557                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1558             b10(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
1559                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
1560             b20(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
1561                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
1562             a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
1563             a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
1564             b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
1565             b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
1566
1567             c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
1568             d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
1569             c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
1570             d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
1571             c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
1572             d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
1573             c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
1574             d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
1575             c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
1576             d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
1577
1578             i = i + inc3
1579             j = j + inc4
1580          END DO
1581          ibase = ibase + inc1
1582          jbase = jbase + inc2
1583       END DO
1584       ia = ia + iink
1585       ib = ib + iink
1586       ic = ic + iink
1587       id = id - iink
1588       ie = ie - iink
1589       jbase = jbase + jump
1590    END DO
1591    IF (ib>id) GO TO 170
1592110 CONTINUE
1593    ibase = 0
1594    DO l = 1, la
1595       i = ibase
1596       j = jbase
1597!DIR$ IVDEP
1598!CDIR NODEP
1599!OCL NOVREC
1600       DO ijk = 1, lot
1601          c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
1602          c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1603               &          (sin36*b(ia+i)+sin72*b(ib+i))
1604          c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1605               &          (sin36*b(ia+i)+sin72*b(ib+i))
1606          c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1607               &          (sin72*b(ia+i)-sin36*b(ib+i))
1608          c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1609               &          (sin72*b(ia+i)-sin36*b(ib+i))
1610          i = i + inc3
1611          j = j + inc4
1612       END DO
1613       ibase = ibase + inc1
1614       jbase = jbase + inc2
1615    END DO
1616
1617    GO TO 170
1618120 CONTINUE
1619    qqrt5 = 2.0*qrt5
1620    ssin36 = 2.0*sin36
1621    ssin72 = 2.0*sin72
1622    DO l = 1, la
1623       i = ibase
1624       j = jbase
1625!DIR$ IVDEP
1626!CDIR NODEP
1627!OCL NOVREC
1628       DO ijk = 1, lot
1629          c(ja+j) = 2.0*(a(ia+i)+(a(ib+i)+a(ic+i)))
1630          c(jb+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
1631               &          i))) - (ssin72*b(ib+i)+ssin36*b(ic+i))
1632          c(jc+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
1633               &          i))) - (ssin36*b(ib+i)-ssin72*b(ic+i))
1634          c(jd+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
1635               &          i))) + (ssin36*b(ib+i)-ssin72*b(ic+i))
1636          c(je+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
1637               &          i))) + (ssin72*b(ib+i)+ssin36*b(ic+i))
1638          i = i + inc3
1639          j = j + inc4
1640       END DO
1641       ibase = ibase + inc1
1642       jbase = jbase + inc2
1643    END DO
1644    GO TO 170
1645
1646    ! Coding for factor 6
1647
1648130 CONTINUE
1649    ia = 1
1650    ib = ia + (2*m-la)*inc1
1651    ic = ib + 2*m*inc1
1652    id = ic + 2*m*inc1
1653    ie = ic
1654    if = ib
1655    ja = 1
1656    jb = ja + jink
1657    jc = jb + jink
1658    jd = jc + jink
1659    je = jd + jink
1660    jf = je + jink
1661
1662    IF (la==m) GO TO 150
1663
1664    DO l = 1, la
1665       i = ibase
1666       j = jbase
1667!DIR$ IVDEP
1668!CDIR NODEP
1669!OCL NOVREC
1670       DO ijk = 1, lot
1671          c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
1672          c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
1673          c(jb+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
1674               &          i)+b(ic+i)))
1675          c(jf+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
1676               &          i)+b(ic+i)))
1677          c(jc+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
1678               &          i)-b(ic+i)))
1679          c(je+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
1680               &          i)-b(ic+i)))
1681          i = i + inc3
1682          j = j + inc4
1683       END DO
1684       ibase = ibase + inc1
1685       jbase = jbase + inc2
1686    END DO
1687    ia = ia + iink
1688    iink = 2*iink
1689    ib = ib + iink
1690    ic = ic + iink
1691    id = id - iink
1692    ie = ie - iink
1693    if = if - iink
1694    jbase = jbase + jump
1695    jump = 2*jump + jink
1696    IF (ic==id) GO TO 140
1697    DO k = la, kstop, la
1698       kb = k + k
1699       kc = kb + kb
1700       kd = kc + kb
1701       ke = kd + kb
1702       kf = ke + kb
1703       c1 = trigs(kb+1)
1704       s1 = trigs(kb+2)
1705       c2 = trigs(kc+1)
1706       s2 = trigs(kc+2)
1707       c3 = trigs(kd+1)
1708       s3 = trigs(kd+2)
1709       c4 = trigs(ke+1)
1710       s4 = trigs(ke+2)
1711       c5 = trigs(kf+1)
1712       s5 = trigs(kf+2)
1713       ibase = 0
1714       DO l = 1, la
1715          i = ibase
1716          j = jbase
1717!DIR$ IVDEP
1718!CDIR NODEP
1719!OCL NOVREC
1720          DO ijk = 1, lot
1721
1722             a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
1723             a20(ijk) = (a(ia+i)+a(id+i)) - 0.5*a11(ijk)
1724             a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
1725             b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
1726             b20(ijk) = (b(ia+i)-b(id+i)) - 0.5*b11(ijk)
1727             b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
1728
1729             c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
1730             d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
1731             c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
1732             d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
1733             c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
1734             d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
1735
1736             a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
1737             b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
1738             a20(ijk) = (a(ia+i)-a(id+i)) - 0.5*a11(ijk)
1739             a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
1740             b20(ijk) = (b(ia+i)+b(id+i)) + 0.5*b11(ijk)
1741             b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
1742
1743             c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+ &
1744                  &            i))-b11(ijk))
1745             d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+ &
1746                  &            i))-b11(ijk))
1747             c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
1748             d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
1749             c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
1750             d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
1751
1752             i = i + inc3
1753             j = j + inc4
1754          END DO
1755          ibase = ibase + inc1
1756          jbase = jbase + inc2
1757       END DO
1758       ia = ia + iink
1759       ib = ib + iink
1760       ic = ic + iink
1761       id = id - iink
1762       ie = ie - iink
1763       if = if - iink
1764       jbase = jbase + jump
1765    END DO
1766    IF (ic>id) GO TO 170
1767140 CONTINUE
1768    ibase = 0
1769    DO l = 1, la
1770       i = ibase
1771       j = jbase
1772!DIR$ IVDEP
1773!CDIR NODEP
1774!OCL NOVREC
1775       DO ijk = 1, lot
1776          c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
1777          c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
1778          c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
1779          c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
1780          c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
1781          c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
1782          i = i + inc3
1783          j = j + inc4
1784       END DO
1785       ibase = ibase + inc1
1786       jbase = jbase + inc2
1787    END DO
1788
1789    GO TO 170
1790150 CONTINUE
1791    ssin60 = 2.0*sin60
1792    DO l = 1, la
1793       i = ibase
1794       j = jbase
1795!DIR$ IVDEP
1796!CDIR NODEP
1797!OCL NOVREC
1798       DO ijk = 1, lot
1799          c(ja+j) = (2.0*(a(ia+i)+a(id+i))) + (2.0*(a(ib+i)+a(ic+i)))
1800          c(jd+j) = (2.0*(a(ia+i)-a(id+i))) - (2.0*(a(ib+i)-a(ic+i)))
1801          c(jb+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
1802               &          i)+b(ic+i)))
1803          c(jf+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
1804               &          i)+b(ic+i)))
1805          c(jc+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
1806               &          i)-b(ic+i)))
1807          c(je+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
1808               &          i)-b(ic+i)))
1809          i = i + inc3
1810          j = j + inc4
1811       END DO
1812       ibase = ibase + inc1
1813       jbase = jbase + inc2
1814    END DO
1815    GO TO 170
1816
1817    ! Coding for factor 8
1818
1819160 CONTINUE
1820    ibad = 3
1821    IF (la/=m) GO TO 180
1822    ia = 1
1823    ib = ia + la*inc1
1824    ic = ib + 2*la*inc1
1825    id = ic + 2*la*inc1
1826    ie = id + 2*la*inc1
1827    ja = 1
1828    jb = ja + jink
1829    jc = jb + jink
1830    jd = jc + jink
1831    je = jd + jink
1832    jf = je + jink
1833    jg = jf + jink
1834    jh = jg + jink
1835    ssin45 = SQRT(2.0)
1836
1837    DO l = 1, la
1838       i = ibase
1839       j = jbase
1840!DIR$ IVDEP
1841!CDIR NODEP
1842!OCL NOVREC
1843       DO ijk = 1, lot
1844          c(ja+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
1845          c(je+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
1846          c(jc+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
1847          c(jg+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
1848          c(jb+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
1849               &          i))-(b(ib+i)+b(id+i)))
1850          c(jf+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
1851               &          i))-(b(ib+i)+b(id+i)))
1852          c(jd+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
1853               &          i))+(b(ib+i)+b(id+i)))
1854          c(jh+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
1855               &          i))+(b(ib+i)+b(id+i)))
1856          i = i + inc3
1857          j = j + inc4
1858       END DO
1859       ibase = ibase + inc1
1860       jbase = jbase + inc2
1861    END DO
1862
1863    ! Return
1864
1865170 CONTINUE
1866    ibad = 0
1867180 CONTINUE
1868    ierr = ibad
1869    RETURN
1870  END SUBROUTINE rpassm
1871
1872  SUBROUTINE set99(trigs,ifax,n)
1873
1874    ! Description:
1875    !
1876    ! Computes factors of n & trigonometric functins required by fft99 & fft991cy
1877    !
1878    ! Method:
1879    !
1880    ! Dimension trigs(n),ifax(1),jfax(10),lfax(7)
1881    !
1882    ! subroutine 'set99' - computes factors of n & trigonometric
1883    ! functins required by fft99 & fft991cy
1884
1885
1886    USE pegrid
1887
1888    IMPLICIT NONE
1889
1890    !  Scalar arguments
1891    INTEGER :: n
1892
1893    !  Array arguments
1894    REAL :: trigs(*)
1895    INTEGER :: ifax(*)
1896
1897    !  Local scalars:
1898    REAL :: angle, del
1899    INTEGER :: i, ifac, ixxx, k, l, nfax, nhl, nil, nu
1900
1901    !  Local arrays:
1902    INTEGER :: jfax(10), lfax(7)
1903
1904    !  Intrinsic functions
1905    INTRINSIC ASIN, COS, MOD, REAL, SIN
1906
1907    !  Data statements
1908    DATA lfax/6, 8, 5, 4, 3, 2, 1/
1909
1910
1911    !  Executable statements
1912    ixxx = 1
1913
1914    del = 4.0*ASIN(1.0)/REAL(n)
1915    nil = 0
1916    nhl = (n/2) - 1
1917    DO k = nil, nhl
1918       angle = REAL(k)*del
1919       trigs(2*k+1) = COS(angle)
1920       trigs(2*k+2) = SIN(angle)
1921    END DO
1922
1923    ! Find factors of n (8,6,5,4,3,2; only one 8 allowed)
1924    ! Look for sixes first, store factors in descending order
1925    nu = n
1926    ifac = 6
1927    k = 0
1928    l = 1
192910  CONTINUE
1930    IF (MOD(nu,ifac)/=0) GO TO 30
1931    k = k + 1
1932    jfax(k) = ifac
1933    IF (ifac/=8) GO TO 20
1934    IF (k==1) GO TO 20
1935    jfax(1) = 8
1936    jfax(k) = 6
193720  CONTINUE
1938    nu = nu/ifac
1939    IF (nu==1) GO TO 40
1940    IF (ifac/=8) GO TO 10
194130  CONTINUE
1942    l = l + 1
1943    ifac = lfax(l)
1944    IF (ifac>1) GO TO 10
1945
1946!    WRITE (nout,'(A,I4,A)') ' n =',n,' - Contains illegal factors'
1947    IF ( myid == 0 )  THEN
1948       PRINT*, '+++ temperton_fft:'
1949       PRINT*, '    number of gridpoints along x or/and y contain illegal ', &
1950                    'factors'
1951       PRINT*, '    only factors 8,6,5,4,3,2 are allowed'
1952    ENDIF
1953    CALL local_stop
1954
1955
1956    RETURN
1957
1958    ! Now reverse order of factors
195940  CONTINUE
1960    nfax = k
1961    ifax(1) = nfax
1962    DO i = 1, nfax
1963       ifax(nfax+2-i) = jfax(i)
1964    END DO
1965    ifax(10) = n
1966    RETURN
1967  END SUBROUTINE set99
1968
1969 END MODULE temperton_fft
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |