Home


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

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

Last change on this file since 3352 was 3241, checked in by raasch, 6 years ago

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