Home


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

source: palm/trunk/SOURCE/temperton_fft_mod.f90 @ 4180

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago

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