V3FIT
fft991.f
1  SUBROUTINE fft991(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(13) :: ifax
9  REAL(rprec), DIMENSION(*) :: 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::nfax,nx,nh,ink,igo,ibase,jbase,l,i,j,m,ia,la,k,ib
14 C-----------------------------------------------
15 c
16 c SUBROUTINE 'fft991' - multiple REAL/half-complex periodic
17 c fast fourier transform
18 c
19 c same as fft99 except that ordering of data corresponds to
20 c that in mrfft2
21 c
22 c PROCEDURE used to convert to half-length complex transform
23 c is given by cooley, lewis and welch (j. sound vib., vol. 12
24 c (1970), 315-337)
25 c
26 c a is the array containing input and output data
27 c work is an area of SIZE (n+1)*lot
28 c trigs is a previously prepared list of trig FUNCTION values
29 c ifax is a previously prepared list of factors of n/2
30 c inc is the increment within each data "vector"
31 c (e.g. inc=1 for consecutively stored data)
32 c jump is the increment between the start of each data vector
33 c n is the length of the data vectors
34 c lot is the number of data vectors
35 c isign = +1 for transform from spectral to gridpoint
36 c = -1 for transform from gridpoint to spectral
37 c
38 c ordering of coefficients:
39 c a(0),b(0),a(1),b(1),a(2),b(2),...,a(n/2),b(n/2)
40 c WHERE b(0)=b(n/2)=0; (n+2) locations required
41 c
42 c ordering of data:
43 c x(0),x(1),x(2),...,x(n-1)
44 c
45 c vectorization is achieved on cray by DOing the transforms in
46 c parallel
47 c
48 c *** n.b. n is assumed to be an even number
49 c
50 c definition of transforms:
51 c -------------------------
52 c
53 c isign=+1: x(j)=SUM(k=0,...,n-1)(c(k)*EXP(2*i*j*k*pi/n))
54 c WHERE c(k)=a(k)+i*b(k) and c(n-k)=a(k)-i*b(k)
55 c
56 c isign=-1: a(k)=(1/n)*SUM(j=0,...,n-1)(x(j)*COS(2*j*k*pi/n))
57 c b(k)=-(1/n)*SUM(j=0,...,n-1)(x(j)*SIN(2*j*k*pi/n))
58 c
59  nfax=ifax(1)
60  nx=n+1
61  nh=n/2
62  ink=inc+inc
63  IF (isign.eq.+1) GOTO 30
64 c
65 c IF necessary, transfer data to work area
66  igo=50
67  IF (mod(nfax,2).eq.1) GOTO 40
68  ibase=1
69  jbase=1
70  DO 20 l=1,lot
71  i=ibase
72  j=jbase
73 cdir$ ivdep
74  DO 10 m=1,n
75  work(j)=a(i)
76  i=i+inc
77  j=j+1
78  10 CONTINUE
79  ibase=ibase+jump
80  jbase=jbase+nx
81  20 CONTINUE
82 c
83  igo=60
84  GOTO 40
85 c
86 c preprocessing (isign=+1)
87 c ------------------------
88 c
89  30 CONTINUE
90  CALL fft99a(a,work,trigs,inc,jump,n,lot)
91  igo=60
92 c
93 c complex transform
94 c -----------------
95 c
96  40 CONTINUE
97  ia=1
98  la=1
99  DO 80 k=1,nfax
100  IF (igo.eq.60) GOTO 60
101  50 CONTINUE
102  CALL vpassm(a(ia),a(ia+inc),work(1),work(2),trigs,
103  * ink,2,jump,nx,lot,nh,ifax(k+1),la)
104  igo=60
105  GOTO 70
106  60 CONTINUE
107  CALL vpassm(work(1),work(2),a(ia),a(ia+inc),trigs,
108  * 2,ink,nx,jump,lot,nh,ifax(k+1),la)
109  igo=50
110  70 CONTINUE
111  la=la*ifax(k+1)
112  80 CONTINUE
113 c
114  IF (isign.eq.-1) GOTO 130
115 c
116 c IF necessary, transfer data from work area
117  IF (mod(nfax,2).eq.1) GOTO 110
118  ibase=1
119  jbase=1
120  DO 100 l=1,lot
121  i=ibase
122  j=jbase
123 cdir$ ivdep
124  DO 90 m=1,n
125  a(j)=work(i)
126  i=i+1
127  j=j+inc
128  90 CONTINUE
129  ibase=ibase+nx
130  jbase=jbase+jump
131  100 CONTINUE
132 c
133 c fill in zeros at END
134  110 CONTINUE
135  ib=n*inc+1
136 cdir$ ivdep
137  DO 120 l=1,lot
138  a(ib)=0.0
139  a(ib+inc)=0.0
140  ib=ib+jump
141  120 CONTINUE
142  GOTO 140
143 c
144 c postprocessing (isign=-1):
145 c --------------------------
146 c
147  130 CONTINUE
148  CALL fft99b(work,a,trigs,inc,jump,n,lot)
149 c
150  140 CONTINUE
151 
152  END SUBROUTINE fft991
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17