V3FIT
cfft99.f
1  SUBROUTINE cfft99(a, work, trigs, ifax, inc, jump, n, lot, isign)
2  USE stel_kinds
3  IMPLICIT NONE
4 C-----------------------------------------------
5 C D u m m y A r g u m e n t s
6 C-----------------------------------------------
7  INTEGER inc, jump, n, lot, isign
8  INTEGER, DIMENSION(*) :: ifax
9  REAL(rprec) :: a(*), work(*), trigs(*)
10 C-----------------------------------------------
11 C L o c a l V a r i a b l e s
12 C-----------------------------------------------
13  INTEGER :: nn, ink, jum, nfax, jnk, jst, ibase, ilast, nh, l, i1,
14  1 i2, m, jbase, i, j, igo, la, k
15  REAL(rprec) :: hREAL, himag
16 C-----------------------------------------------
17 c
18 c SUBROUTINE 'cfft99' - multiple fast complex fourier transform
19 c
20 c a is the array containing input and output data
21 c work is an area of SIZE n*lot
22 c trigs is a previously prepared list of trig FUNCTION values
23 c ifax is a previously prepared list of factors of n
24 c inc is the increment within each data "vector"
25 c (e.g. inc=1 for consecutively stored data)
26 c jump is the increment between the start of each data vector
27 c n is the length of the data vectors
28 c lot is the number of data vectors
29 c isign = +1 for transform from spectral to gridpoint
30 c = -1 for transform from gridpoint to spectral
31 c
32 c
33 c vectorization is achieved on cray by DOing the transforms in
34 c parallel.
35 c
36 c
37  nn = n+n
38  ink=inc+inc
39  jum = jump+jump
40  nfax=ifax(1)
41  jnk = 2
42  jst = 2
43  IF (isign.ge.0) GOTO 30
44  jnk = -2
45  jst = nn-2
46  IF (mod(nfax,2).eq.1) GOTO 40
47  ibase = 1
48  ilast = (n-1)*ink
49  nh = n/2
50  DO 20 l=1,lot
51  i1 = ibase+ink
52  i2 = ibase+ilast
53 cdir$ ivdep
54  DO 10 m=1,nh
55 c swap REAL and imaginary portions
56  hreal = a(i1)
57  himag = a(i1+1)
58  a(i1) = a(i2)
59  a(i1+1) = a(i2+1)
60  a(i2) = hreal
61  a(i2+1) = himag
62  i1 = i1+ink
63  i2 = i2-ink
64  10 CONTINUE
65  ibase = ibase+jum
66  20 CONTINUE
67  GOTO 100
68 c
69  30 CONTINUE
70  IF (mod(nfax,2).eq.0) GOTO 100
71 c
72  40 CONTINUE
73 c
74 c during the transform process, nfax steps are taken, and the
75 c results are stored alternately in work and in a. IF nfax is
76 c odd, the input data are first moved to work so that the final
77 c result (after nfax steps) is stored in array a.
78 c
79  ibase=1
80  jbase=1
81  DO 60 l=1,lot
82 c move REAL and imaginary portions of element zero
83  work(jbase) = a(ibase)
84  work(jbase+1) = a(ibase+1)
85  i=ibase+ink
86  j=jbase+jst
87 cdir$ ivdep
88  DO 50 m=2,n
89 c move REAL and imaginary portions of other elements (possibly in
90 c reverse order, depending on jst and jnk)
91  work(j) = a(i)
92  work(j+1) = a(i+1)
93  i=i+ink
94  j=j+jnk
95  50 CONTINUE
96  ibase=ibase+jum
97  jbase=jbase+nn
98  60 CONTINUE
99 c
100  100 CONTINUE
101 c
102 c perform the transform passes, one pass for each factor. during
103 c each pass the data are moved from a to work or from work to a.
104 c
105 c for nfax even, the first pass moves from a to work
106  igo = 110
107 c for nfax odd, the first pass moves from work to a
108  IF (mod(nfax,2).eq.1) igo = 120
109  la=1
110  DO 140 k=1,nfax
111  IF (igo.eq.120) GOTO 120
112  110 CONTINUE
113  CALL vpassm(a(1),a(2),work(1),work(2),trigs,
114  * ink,2,jum,nn,lot,n,ifax(k+1),la)
115  igo=120
116  GOTO 130
117  120 CONTINUE
118  CALL vpassm(work(1),work(2),a(1),a(2),trigs,
119  * 2,ink,nn,jum,lot,n,ifax(k+1),la)
120  igo=110
121  130 CONTINUE
122  la=la*ifax(k+1)
123  140 CONTINUE
124 c
125 c at this point the final transform result is stored in a.
126 c
127  END SUBROUTINE cfft99