Home


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

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

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago

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