Home


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

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

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