Home


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

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

Last change on this file since 3726 was 3725, checked in by raasch, 6 years ago

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