Home


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

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

Last change on this file since 1332 was 1323, checked in by raasch, 11 years ago

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