Home


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

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

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

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