Home


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

source: palm/trunk/SOURCE/temperton_fft.f90 @ 1682

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago

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