Home


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

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

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

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