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

Last change on this file since 1321 was 1321, checked in by raasch, 11 years ago

last commit documented

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