Home


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

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

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago

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