Home


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

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

Last change on this file since 357 was 258, checked in by heinze, 16 years ago

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