Home


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

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

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

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