mo_fft.f90
MODULE mo_fft
! ---------------------------------------------------------------
!
5: ! module *mo_fft* - quantities needed for the fast *fourier transforms.
!
! ---------------------------------------------------------------
IMPLICIT NONE
10:
REAL, ALLOCATABLE :: trig(:) ! array used by *fft991*.
INTEGER :: nfax(10) ! array used by *fft991*.
15: ! nfft: maximum length of calls to *fft.
#ifdef _CRAY1
#if (_MAXVL == 64)
INTEGER, PARAMETER :: nfft = 64
20: #endif
#if (_MAXVL == 128)
INTEGER, PARAMETER :: nfft = 128
#endif
#else
25: #ifdef SX
INTEGER, PARAMETER :: nfft = 256
#else
#ifdef sun
INTEGER, PARAMETER :: nfft = 192
30: #else
INTEGER, PARAMETER :: nfft = 2048
#endif
#endif
#endif
35:
CONTAINS
SUBROUTINE inifft(nlon)
40:
! Description:
!
! Preset constants in *mo_fft*.
!
45: ! Method:
!
! *inifft* is called from *setdyn*.
!
! Externals:
50: ! *set99* set up trigonometric tables.
!
! Authors:
!
! M. Jarraud, ECMWF, December 1982, original source
55: ! L. Kornblueh, MPI, May 1998, f90 rewrite
! U. Schulzweida, MPI, May 1998, f90 rewrite
!
! for more details see file AUTHORS
!
60:
USE mo_doctor, ONLY: nout
USE mo_mpi, ONLY: p_pe, p_io
IMPLICIT NONE
65:
INTEGER :: nlon
! Executable statements
70:
IF (.NOT. ALLOCATED(trig)) ALLOCATE(trig(nlon))
! Define parameters depending on length of vectorregisters
75: IF (p_Pe == p_io) THEN
WRITE(nout,'(/,a,i5,/)') ' Maximum ffts to run in parallel:',nfft
END IF
!-- 1. Set up trigonometric tables
80:
CALL set99(trig,nfax,nlon)
END SUBROUTINE inifft
85: SUBROUTINE fft991cy(a,work,trigs,ifax,inc,jump,n,lot,isign)
! Description:
!
! Calls fortran-versions of fft's.
90: !
! Method:
!
! Subroutine 'fft991cy' - multiple fast real periodic transform
! supersedes previous routine 'fft991cy'.
95: !
! Real transform of length n performed by removing redundant
! operations from complex transform of length n.
!
! a is the array containing input & output data.
100: ! work is an area of size (n+1)*min(lot,nfft).
! trigs is a previously prepared list of trig function values.
! ifax is a previously prepared list of factors of n.
! inc is the increment within each data 'vector'
! (e.g. inc=1 for consecutively stored data).
105: ! jump is the increment between the start of each data vector.
! n is the length of the data vectors.
! lot is the number of data vectors.
! isign = +1 for transform from spectral to gridpoint
! = -1 for transform from gridpoint to spectral.
110: !
! ordering of coefficients:
! a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2)
! where b(0)=b(n/2)=0; (n+2) locations required.
!
115: ! ordering of data:
! x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required.
!
! Vectorization is achieved on cray by doing the transforms
! in parallel.
120: !
! n must be composed of factors 2,3 & 5 but does not have to be even.
!
! definition of transforms:
!
125: ! isign=+1: x(j)=sum(k=0,.,n-1)(c(k)*exp(2*i*j*k*pi/n))
! where c(k)=a(k)+i*b(k) and c(n-k)=a(k)-i*b(k)
! isign=-1: a(k)=(1/n)*sum(j=0,.,n-1)(x(j)*cos(2*j*k*pi/n))
! b(k)=-(1/n)*sum(j=0,.,n-1)(x(j)*sin(2*j*k*pi/n))
130:
! calls fortran-versions of fft's !!!
! dimension a(n),work(n),trigs(n),ifax(1)
USE mo_doctor, ONLY: nout
135:
IMPLICIT NONE
! Scalar arguments
INTEGER :: inc, isign, jump, lot, n
140:
! Array arguments
REAL :: a(*), trigs(*), work(*)
INTEGER :: ifax(*)
145: ! Local scalars:
INTEGER :: i, ia, ibase, ierr, ifac, igo, ii, istart, ix, iz, j, jbase, jj, &
& k, la, nb, nblox, nfax, nvex, nx
! Intrinsic functions
150: INTRINSIC MOD
! Executable statements
155: IF (ifax(10)/=n) CALL set99(trigs,ifax,n)
nfax = ifax(1)
nx = n + 1
IF (MOD(n,2)==1) nx = n
nblox = 1 + (lot-1)/nfft
160: nvex = lot - (nblox-1)*nfft
IF (isign==-1) GO TO 50
! isign=+1, spectral to gridpoint transform
165: istart = 1
DO nb = 1, nblox
ia = istart
i = istart
!DIR$ IVDEP
170: !CDIR NODEP
!OCL NOVREC
DO j = 1, nvex
a(i+inc) = 0.5*a(i)
i = i + jump
175: END DO
IF (MOD(n,2)==1) GO TO 10
i = istart + n*inc
DO j = 1, nvex
a(i) = 0.5*a(i)
180: i = i + jump
END DO
10 CONTINUE
ia = istart + inc
la = 1
185: igo = + 1
DO k = 1, nfax
ifac = ifax(k+1)
ierr = -1
190: IF (igo==-1) GO TO 20
CALL rpassm(a(ia),a(ia+la*inc),work(1),work(ifac*la+1),trigs,inc,1, &
& jump,nx,nvex,n,ifac,la,ierr)
GO TO 30
20 CONTINUE
195: CALL rpassm(work(1),work(la+1),a(ia),a(ia+ifac*la*inc),trigs,1,inc,nx, &
& jump,nvex,n,ifac,la,ierr)
30 CONTINUE
IF (ierr/=0) GO TO 100
la = ifac*la
200: igo = -igo
ia = istart
END DO
! If necessary, copy results back to a
205:
IF (MOD(nfax,2)==0) GO TO 40
ibase = 1
jbase = ia
DO jj = 1, nvex
210: i = ibase
j = jbase
DO ii = 1, n
a(j) = work(i)
i = i + 1
215: j = j + inc
END DO
ibase = ibase + nx
jbase = jbase + jump
END DO
220: 40 CONTINUE
! Fill in zeros at end
ix = istart + n*inc
225: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO j = 1, nvex
a(ix) = 0.0
230: a(ix+inc) = 0.0
ix = ix + jump
END DO
istart = istart + nvex*jump
235: nvex = nfft
END DO
RETURN
! isign=-1, gridpoint to spectral transform
240:
50 CONTINUE
istart = 1
DO nb = 1, nblox
ia = istart
245: la = n
igo = + 1
DO k = 1, nfax
ifac = ifax(nfax+2-k)
250: la = la/ifac
ierr = -1
IF (igo==-1) GO TO 60
CALL qpassm(a(ia),a(ia+ifac*la*inc),work(1),work(la+1),trigs,inc,1, &
& jump,nx,nvex,n,ifac,la,ierr)
255: GO TO 70
60 CONTINUE
CALL qpassm(work(1),work(ifac*la+1),a(ia),a(ia+la*inc),trigs,1,inc,nx, &
& jump,nvex,n,ifac,la,ierr)
70 CONTINUE
260: IF (ierr/=0) GO TO 100
igo = -igo
ia = istart + inc
END DO
265: ! If necessary, copy results back to a
IF (MOD(nfax,2)==0) GO TO 80
ibase = 1
jbase = ia
270: DO jj = 1, nvex
i = ibase
j = jbase
DO ii = 1, n
a(j) = work(i)
275: i = i + 1
j = j + inc
END DO
ibase = ibase + nx
jbase = jbase + jump
280: END DO
80 CONTINUE
! Shift a(0) & fill in zero imag parts
285: ix = istart
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO j = 1, nvex
290: a(ix) = a(ix+inc)
a(ix+inc) = 0.0
ix = ix + jump
END DO
IF (MOD(n,2)==1) GO TO 90
295: iz = istart + (n+1)*inc
DO j = 1, nvex
a(iz) = 0.0
iz = iz + jump
END DO
300: 90 CONTINUE
istart = istart + nvex*jump
nvex = nfft
END DO
305: RETURN
! Error messages
100 CONTINUE
310:
SELECT CASE (ierr)
CASE (:-1)
WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
CASE (0)
315: WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
CASE (1:)
WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
END SELECT
320: RETURN
END SUBROUTINE fft991cy
SUBROUTINE qpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
325: ! Description:
!
! Performs one pass through data as part of
! multiple real fft (fourier analysis) routine.
!
330: ! Method:
!
! a is first real input vector
! equivalence b(1) with a(ifac*la*inc1+1)
! c is first real output vector
335: ! equivalence d(1) with c(la*inc2+1)
! trigs is a precalculated list of sines & cosines
! inc1 is the addressing increment for a
! inc2 is the addressing increment for c
! inc3 is the increment between input vectors a
340: ! inc4 is the increment between output vectors c
! lot is the number of vectors
! n is the length of the vectors
! ifac is the current factor of n
! la = n/(product of factors used so far)
345: ! ierr is an error indicator:
! 0 - pass completed without error
! 1 - lot greater than nfft
! 2 - ifac not catered for
! 3 - ifac only catered for if la=n/ifac
350: !
IMPLICIT NONE
! Scalar arguments
355: INTEGER :: ierr, ifac, inc1, inc2, inc3, inc4, la, lot, n
! Array arguments
! REAL :: a(n),b(n),c(n),d(n),trigs(n)
REAL :: a(*), b(*), c(*), d(*), trigs(*)
360:
! Local scalars:
REAL :: a0, a1, a10, a11, a2, a20, a21, a3, a4, a5, a6, b0, b1, b10, b11, &
& b2, b20, b21, b3, b4, b5, b6, c1, c2, c3, c4, c5, qrt5, s1, s2, s3, s4, &
& s5, sin36, sin45, sin60, sin72, z, zqrt5, zsin36, zsin45, zsin60, &
365: & zsin72
INTEGER :: i, ia, ib, ibad, ibase, ic, id, ie, if, ig, igo, ih, iink, ijk, &
& ijump, j, ja, jb, jbase, jc, jd, je, jf, jink, k, kb, kc, kd, ke, kf, &
& kstop, l, m
370: ! Intrinsic functions
INTRINSIC REAL, SQRT
! Data statements
DATA sin36/0.587785252292473/, sin72/0.951056516295154/, &
375: & qrt5/0.559016994374947/, sin60/0.866025403784437/
! Executable statements
380: m = n/ifac
iink = la*inc1
jink = la*inc2
ijump = (ifac-1)*iink
kstop = (n-ifac)/(2*ifac)
385:
ibad = 1
IF (lot>nfft) GO TO 180
ibase = 0
jbase = 0
390: igo = ifac - 1
IF (igo==7) igo = 6
ibad = 2
IF (igo<1 .OR. igo>6) GO TO 180
GO TO (10,40,70,100,130,160) igo
395:
! Coding for factor 2
10 CONTINUE
ia = 1
400: ib = ia + iink
ja = 1
jb = ja + (2*m-la)*inc2
IF (la==m) GO TO 30
405:
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
410: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i) + a(ib+i)
c(jb+j) = a(ia+i) - a(ib+i)
415: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
420: END DO
ja = ja + jink
jink = 2*jink
jb = jb - jink
ibase = ibase + ijump
425: ijump = 2*ijump + iink
IF (ja==jb) GO TO 20
DO k = la, kstop, la
kb = k + k
c1 = trigs(kb+1)
430: s1 = trigs(kb+2)
jbase = 0
DO l = 1, la
i = ibase
j = jbase
435: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
440: c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
i = i + inc3
j = j + inc4
445: END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
ibase = ibase + ijump
450: ja = ja + jink
jb = jb - jink
END DO
IF (ja>jb) GO TO 170
20 CONTINUE
455: jbase = 0
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
460: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i)
d(ja+j) = -a(ib+i)
465: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
470: END DO
GO TO 170
30 CONTINUE
z = 1.0/REAL(n)
475: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
480: !OCL NOVREC
DO ijk = 1, lot
c(ja+j) = z*(a(ia+i)+a(ib+i))
c(jb+j) = z*(a(ia+i)-a(ib+i))
i = i + inc3
485: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
490: GO TO 170
! Coding for factor 3
40 CONTINUE
495: ia = 1
ib = ia + iink
ic = ib + iink
ja = 1
jb = ja + (2*m-la)*inc2
500: jc = jb
IF (la==m) GO TO 60
DO l = 1, la
505: i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
510: DO ijk = 1, lot
c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
c(jb+j) = a(ia+i) - 0.5*(a(ib+i)+a(ic+i))
d(jb+j) = sin60*(a(ic+i)-a(ib+i))
i = i + inc3
515: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
520: ja = ja + jink
jink = 2*jink
jb = jb + jink
jc = jc - jink
ibase = ibase + ijump
525: ijump = 2*ijump + iink
IF (ja==jc) GO TO 50
DO k = la, kstop, la
kb = k + k
kc = kb + kb
530: c1 = trigs(kb+1)
s1 = trigs(kb+2)
c2 = trigs(kc+1)
s2 = trigs(kc+2)
jbase = 0
535: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
540: !OCL NOVREC
DO ijk = 1, lot
a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
a2 = a(ia+i) - 0.5*a1
545: b2 = b(ia+i) - 0.5*b1
a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
c(ja+j) = a(ia+i) + a1
d(ja+j) = b(ia+i) + b1
550: c(jb+j) = a2 + b3
d(jb+j) = b2 - a3
c(jc+j) = a2 - b3
d(jc+j) = -(b2+a3)
i = i + inc3
555: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
560: ibase = ibase + ijump
ja = ja + jink
jb = jb + jink
jc = jc - jink
END DO
565: IF (ja>jc) GO TO 170
50 CONTINUE
jbase = 0
DO l = 1, la
i = ibase
570: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
575: c(ja+j) = a(ia+i) + 0.5*(a(ib+i)-a(ic+i))
d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
i = i + inc3
j = j + inc4
580: END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
585: GO TO 170
60 CONTINUE
z = 1.0/REAL(n)
zsin60 = z*sin60
DO l = 1, la
590: i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
595: DO ijk = 1, lot
c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
c(jb+j) = z*(a(ia+i)-0.5*(a(ib+i)+a(ic+i)))
d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
i = i + inc3
600: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
605: GO TO 170
! Coding for factor 4
70 CONTINUE
610: ia = 1
ib = ia + iink
ic = ib + iink
id = ic + iink
ja = 1
615: jb = ja + (2*m-la)*inc2
jc = jb + 2*m*inc2
jd = jb
IF (la==m) GO TO 90
620:
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
625: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
630: c(jb+j) = a(ia+i) - a(ic+i)
d(jb+j) = a(id+i) - a(ib+i)
i = i + inc3
j = j + inc4
END DO
635: ibase = ibase + inc1
jbase = jbase + inc2
END DO
ja = ja + jink
jink = 2*jink
640: jb = jb + jink
jc = jc - jink
jd = jd - jink
ibase = ibase + ijump
ijump = 2*ijump + iink
645: IF (jb==jc) GO TO 80
DO k = la, kstop, la
kb = k + k
kc = kb + kb
kd = kc + kb
650: c1 = trigs(kb+1)
s1 = trigs(kb+2)
c2 = trigs(kc+1)
s2 = trigs(kc+2)
c3 = trigs(kd+1)
655: s3 = trigs(kd+2)
jbase = 0
DO l = 1, la
i = ibase
j = jbase
660: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
665: a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
670: b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
c(ja+j) = a0 + a1
c(jc+j) = a0 - a1
d(ja+j) = b0 + b1
675: d(jc+j) = b1 - b0
c(jb+j) = a2 + b3
c(jd+j) = a2 - b3
d(jb+j) = b2 - a3
d(jd+j) = -(b2+a3)
680: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
685: END DO
ibase = ibase + ijump
ja = ja + jink
jb = jb + jink
jc = jc - jink
690: jd = jd - jink
END DO
IF (jb>jc) GO TO 170
80 CONTINUE
sin45 = SQRT(0.5)
695: jbase = 0
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
700: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
705: d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
i = i + inc3
j = j + inc4
END DO
710: ibase = ibase + inc1
jbase = jbase + inc2
END DO
GO TO 170
715: 90 CONTINUE
z = 1.0/REAL(n)
DO l = 1, la
i = ibase
j = jbase
720: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
725: c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
c(jb+j) = z*(a(ia+i)-a(ic+i))
d(jb+j) = z*(a(id+i)-a(ib+i))
i = i + inc3
j = j + inc4
730: END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
GO TO 170
735:
! Coding for factor 5
100 CONTINUE
ia = 1
740: ib = ia + iink
ic = ib + iink
id = ic + iink
ie = id + iink
ja = 1
745: jb = ja + (2*m-la)*inc2
jc = jb + 2*m*inc2
jd = jc
je = jb
750: IF (la==m) GO TO 120
DO l = 1, la
i = ibase
j = jbase
755: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
a1 = a(ib+i) + a(ie+i)
760: a3 = a(ib+i) - a(ie+i)
a2 = a(ic+i) + a(id+i)
a4 = a(ic+i) - a(id+i)
a5 = a(ia+i) - 0.25*(a1+a2)
a6 = qrt5*(a1-a2)
765: c(ja+j) = a(ia+i) + (a1+a2)
c(jb+j) = a5 + a6
c(jc+j) = a5 - a6
d(jb+j) = -sin72*a3 - sin36*a4
d(jc+j) = -sin36*a3 + sin72*a4
770: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
775: END DO
ja = ja + jink
jink = 2*jink
jb = jb + jink
jc = jc + jink
780: jd = jd - jink
je = je - jink
ibase = ibase + ijump
ijump = 2*ijump + iink
IF (jb==jd) GO TO 110
785: DO k = la, kstop, la
kb = k + k
kc = kb + kb
kd = kc + kb
ke = kd + kb
790: c1 = trigs(kb+1)
s1 = trigs(kb+2)
c2 = trigs(kc+1)
s2 = trigs(kc+2)
c3 = trigs(kd+1)
795: s3 = trigs(kd+2)
c4 = trigs(ke+1)
s4 = trigs(ke+2)
jbase = 0
DO l = 1, la
800: i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
805: DO ijk = 1, lot
a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
810: b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
a5 = a(ia+i) - 0.25*(a1+a2)
815: a6 = qrt5*(a1-a2)
b5 = b(ia+i) - 0.25*(b1+b2)
b6 = qrt5*(b1-b2)
a10 = a5 + a6
a20 = a5 - a6
820: b10 = b5 + b6
b20 = b5 - b6
a11 = sin72*b3 + sin36*b4
a21 = sin36*b3 - sin72*b4
b11 = sin72*a3 + sin36*a4
825: b21 = sin36*a3 - sin72*a4
c(ja+j) = a(ia+i) + (a1+a2)
c(jb+j) = a10 + a11
c(je+j) = a10 - a11
c(jc+j) = a20 + a21
830: c(jd+j) = a20 - a21
d(ja+j) = b(ia+i) + (b1+b2)
d(jb+j) = b10 - b11
d(je+j) = -(b10+b11)
d(jc+j) = b20 - b21
835: d(jd+j) = -(b20+b21)
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
840: jbase = jbase + inc2
END DO
ibase = ibase + ijump
ja = ja + jink
jb = jb + jink
845: jc = jc + jink
jd = jd - jink
je = je - jink
END DO
IF (jb>jd) GO TO 170
850: 110 CONTINUE
jbase = 0
DO l = 1, la
i = ibase
j = jbase
855: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
a1 = a(ib+i) + a(ie+i)
860: a3 = a(ib+i) - a(ie+i)
a2 = a(ic+i) + a(id+i)
a4 = a(ic+i) - a(id+i)
a5 = a(ia+i) + 0.25*(a3-a4)
a6 = qrt5*(a3+a4)
865: c(ja+j) = a5 + a6
c(jb+j) = a5 - a6
c(jc+j) = a(ia+i) - (a3-a4)
d(ja+j) = -sin36*a1 - sin72*a2
d(jb+j) = -sin72*a1 + sin36*a2
870: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
875: END DO
GO TO 170
120 CONTINUE
z = 1.0/REAL(n)
880: zqrt5 = z*qrt5
zsin36 = z*sin36
zsin72 = z*sin72
DO l = 1, la
i = ibase
885: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
890: a1 = a(ib+i) + a(ie+i)
a3 = a(ib+i) - a(ie+i)
a2 = a(ic+i) + a(id+i)
a4 = a(ic+i) - a(id+i)
a5 = z*(a(ia+i)-0.25*(a1+a2))
895: a6 = zqrt5*(a1-a2)
c(ja+j) = z*(a(ia+i)+(a1+a2))
c(jb+j) = a5 + a6
c(jc+j) = a5 - a6
d(jb+j) = -zsin72*a3 - zsin36*a4
900: d(jc+j) = -zsin36*a3 + zsin72*a4
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
905: jbase = jbase + inc2
END DO
GO TO 170
! Coding for factor 6
910:
130 CONTINUE
ia = 1
ib = ia + iink
ic = ib + iink
915: id = ic + iink
ie = id + iink
if = ie + iink
ja = 1
jb = ja + (2*m-la)*inc2
920: jc = jb + 2*m*inc2
jd = jc + 2*m*inc2
je = jc
jf = jb
925: IF (la==m) GO TO 150
DO l = 1, la
i = ibase
j = jbase
930: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
935: c(ja+j) = (a(ia+i)+a(id+i)) + a11
c(jc+j) = (a(ia+i)+a(id+i)-0.5*a11)
d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
c(jb+j) = (a(ia+i)-a(id+i)) - 0.5*a11
940: d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
c(jd+j) = (a(ia+i)-a(id+i)) + a11
i = i + inc3
j = j + inc4
END DO
945: ibase = ibase + inc1
jbase = jbase + inc2
END DO
ja = ja + jink
jink = 2*jink
950: jb = jb + jink
jc = jc + jink
jd = jd - jink
je = je - jink
jf = jf - jink
955: ibase = ibase + ijump
ijump = 2*ijump + iink
IF (jc==jd) GO TO 140
DO k = la, kstop, la
kb = k + k
960: kc = kb + kb
kd = kc + kb
ke = kd + kb
kf = ke + kb
c1 = trigs(kb+1)
965: s1 = trigs(kb+2)
c2 = trigs(kc+1)
s2 = trigs(kc+2)
c3 = trigs(kd+1)
s3 = trigs(kd+2)
970: c4 = trigs(ke+1)
s4 = trigs(ke+2)
c5 = trigs(kf+1)
s5 = trigs(kf+2)
jbase = 0
975: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
980: !OCL NOVREC
DO ijk = 1, lot
a1 = c1*a(ib+i) + s1*b(ib+i)
b1 = c1*b(ib+i) - s1*a(ib+i)
a2 = c2*a(ic+i) + s2*b(ic+i)
985: b2 = c2*b(ic+i) - s2*a(ic+i)
a3 = c3*a(id+i) + s3*b(id+i)
b3 = c3*b(id+i) - s3*a(id+i)
a4 = c4*a(ie+i) + s4*b(ie+i)
b4 = c4*b(ie+i) - s4*a(ie+i)
990: a5 = c5*a(if+i) + s5*b(if+i)
b5 = c5*b(if+i) - s5*a(if+i)
a11 = (a2+a5) + (a1+a4)
a20 = (a(ia+i)+a3) - 0.5*a11
a21 = sin60*((a2+a5)-(a1+a4))
995: b11 = (b2+b5) + (b1+b4)
b20 = (b(ia+i)+b3) - 0.5*b11
b21 = sin60*((b2+b5)-(b1+b4))
c(ja+j) = (a(ia+i)+a3) + a11
d(ja+j) = (b(ia+i)+b3) + b11
1000: c(jc+j) = a20 - b21
d(jc+j) = a21 + b20
c(je+j) = a20 + b21
d(je+j) = a21 - b20
a11 = (a2-a5) + (a4-a1)
1005: a20 = (a(ia+i)-a3) - 0.5*a11
a21 = sin60*((a4-a1)-(a2-a5))
b11 = (b5-b2) - (b4-b1)
b20 = (b3-b(ia+i)) - 0.5*b11
b21 = sin60*((b5-b2)+(b4-b1))
1010: c(jb+j) = a20 - b21
d(jb+j) = a21 - b20
c(jd+j) = a11 + (a(ia+i)-a3)
d(jd+j) = b11 + (b3-b(ia+i))
c(jf+j) = a20 + b21
1015: d(jf+j) = a21 + b20
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1020: jbase = jbase + inc2
END DO
ibase = ibase + ijump
ja = ja + jink
jb = jb + jink
1025: jc = jc + jink
jd = jd - jink
je = je - jink
jf = jf - jink
END DO
1030: IF (jc>jd) GO TO 170
140 CONTINUE
jbase = 0
DO l = 1, la
i = ibase
1035: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
1040: c(ja+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
d(ja+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
c(jc+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
1045: d(jc+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1050: jbase = jbase + inc2
END DO
GO TO 170
150 CONTINUE
1055: z = 1.0/REAL(n)
zsin60 = z*sin60
DO l = 1, la
i = ibase
j = jbase
1060: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
1065: c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
c(jc+j) = z*((a(ia+i)+a(id+i))-0.5*a11)
d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
c(jb+j) = z*((a(ia+i)-a(id+i))-0.5*a11)
1070: d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
i = i + inc3
j = j + inc4
END DO
1075: ibase = ibase + inc1
jbase = jbase + inc2
END DO
GO TO 170
1080: ! Coding for factor 8
160 CONTINUE
ibad = 3
IF (la/=m) GO TO 180
1085: ia = 1
ib = ia + iink
ic = ib + iink
id = ic + iink
ie = id + iink
1090: if = ie + iink
ig = if + iink
ih = ig + iink
ja = 1
jb = ja + la*inc2
1095: jc = jb + 2*m*inc2
jd = jc + 2*m*inc2
je = jd + 2*m*inc2
z = 1.0/REAL(n)
zsin45 = z*SQRT(0.5)
1100:
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
1105: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
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))))
1110: 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))))
c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+ &
1115: & i)-a(ib+i)))
c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+ &
& i)-a(ib+i)))
d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + &
& z*(a(ig+i)-a(ic+i))
1120: d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - &
& z*(a(ig+i)-a(ic+i))
i = i + inc3
j = j + inc4
END DO
1125: ibase = ibase + inc1
jbase = jbase + inc2
END DO
! Return
1130:
170 CONTINUE
ibad = 0
180 CONTINUE
ierr = ibad
1135: RETURN
END SUBROUTINE qpassm
SUBROUTINE rpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
! Dimension a(n),b(n),c(n),d(n),trigs(n)
1140:
IMPLICIT NONE
! Scalar arguments
INTEGER :: ierr, ifac, inc1, inc2, inc3, inc4, la, lot, n
1145:
! Array arguments
REAL :: a(*), b(*), c(*), d(*), trigs(*)
! Local scalars:
1150: REAL :: c1, c2, c3, c4, c5, qqrt5, qrt5, s1, s2, s3, s4, s5, sin36, sin45, &
& sin60, sin72, ssin36, ssin45, ssin60, ssin72
INTEGER :: i, ia, ib, ibad, ibase, ic, id, ie, if, igo, iink, ijk, j, ja, &
& jb, jbase, jc, jd, je, jf, jg, jh, jink, jump, k, kb, kc, kd, ke, kf, &
& kstop, l, m
1155:
! Local arrays:
REAL :: a10(nfft), a11(nfft), a20(nfft), a21(nfft), b10(nfft), b11(nfft), b20(nfft), &
& b21(nfft)
1160: ! Intrinsic functions
INTRINSIC SQRT
! Data statements
DATA sin36/0.587785252292473/, sin72/0.951056516295154/, &
1165: & qrt5/0.559016994374947/, sin60/0.866025403784437/
! Executable statements
1170: m = n/ifac
iink = la*inc1
jink = la*inc2
jump = (ifac-1)*jink
kstop = (n-ifac)/(2*ifac)
1175:
ibad = 1
IF (lot>nfft) GO TO 180
ibase = 0
jbase = 0
1180: igo = ifac - 1
IF (igo==7) igo = 6
ibad = 2
IF (igo<1 .OR. igo>6) GO TO 180
GO TO (10,40,70,100,130,160) igo
1185:
! Coding for factor 2
10 CONTINUE
ia = 1
1190: ib = ia + (2*m-la)*inc1
ja = 1
jb = ja + jink
IF (la==m) GO TO 30
1195:
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
1200: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i) + a(ib+i)
c(jb+j) = a(ia+i) - a(ib+i)
1205: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
1210: END DO
ia = ia + iink
iink = 2*iink
ib = ib - iink
ibase = 0
1215: jbase = jbase + jump
jump = 2*jump + jink
IF (ia==ib) GO TO 20
DO k = la, kstop, la
kb = k + k
1220: c1 = trigs(kb+1)
s1 = trigs(kb+2)
ibase = 0
DO l = 1, la
i = ibase
1225: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
1230: c(ja+j) = a(ia+i) + a(ib+i)
d(ja+j) = b(ia+i) - b(ib+i)
c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i))
d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i))
i = i + inc3
1235: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
1240: ia = ia + iink
ib = ib - iink
jbase = jbase + jump
END DO
IF (ia>ib) GO TO 170
1245: 20 CONTINUE
ibase = 0
DO l = 1, la
i = ibase
j = jbase
1250: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i)
1255: c(jb+j) = -b(ia+i)
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1260: jbase = jbase + inc2
END DO
GO TO 170
30 CONTINUE
1265: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
1270: !OCL NOVREC
DO ijk = 1, lot
c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
c(jb+j) = 2.0*(a(ia+i)-a(ib+i))
i = i + inc3
1275: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
1280: GO TO 170
! Coding for factor 3
40 CONTINUE
1285: ia = 1
ib = ia + (2*m-la)*inc1
ic = ib
ja = 1
jb = ja + jink
1290: jc = jb + jink
IF (la==m) GO TO 60
DO l = 1, la
1295: i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
1300: DO ijk = 1, lot
c(ja+j) = a(ia+i) + a(ib+i)
c(jb+j) = (a(ia+i)-0.5*a(ib+i)) - (sin60*(b(ib+i)))
c(jc+j) = (a(ia+i)-0.5*a(ib+i)) + (sin60*(b(ib+i)))
i = i + inc3
1305: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
1310: ia = ia + iink
iink = 2*iink
ib = ib + iink
ic = ic - iink
jbase = jbase + jump
1315: jump = 2*jump + jink
IF (ia==ic) GO TO 50
DO k = la, kstop, la
kb = k + k
kc = kb + kb
1320: c1 = trigs(kb+1)
s1 = trigs(kb+2)
c2 = trigs(kc+1)
s2 = trigs(kc+2)
ibase = 0
1325: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
1330: !OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
c(jb+j) = c1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
1335: & b(ic+i)))) - s1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
& a(ic+i))))
d(jb+j) = s1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
& b(ic+i)))) + c1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
& a(ic+i))))
1340: c(jc+j) = c2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
& b(ic+i)))) - s2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
& a(ic+i))))
d(jc+j) = s2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
& b(ic+i)))) + c2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
1345: & a(ic+i))))
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1350: jbase = jbase + inc2
END DO
ia = ia + iink
ib = ib + iink
ic = ic - iink
1355: jbase = jbase + jump
END DO
IF (ia>ic) GO TO 170
50 CONTINUE
ibase = 0
1360: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
1365: !OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i) + a(ib+i)
c(jb+j) = (0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
c(jc+j) = -(0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
1370: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
1375: END DO
GO TO 170
60 CONTINUE
ssin60 = 2.0*sin60
1380: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
1385: !OCL NOVREC
DO ijk = 1, lot
c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
c(jb+j) = (2.0*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
c(jc+j) = (2.0*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
1390: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
1395: END DO
GO TO 170
! Coding for factor 4
1400: 70 CONTINUE
ia = 1
ib = ia + (2*m-la)*inc1
ic = ib + 2*m*inc1
id = ib
1405: ja = 1
jb = ja + jink
jc = jb + jink
jd = jc + jink
1410: IF (la==m) GO TO 90
DO l = 1, la
i = ibase
j = jbase
1415: !DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
1420: c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
i = i + inc3
j = j + inc4
1425: END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
ia = ia + iink
1430: iink = 2*iink
ib = ib + iink
ic = ic - iink
id = id - iink
jbase = jbase + jump
1435: jump = 2*jump + jink
IF (ib==ic) GO TO 80
DO k = la, kstop, la
kb = k + k
kc = kb + kb
1440: kd = kc + kb
c1 = trigs(kb+1)
s1 = trigs(kb+2)
c2 = trigs(kc+1)
s2 = trigs(kc+2)
1445: c3 = trigs(kd+1)
s3 = trigs(kd+2)
ibase = 0
DO l = 1, la
i = ibase
1450: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
1455: c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+ &
& i)-b(ic+i))-(b(ib+i)-b(id+i)))
d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+ &
1460: & i)-b(ic+i))-(b(ib+i)-b(id+i)))
c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+ &
& i)+b(ic+i))+(a(ib+i)-a(id+i)))
d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+ &
& i)+b(ic+i))+(a(ib+i)-a(id+i)))
1465: c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+ &
& i)+b(ic+i))-(a(ib+i)-a(id+i)))
d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+ &
& i)+b(ic+i))-(a(ib+i)-a(id+i)))
i = i + inc3
1470: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
1475: ia = ia + iink
ib = ib + iink
ic = ic - iink
id = id - iink
jbase = jbase + jump
1480: END DO
IF (ib>ic) GO TO 170
80 CONTINUE
ibase = 0
sin45 = SQRT(0.5)
1485: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
1490: !OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i) + a(ib+i)
c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
c(jc+j) = b(ib+i) - b(ia+i)
1495: c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1500: jbase = jbase + inc2
END DO
GO TO 170
90 CONTINUE
1505: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
1510: !OCL NOVREC
DO ijk = 1, lot
c(ja+j) = 2.0*((a(ia+i)+a(ic+i))+a(ib+i))
c(jb+j) = 2.0*((a(ia+i)-a(ic+i))-b(ib+i))
c(jc+j) = 2.0*((a(ia+i)+a(ic+i))-a(ib+i))
1515: c(jd+j) = 2.0*((a(ia+i)-a(ic+i))+b(ib+i))
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1520: jbase = jbase + inc2
END DO
! Coding for factor 5
1525: GO TO 170
100 CONTINUE
ia = 1
ib = ia + (2*m-la)*inc1
ic = ib + 2*m*inc1
1530: id = ic
ie = ib
ja = 1
jb = ja + jink
jc = jb + jink
1535: jd = jc + jink
je = jd + jink
IF (la==m) GO TO 120
1540: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
1545: !OCL NOVREC
DO ijk = 1, lot
c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
c(jb+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
& (sin72*b(ib+i)+sin36*b(ic+i))
1550: c(jc+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
& (sin36*b(ib+i)-sin72*b(ic+i))
c(jd+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
& (sin36*b(ib+i)-sin72*b(ic+i))
c(je+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
1555: & (sin72*b(ib+i)+sin36*b(ic+i))
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1560: jbase = jbase + inc2
END DO
ia = ia + iink
iink = 2*iink
ib = ib + iink
1565: ic = ic + iink
id = id - iink
ie = ie - iink
jbase = jbase + jump
jump = 2*jump + jink
1570: IF (ib==id) GO TO 110
DO k = la, kstop, la
kb = k + k
kc = kb + kb
kd = kc + kb
1575: ke = kd + kb
c1 = trigs(kb+1)
s1 = trigs(kb+2)
c2 = trigs(kc+1)
s2 = trigs(kc+2)
1580: c3 = trigs(kd+1)
s3 = trigs(kd+2)
c4 = trigs(ke+1)
s4 = trigs(ke+2)
ibase = 0
1585: DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
!CDIR NODEP
1590: !OCL NOVREC
DO ijk = 1, lot
a10(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
& qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1595: a20(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
& qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
b10(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
& qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
b20(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
1600: & qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
1605:
c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
1610: c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
1615: d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
i = i + inc3
j = j + inc4
END DO
1620: ibase = ibase + inc1
jbase = jbase + inc2
END DO
ia = ia + iink
ib = ib + iink
1625: ic = ic + iink
id = id - iink
ie = ie - iink
jbase = jbase + jump
END DO
1630: IF (ib>id) GO TO 170
110 CONTINUE
ibase = 0
DO l = 1, la
i = ibase
1635: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
1640: c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
& (sin36*b(ia+i)+sin72*b(ib+i))
c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
& (sin36*b(ia+i)+sin72*b(ib+i))
1645: c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
& (sin72*b(ia+i)-sin36*b(ib+i))
c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
& (sin72*b(ia+i)-sin36*b(ib+i))
i = i + inc3
1650: j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
1655:
GO TO 170
120 CONTINUE
qqrt5 = 2.0*qrt5
ssin36 = 2.0*sin36
1660: ssin72 = 2.0*sin72
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
1665: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = 2.0*(a(ia+i)+(a(ib+i)+a(ic+i)))
c(jb+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
1670: & i))) - (ssin72*b(ib+i)+ssin36*b(ic+i))
c(jc+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
& i))) - (ssin36*b(ib+i)-ssin72*b(ic+i))
c(jd+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
& i))) + (ssin36*b(ib+i)-ssin72*b(ic+i))
1675: c(je+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
& i))) + (ssin72*b(ib+i)+ssin36*b(ic+i))
i = i + inc3
j = j + inc4
END DO
1680: ibase = ibase + inc1
jbase = jbase + inc2
END DO
GO TO 170
1685: ! Coding for factor 6
130 CONTINUE
ia = 1
ib = ia + (2*m-la)*inc1
1690: ic = ib + 2*m*inc1
id = ic + 2*m*inc1
ie = ic
if = ib
ja = 1
1695: jb = ja + jink
jc = jb + jink
jd = jc + jink
je = jd + jink
jf = je + jink
1700:
IF (la==m) GO TO 150
DO l = 1, la
i = ibase
1705: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
1710: c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
c(jb+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
& i)+b(ic+i)))
c(jf+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
1715: & i)+b(ic+i)))
c(jc+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
& i)-b(ic+i)))
c(je+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
& i)-b(ic+i)))
1720: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
1725: END DO
ia = ia + iink
iink = 2*iink
ib = ib + iink
ic = ic + iink
1730: id = id - iink
ie = ie - iink
if = if - iink
jbase = jbase + jump
jump = 2*jump + jink
1735: IF (ic==id) GO TO 140
DO k = la, kstop, la
kb = k + k
kc = kb + kb
kd = kc + kb
1740: ke = kd + kb
kf = ke + kb
c1 = trigs(kb+1)
s1 = trigs(kb+2)
c2 = trigs(kc+1)
1745: s2 = trigs(kc+2)
c3 = trigs(kd+1)
s3 = trigs(kd+2)
c4 = trigs(ke+1)
s4 = trigs(ke+2)
1750: c5 = trigs(kf+1)
s5 = trigs(kf+2)
ibase = 0
DO l = 1, la
i = ibase
1755: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
1760:
a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
a20(ijk) = (a(ia+i)+a(id+i)) - 0.5*a11(ijk)
a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
1765: b20(ijk) = (b(ia+i)-b(id+i)) - 0.5*b11(ijk)
b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
1770: c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
1775: a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
a20(ijk) = (a(ia+i)-a(id+i)) - 0.5*a11(ijk)
a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
b20(ijk) = (b(ia+i)+b(id+i)) + 0.5*b11(ijk)
1780: b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+ &
& i))-b11(ijk))
d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+ &
1785: & i))-b11(ijk))
c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
1790:
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1795: jbase = jbase + inc2
END DO
ia = ia + iink
ib = ib + iink
ic = ic + iink
1800: id = id - iink
ie = ie - iink
if = if - iink
jbase = jbase + jump
END DO
1805: IF (ic>id) GO TO 170
140 CONTINUE
ibase = 0
DO l = 1, la
i = ibase
1810: j = jbase
!DIR$ IVDEP
!CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
1815: c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
1820: c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
1825: jbase = jbase + inc2
END DO
GO TO 170
150 CONTINUE
1830: ssin60 = 2.0*sin60
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
1835: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = (2.0*(a(ia+i)+a(id+i))) + (2.0*(a(ib+i)+a(ic+i)))
c(jd+j) = (2.0*(a(ia+i)-a(id+i))) - (2.0*(a(ib+i)-a(ic+i)))
1840: c(jb+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
& i)+b(ic+i)))
c(jf+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
& i)+b(ic+i)))
c(jc+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
1845: & i)-b(ic+i)))
c(je+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
& i)-b(ic+i)))
i = i + inc3
j = j + inc4
1850: END DO
ibase = ibase + inc1
jbase = jbase + inc2
END DO
GO TO 170
1855:
! Coding for factor 8
160 CONTINUE
ibad = 3
1860: IF (la/=m) GO TO 180
ia = 1
ib = ia + la*inc1
ic = ib + 2*la*inc1
id = ic + 2*la*inc1
1865: ie = id + 2*la*inc1
ja = 1
jb = ja + jink
jc = jb + jink
jd = jc + jink
1870: je = jd + jink
jf = je + jink
jg = jf + jink
jh = jg + jink
ssin45 = SQRT(2.0)
1875:
DO l = 1, la
i = ibase
j = jbase
!DIR$ IVDEP
1880: !CDIR NODEP
!OCL NOVREC
DO ijk = 1, lot
c(ja+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
c(je+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
1885: c(jc+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
c(jg+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
c(jb+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
& i))-(b(ib+i)+b(id+i)))
c(jf+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
1890: & i))-(b(ib+i)+b(id+i)))
c(jd+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
& i))+(b(ib+i)+b(id+i)))
c(jh+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
& i))+(b(ib+i)+b(id+i)))
1895: i = i + inc3
j = j + inc4
END DO
ibase = ibase + inc1
jbase = jbase + inc2
1900: END DO
! Return
170 CONTINUE
1905: ibad = 0
180 CONTINUE
ierr = ibad
RETURN
END SUBROUTINE rpassm
1910:
SUBROUTINE set99(trigs,ifax,n)
! Description:
!
1915: ! Computes factors of n & trigonometric functins required by fft99 & fft991cy
!
! Method:
!
! Dimension trigs(n),ifax(1),jfax(10),lfax(7)
1920: !
! subroutine 'set99' - computes factors of n & trigonometric
! functins required by fft99 & fft991cy
1925: USE mo_doctor, ONLY: nout
IMPLICIT NONE
! Scalar arguments
1930: INTEGER :: n
! Array arguments
REAL :: trigs(*)
INTEGER :: ifax(*)
1935:
! Local scalars:
REAL :: angle, del
INTEGER :: i, ifac, ixxx, k, l, nfax, nhl, nil, nu
1940: ! Local arrays:
INTEGER :: jfax(10), lfax(7)
! Intrinsic functions
INTRINSIC ASIN, COS, MOD, REAL, SIN
1945:
! Data statements
DATA lfax/6, 8, 5, 4, 3, 2, 1/
1950: ! Executable statements
ixxx = 1
del = 4.0*ASIN(1.0)/REAL(n)
nil = 0
1955: nhl = (n/2) - 1
DO k = nil, nhl
angle = REAL(k)*del
trigs(2*k+1) = COS(angle)
trigs(2*k+2) = SIN(angle)
1960: END DO
! Find factors of n (8,6,5,4,3,2; only one 8 allowed)
! Look for sixes first, store factors in descending order
nu = n
1965: ifac = 6
k = 0
l = 1
10 CONTINUE
IF (MOD(nu,ifac)/=0) GO TO 30
1970: k = k + 1
jfax(k) = ifac
IF (ifac/=8) GO TO 20
IF (k==1) GO TO 20
jfax(1) = 8
1975: jfax(k) = 6
20 CONTINUE
nu = nu/ifac
IF (nu==1) GO TO 40
IF (ifac/=8) GO TO 10
1980: 30 CONTINUE
l = l + 1
ifac = lfax(l)
IF (ifac>1) GO TO 10
1985: WRITE (nout,'(A,I4,A)') ' n =',n,' - Contains illegal factors'
RETURN
! Now reverse order of factors
1990: 40 CONTINUE
nfax = k
ifax(1) = nfax
DO i = 1, nfax
ifax(nfax+2-i) = jfax(i)
1995: END DO
ifax(10) = n
RETURN
END SUBROUTINE set99
2000: END MODULE mo_fft
Info Section
uses: mo_doctor, mo_mpi
calls: qpassm, rpassm, set99
back to top
ECHAM 4 vf90 (C) 1998 Max-Planck-Institut für Meteorologie, Hamburg
Wed Nov 24 01:25:21 CST 1999
HTML derived from FORTRAN source by f2html.pl v0.3 (C) 1997,98 Beroud Jean-Marc.