Home


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

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

Last change on this file since 1526 was 1343, checked in by kanani, 11 years ago

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