Home


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

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

Last change on this file since 1342 was 1342, checked in by kanani, 10 years ago

REAL constants defined as wp-kind
  • Property svn:keywords set to Id
File size: 57.1 KB
Line 
1 MODULE temperton_fft
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! REAL constants defined as wp-kind
7!
8! Former revisions:
9! -----------------
10! $Id: temperton_fft.f90 1342 2014-03-26 17:04:47Z kanani $
11!
12! 1322 2014-03-20 16:38:49Z raasch
13! REAL constants defined as wp-kind
14!
15! 1320 2014-03-20 08:40:49Z raasch
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
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
33    USE kinds
34
35    IMPLICIT NONE
36
37    PRIVATE
38
39    PUBLIC set99, fft991cy
40
41
42    INTEGER(iwp)          ::  nfax(10)   !: array used by *fft991*.
43    REAL(wp), ALLOCATABLE ::  trig(:)    !: array used by *fft991*.
44
45!
46!-- nfft: maximum length of calls to *fft.
47#if defined( __nec )
48    INTEGER(iwp), PARAMETER ::  nfft = 256  !:
49#else
50    INTEGER(iwp), PARAMETER ::  nfft =  32  !:
51#endif
52
53    INTEGER(iwp), PARAMETER ::  nout =   6  !: standard output stream
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
106    USE kinds
107
108    IMPLICIT NONE
109
110    !  Scalar arguments
111    INTEGER(iwp) ::  inc   !:
112    INTEGER(iwp) ::  isign !:
113    INTEGER(iwp) ::  jump  !:
114    INTEGER(iwp) ::  lot   !:
115    INTEGER(iwp) ::  n     !:
116
117    !  Array arguments
118    REAL(wp)     ::  a(*)     !:
119    REAL(wp)     ::  trigs(*) !:
120    REAL(wp)     ::  work(*)  !:
121    INTEGER(iwp) ::  ifax(*)  !:
122
123    !  Local scalars:
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     !:
144
145    !  Intrinsic functions
146!    INTRINSIC MOD
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_wp*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_wp*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_wp
226          a(ix+inc) = 0.0_wp
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_wp
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_wp
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
348    USE kinds
349
350    IMPLICIT NONE 
351
352    !  Scalar arguments
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    !:
362
363    !  Array arguments
364    ! REAL :: a(n),b(n),c(n),d(n),trigs(n)
365    REAL(wp) ::  a(*)     !:
366    REAL(wp) ::  b(*)     !:
367    REAL(wp) ::  c(*)     !:
368    REAL(wp) ::  d(*)     !:
369    REAL(wp) ::  trigs(*) !:
370 
371    !  Local scalars:
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 !:
415
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
450    !  Intrinsic functions
451!    INTRINSIC REAL, SQRT
452
453    !  Data statements
454    DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
455         &      qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
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
554    z = 1.0_wp/REAL(n,KIND=wp)
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_wp*(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_wp*a1
625             b2 = b(ia+i) - 0.5_wp*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_wp*(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
667    z = 1.0_wp/REAL(n,KIND=wp)
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_wp*(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
774    sin45 = SQRT(0.5_wp)
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
796    z = 1.0_wp/REAL(n,KIND=wp)
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_wp*(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_wp*(a1+a2)
895             a6 = qrt5*(a1-a2)
896             b5 = b(ia+i) - 0.25_wp*(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_wp*(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
959    z = 1.0_wp/REAL(n,KIND=wp)
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_wp*(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_wp*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_wp*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_wp*a11
1074             a21 = sin60*((a2+a5)-(a1+a4))
1075             b11 = (b2+b5) + (b1+b4)
1076             b20 = (b(ia+i)+b3) - 0.5_wp*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_wp*a11
1086             a21 = sin60*((a4-a1)-(a2-a5))
1087             b11 = (b5-b2) - (b4-b1)
1088             b20 = (b3-b(ia+i)) - 0.5_wp*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_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
1121          d(ja+j) = -(a(id+i)+0.5_wp*(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_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
1125          d(jc+j) = -(a(id+i)+0.5_wp*(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
1135    z = 1.0_wp/REAL(n,KIND=wp)
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_wp*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_wp*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
1178    z = 1.0_wp/REAL(n,KIND=wp)
1179    zsin45 = z*SQRT(0.5_wp)
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
1221    USE kinds
1222
1223    IMPLICIT NONE
1224
1225    !  Scalar arguments
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    !:
1235
1236    !  Array arguments
1237    REAL(wp) ::  a(*)     !:
1238    REAL(wp) ::  b(*)     !:
1239    REAL(wp) ::  c(*)     !:
1240    REAL(wp) ::  d(*)     !:
1241    REAL(wp) ::  trigs(*) !:
1242
1243    !  Local scalars:
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 !:
1264
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
1299    !  Local arrays:
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) !:
1308
1309    !  Intrinsic functions
1310!    INTRINSIC SQRT
1311
1312    !  Data statements
1313    DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
1314         &      qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
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_wp*(a(ia+i)+a(ib+i))
1422          c(jb+j) = 2.0_wp*(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_wp*a(ib+i)) - (sin60*(b(ib+i)))
1452          c(jc+j) = (a(ia+i)-0.5_wp*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_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
1484                  &            b(ic+i))))                                      &
1485                  &    - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
1486                  &            a(ic+i))))
1487             d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
1488                  &            b(ic+i))))                                      &
1489                  &    + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
1490                  &            a(ic+i))))
1491             c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
1492                  &            b(ic+i))))                                      &
1493                  &    - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
1494                  &            a(ic+i))))
1495             d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
1496                  &            b(ic+i))))                                      &
1497                  &    + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
1498                  &            a(ic+i))))
1499             i = i + inc3
1500             j = j + inc4
1501          END DO
1502          ibase = ibase + inc1
1503          jbase = jbase + inc2
1504       END DO
1505       ia = ia + iink
1506       ib = ib + iink
1507       ic = ic - iink
1508       jbase = jbase + jump
1509    END DO
1510    IF (ia>ic) GO TO 170
151150  CONTINUE
1512    ibase = 0
1513    DO l = 1, la
1514       i = ibase
1515       j = jbase
1516!DIR$ IVDEP
1517!CDIR NODEP
1518!OCL NOVREC
1519       DO ijk = 1, lot
1520          c(ja+j) = a(ia+i) + a(ib+i)
1521          c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
1522          c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
1523          i = i + inc3
1524          j = j + inc4
1525       END DO
1526       ibase = ibase + inc1
1527       jbase = jbase + inc2
1528    END DO
1529
1530    GO TO 170
153160  CONTINUE
1532    ssin60 = 2.0_wp*sin60
1533    DO l = 1, la
1534       i = ibase
1535       j = jbase
1536!DIR$ IVDEP
1537!CDIR NODEP
1538!OCL NOVREC
1539       DO ijk = 1, lot
1540          c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
1541          c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
1542          c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
1543          i = i + inc3
1544          j = j + inc4
1545       END DO
1546       ibase = ibase + inc1
1547       jbase = jbase + inc2
1548    END DO
1549    GO TO 170
1550
1551    ! Coding for factor 4
1552
155370  CONTINUE
1554    ia = 1
1555    ib = ia + (2*m-la)*inc1
1556    ic = ib + 2*m*inc1
1557    id = ib
1558    ja = 1
1559    jb = ja + jink
1560    jc = jb + jink
1561    jd = jc + jink
1562
1563    IF (la==m) GO TO 90
1564
1565    DO l = 1, la
1566       i = ibase
1567       j = jbase
1568!DIR$ IVDEP
1569!CDIR NODEP
1570!OCL NOVREC
1571       DO ijk = 1, lot
1572          c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
1573          c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
1574          c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
1575          c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
1576          i = i + inc3
1577          j = j + inc4
1578       END DO
1579       ibase = ibase + inc1
1580       jbase = jbase + inc2
1581    END DO
1582    ia = ia + iink
1583    iink = 2*iink
1584    ib = ib + iink
1585    ic = ic - iink
1586    id = id - iink
1587    jbase = jbase + jump
1588    jump = 2*jump + jink
1589    IF (ib==ic) GO TO 80
1590    DO k = la, kstop, la
1591       kb = k + k
1592       kc = kb + kb
1593       kd = kc + kb
1594       c1 = trigs(kb+1)
1595       s1 = trigs(kb+2)
1596       c2 = trigs(kc+1)
1597       s2 = trigs(kc+2)
1598       c3 = trigs(kd+1)
1599       s3 = trigs(kd+2)
1600       ibase = 0
1601       DO l = 1, la
1602          i = ibase
1603          j = jbase
1604!DIR$ IVDEP
1605!CDIR NODEP
1606!OCL NOVREC
1607          DO ijk = 1, lot
1608             c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
1609             d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
1610             c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+ &
1611                  &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
1612             d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+ &
1613                  &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
1614             c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+ &
1615                  &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
1616             d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+ &
1617                  &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
1618             c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+ &
1619                  &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
1620             d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+ &
1621                  &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
1622             i = i + inc3
1623             j = j + inc4
1624          END DO
1625          ibase = ibase + inc1
1626          jbase = jbase + inc2
1627       END DO
1628       ia = ia + iink
1629       ib = ib + iink
1630       ic = ic - iink
1631       id = id - iink
1632       jbase = jbase + jump
1633    END DO
1634    IF (ib>ic) GO TO 170
163580  CONTINUE
1636    ibase = 0
1637    sin45 = SQRT(0.5_wp)
1638    DO l = 1, la
1639       i = ibase
1640       j = jbase
1641!DIR$ IVDEP
1642!CDIR NODEP
1643!OCL NOVREC
1644       DO ijk = 1, lot
1645          c(ja+j) = a(ia+i) + a(ib+i)
1646          c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
1647          c(jc+j) = b(ib+i) - b(ia+i)
1648          c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
1649          i = i + inc3
1650          j = j + inc4
1651       END DO
1652       ibase = ibase + inc1
1653       jbase = jbase + inc2
1654    END DO
1655
1656    GO TO 170
165790  CONTINUE
1658    DO l = 1, la
1659       i = ibase
1660       j = jbase
1661!DIR$ IVDEP
1662!CDIR NODEP
1663!OCL NOVREC
1664       DO ijk = 1, lot
1665          c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
1666          c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
1667          c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
1668          c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
1669          i = i + inc3
1670          j = j + inc4
1671       END DO
1672       ibase = ibase + inc1
1673       jbase = jbase + inc2
1674    END DO
1675
1676    ! Coding for factor 5
1677
1678    GO TO 170
1679100 CONTINUE
1680    ia = 1
1681    ib = ia + (2*m-la)*inc1
1682    ic = ib + 2*m*inc1
1683    id = ic
1684    ie = ib
1685    ja = 1
1686    jb = ja + jink
1687    jc = jb + jink
1688    jd = jc + jink
1689    je = jd + jink
1690
1691    IF (la==m) GO TO 120
1692
1693    DO l = 1, la
1694       i = ibase
1695       j = jbase
1696!DIR$ IVDEP
1697!CDIR NODEP
1698!OCL NOVREC
1699       DO ijk = 1, lot
1700          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
1701          c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
1702               &          (sin72*b(ib+i)+sin36*b(ic+i))
1703          c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
1704               &          (sin36*b(ib+i)-sin72*b(ic+i))
1705          c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
1706               &          (sin36*b(ib+i)-sin72*b(ic+i))
1707          c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
1708               &          (sin72*b(ib+i)+sin36*b(ic+i))
1709          i = i + inc3
1710          j = j + inc4
1711       END DO
1712       ibase = ibase + inc1
1713       jbase = jbase + inc2
1714    END DO
1715    ia = ia + iink
1716    iink = 2*iink
1717    ib = ib + iink
1718    ic = ic + iink
1719    id = id - iink
1720    ie = ie - iink
1721    jbase = jbase + jump
1722    jump = 2*jump + jink
1723    IF (ib==id) GO TO 110
1724    DO k = la, kstop, la
1725       kb = k + k
1726       kc = kb + kb
1727       kd = kc + kb
1728       ke = kd + kb
1729       c1 = trigs(kb+1)
1730       s1 = trigs(kb+2)
1731       c2 = trigs(kc+1)
1732       s2 = trigs(kc+2)
1733       c3 = trigs(kd+1)
1734       s3 = trigs(kd+2)
1735       c4 = trigs(ke+1)
1736       s4 = trigs(ke+2)
1737       ibase = 0
1738       DO l = 1, la
1739          i = ibase
1740          j = jbase
1741!DIR$ IVDEP
1742!CDIR NODEP
1743!OCL NOVREC
1744          DO ijk = 1, lot
1745
1746             a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
1747                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1748             a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
1749                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1750             b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
1751                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
1752             b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
1753                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
1754             a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
1755             a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
1756             b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
1757             b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
1758
1759             c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
1760             d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
1761             c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
1762             d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
1763             c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
1764             d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
1765             c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
1766             d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
1767             c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
1768             d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
1769
1770             i = i + inc3
1771             j = j + inc4
1772          END DO
1773          ibase = ibase + inc1
1774          jbase = jbase + inc2
1775       END DO
1776       ia = ia + iink
1777       ib = ib + iink
1778       ic = ic + iink
1779       id = id - iink
1780       ie = ie - iink
1781       jbase = jbase + jump
1782    END DO
1783    IF (ib>id) GO TO 170
1784110 CONTINUE
1785    ibase = 0
1786    DO l = 1, la
1787       i = ibase
1788       j = jbase
1789!DIR$ IVDEP
1790!CDIR NODEP
1791!OCL NOVREC
1792       DO ijk = 1, lot
1793          c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
1794          c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1795               &          (sin36*b(ia+i)+sin72*b(ib+i))
1796          c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1797               &          (sin36*b(ia+i)+sin72*b(ib+i))
1798          c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1799               &          (sin72*b(ia+i)-sin36*b(ib+i))
1800          c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1801               &          (sin72*b(ia+i)-sin36*b(ib+i))
1802          i = i + inc3
1803          j = j + inc4
1804       END DO
1805       ibase = ibase + inc1
1806       jbase = jbase + inc2
1807    END DO
1808
1809    GO TO 170
1810120 CONTINUE
1811    qqrt5 = 2.0_wp*qrt5
1812    ssin36 = 2.0_wp*sin36
1813    ssin72 = 2.0_wp*sin72
1814    DO l = 1, la
1815       i = ibase
1816       j = jbase
1817!DIR$ IVDEP
1818!CDIR NODEP
1819!OCL NOVREC
1820       DO ijk = 1, lot
1821          c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
1822          c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
1823               &          i))) - (ssin72*b(ib+i)+ssin36*b(ic+i))
1824          c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
1825               &          i))) - (ssin36*b(ib+i)-ssin72*b(ic+i))
1826          c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
1827               &          i))) + (ssin36*b(ib+i)-ssin72*b(ic+i))
1828          c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
1829               &          i))) + (ssin72*b(ib+i)+ssin36*b(ic+i))
1830          i = i + inc3
1831          j = j + inc4
1832       END DO
1833       ibase = ibase + inc1
1834       jbase = jbase + inc2
1835    END DO
1836    GO TO 170
1837
1838    ! Coding for factor 6
1839
1840130 CONTINUE
1841    ia = 1
1842    ib = ia + (2*m-la)*inc1
1843    ic = ib + 2*m*inc1
1844    id = ic + 2*m*inc1
1845    ie = ic
1846    if = ib
1847    ja = 1
1848    jb = ja + jink
1849    jc = jb + jink
1850    jd = jc + jink
1851    je = jd + jink
1852    jf = je + jink
1853
1854    IF (la==m) GO TO 150
1855
1856    DO l = 1, la
1857       i = ibase
1858       j = jbase
1859!DIR$ IVDEP
1860!CDIR NODEP
1861!OCL NOVREC
1862       DO ijk = 1, lot
1863          c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
1864          c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
1865          c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
1866               &          i)+b(ic+i)))
1867          c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
1868               &          i)+b(ic+i)))
1869          c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
1870               &          i)-b(ic+i)))
1871          c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
1872               &          i)-b(ic+i)))
1873          i = i + inc3
1874          j = j + inc4
1875       END DO
1876       ibase = ibase + inc1
1877       jbase = jbase + inc2
1878    END DO
1879    ia = ia + iink
1880    iink = 2*iink
1881    ib = ib + iink
1882    ic = ic + iink
1883    id = id - iink
1884    ie = ie - iink
1885    if = if - iink
1886    jbase = jbase + jump
1887    jump = 2*jump + jink
1888    IF (ic==id) GO TO 140
1889    DO k = la, kstop, la
1890       kb = k + k
1891       kc = kb + kb
1892       kd = kc + kb
1893       ke = kd + kb
1894       kf = ke + kb
1895       c1 = trigs(kb+1)
1896       s1 = trigs(kb+2)
1897       c2 = trigs(kc+1)
1898       s2 = trigs(kc+2)
1899       c3 = trigs(kd+1)
1900       s3 = trigs(kd+2)
1901       c4 = trigs(ke+1)
1902       s4 = trigs(ke+2)
1903       c5 = trigs(kf+1)
1904       s5 = trigs(kf+2)
1905       ibase = 0
1906       DO l = 1, la
1907          i = ibase
1908          j = jbase
1909!DIR$ IVDEP
1910!CDIR NODEP
1911!OCL NOVREC
1912          DO ijk = 1, lot
1913
1914             a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
1915             a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
1916             a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
1917             b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
1918             b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
1919             b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
1920
1921             c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
1922             d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
1923             c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
1924             d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
1925             c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
1926             d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
1927
1928             a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
1929             b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
1930             a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
1931             a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
1932             b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
1933             b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
1934
1935             c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+ &
1936                  &            i))-b11(ijk))
1937             d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+ &
1938                  &            i))-b11(ijk))
1939             c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
1940             d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
1941             c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
1942             d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
1943
1944             i = i + inc3
1945             j = j + inc4
1946          END DO
1947          ibase = ibase + inc1
1948          jbase = jbase + inc2
1949       END DO
1950       ia = ia + iink
1951       ib = ib + iink
1952       ic = ic + iink
1953       id = id - iink
1954       ie = ie - iink
1955       if = if - iink
1956       jbase = jbase + jump
1957    END DO
1958    IF (ic>id) GO TO 170
1959140 CONTINUE
1960    ibase = 0
1961    DO l = 1, la
1962       i = ibase
1963       j = jbase
1964!DIR$ IVDEP
1965!CDIR NODEP
1966!OCL NOVREC
1967       DO ijk = 1, lot
1968          c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
1969          c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
1970          c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
1971          c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
1972          c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
1973          c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
1974          i = i + inc3
1975          j = j + inc4
1976       END DO
1977       ibase = ibase + inc1
1978       jbase = jbase + inc2
1979    END DO
1980
1981    GO TO 170
1982150 CONTINUE
1983    ssin60 = 2.0_wp*sin60
1984    DO l = 1, la
1985       i = ibase
1986       j = jbase
1987!DIR$ IVDEP
1988!CDIR NODEP
1989!OCL NOVREC
1990       DO ijk = 1, lot
1991          c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
1992          c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
1993          c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
1994               &          i)+b(ic+i)))
1995          c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
1996               &          i)+b(ic+i)))
1997          c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
1998               &          i)-b(ic+i)))
1999          c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
2000               &          i)-b(ic+i)))
2001          i = i + inc3
2002          j = j + inc4
2003       END DO
2004       ibase = ibase + inc1
2005       jbase = jbase + inc2
2006    END DO
2007    GO TO 170
2008
2009    ! Coding for factor 8
2010
2011160 CONTINUE
2012    ibad = 3
2013    IF (la/=m) GO TO 180
2014    ia = 1
2015    ib = ia + la*inc1
2016    ic = ib + 2*la*inc1
2017    id = ic + 2*la*inc1
2018    ie = id + 2*la*inc1
2019    ja = 1
2020    jb = ja + jink
2021    jc = jb + jink
2022    jd = jc + jink
2023    je = jd + jink
2024    jf = je + jink
2025    jg = jf + jink
2026    jh = jg + jink
2027    ssin45 = SQRT(2.0_wp)
2028
2029    DO l = 1, la
2030       i = ibase
2031       j = jbase
2032!DIR$ IVDEP
2033!CDIR NODEP
2034!OCL NOVREC
2035       DO ijk = 1, lot
2036          c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
2037          c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
2038          c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
2039          c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
2040          c(jb+j) = 2.0_wp*((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(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
2043               &              i))-(b(ib+i)+b(id+i)))
2044          c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
2045               &              i))+(b(ib+i)+b(id+i)))
2046          c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
2047               &              i))+(b(ib+i)+b(id+i)))
2048          i = i + inc3
2049          j = j + inc4
2050       END DO
2051       ibase = ibase + inc1
2052       jbase = jbase + inc2
2053    END DO
2054
2055    ! Return
2056
2057170 CONTINUE
2058    ibad = 0
2059180 CONTINUE
2060    ierr = ibad
2061    RETURN
2062  END SUBROUTINE rpassm
2063
2064  SUBROUTINE set99(trigs,ifax,n)
2065
2066    ! Description:
2067    !
2068    ! Computes factors of n & trigonometric functins required by fft99 & fft991cy
2069    !
2070    ! Method:
2071    !
2072    ! Dimension trigs(n),ifax(1),jfax(10),lfax(7)
2073    !
2074    ! subroutine 'set99' - computes factors of n & trigonometric
2075    ! functins required by fft99 & fft991cy
2076
2077
2078    USE control_parameters,                                                    &
2079        ONLY:  message_string
2080
2081    USE kinds
2082
2083    IMPLICIT NONE
2084
2085    !  Scalar arguments
2086    INTEGER(iwp) ::  n !:
2087
2088    !  Array arguments
2089    INTEGER(iwp) ::  ifax(*)  !:
2090    REAL(wp)     ::  trigs(*) !:
2091
2092
2093    !  Local scalars:
2094    REAL(wp) ::  angle    !:
2095    REAL(wp) ::  del      !:
2096    INTEGER(iwp) ::  i    !:
2097    INTEGER(iwp) ::  ifac !:
2098    INTEGER(iwp) ::  ixxx !:
2099    INTEGER(iwp) ::  k    !:
2100    INTEGER(iwp) ::  l    !:
2101    INTEGER(iwp) ::  nfax !:
2102    INTEGER(iwp) ::  nhl  !:
2103    INTEGER(iwp) ::  nil  !:
2104    INTEGER(iwp) ::  nu   !:
2105
2106    !  Local arrays:
2107    INTEGER(iwp) ::  jfax(10) !:
2108    INTEGER(iwp) ::  lfax(7)  !:
2109
2110    !  Intrinsic functions
2111!    INTRINSIC ASIN, COS, MOD, REAL, SIN
2112
2113    !  Data statements
2114    DATA lfax/6, 8, 5, 4, 3, 2, 1/
2115
2116
2117    !  Executable statements
2118    ixxx = 1
2119
2120    del = 4.0_wp*ASIN(1.0_wp)/REAL(n,KIND=wp)
2121    nil = 0
2122    nhl = (n/2) - 1
2123    DO k = nil, nhl
2124       angle = REAL(k,KIND=wp)*del
2125       trigs(2*k+1) = COS(angle)
2126       trigs(2*k+2) = SIN(angle)
2127    END DO
2128
2129    ! Find factors of n (8,6,5,4,3,2; only one 8 allowed)
2130    ! Look for sixes first, store factors in descending order
2131    nu = n
2132    ifac = 6
2133    k = 0
2134    l = 1
213510  CONTINUE
2136    IF (MOD(nu,ifac)/=0) GO TO 30
2137    k = k + 1
2138    jfax(k) = ifac
2139    IF (ifac/=8) GO TO 20
2140    IF (k==1) GO TO 20
2141    jfax(1) = 8
2142    jfax(k) = 6
214320  CONTINUE
2144    nu = nu/ifac
2145    IF (nu==1) GO TO 40
2146    IF (ifac/=8) GO TO 10
214730  CONTINUE
2148    l = l + 1
2149    ifac = lfax(l)
2150    IF (ifac>1) GO TO 10
2151
2152!    WRITE (nout,'(A,I4,A)') ' n =',n,' - Contains illegal factors'
2153    message_string = 'number of gridpoints along x or/and y ' // &
2154                     'contain illegal  factors' //               &
2155                     '&only factors 8,6,5,4,3,2 are allowed' 
2156    CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
2157
2158    RETURN
2159
2160    ! Now reverse order of factors
216140  CONTINUE
2162    nfax = k
2163    ifax(1) = nfax
2164    DO i = 1, nfax
2165       ifax(nfax+2-i) = jfax(i)
2166    END DO
2167    ifax(10) = n
2168    RETURN
2169  END SUBROUTINE set99
2170
2171 END MODULE temperton_fft
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |