Home


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

source: palm/trunk/SOURCE/temperton_fft_mod.f90 @ 3045

Last change on this file since 3045 was 3045, checked in by Giersch, 6 years ago

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