Home


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

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

Last change on this file since 2284 was 1851, checked in by maronga, 9 years ago

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