V3FIT
vpassm.f
1  SUBROUTINE vpassm(a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot,
2  1 n, ifac, la)
3  USE stel_kinds
4  IMPLICIT NONE
5 C-----------------------------------------------
6 C L o c a l P a r a m e t e r s
7 C-----------------------------------------------
8  REAL(rprec), PARAMETER :: SIN36=0.587785252292473_dp
9  REAL(rprec), PARAMETER :: COS36=0.809016994374947_dp
10  REAL(rprec), PARAMETER :: SIN72=0.951056516295154_dp
11  REAL(rprec), PARAMETER :: COS72=0.309016994374947_dp
12  REAL(rprec), PARAMETER :: SIN60=0.866025403784437_dp
13  REAL(rprec), PARAMETER :: p5 = 0.5_dp
14 C-----------------------------------------------
15 C D u m m y A r g u m e n t s
16 C-----------------------------------------------
17  INTEGER inc1, inc2, inc3, inc4, lot, n, ifac, la
18  REAL(rprec), DIMENSION(*) :: a, b, c, d
19  REAL(rprec), DIMENSION(*) :: trigs
20 C-----------------------------------------------
21 C L o c a l V a r i a b l e s
22 C-----------------------------------------------
23  INTEGER :: m, iink, jink, jump, ibase, jbase, igo, ia, ja, ib, jb
24  1 , l, i, j, ijk, la1, k, kb, ic, jc, kc, id, jd, kd, ie, je, ke
25  REAL(rprec) :: c1,s1,c2,s2,c3,s3,c4,s4
26 C-----------------------------------------------
27 c
28 c SUBROUTINE 'vpassm' - multiple version of 'vpassa'
29 c performs one pass through data
30 c as part of multiple complex (inverse) fft routine
31 c a is first real input vector
32 c b is first imaginary input vector
33 c c is first real output vector
34 c d is first imaginary output vector
35 c trigs is precalculated table of sines & cosines
36 c inc1 is addressing increment for a and b
37 c inc2 is addressing increment for c and d
38 c inc3 is addressing increment between a's & b's
39 c inc4 is addressing increment between c's & d's
40 c lot is the number of vectors
41 c n is length of vectors
42 c ifac is current factor of n
43 c la is product of previous factors
44 c
45 c
46  m=n/ifac
47  iink=m*inc1
48  jink=la*inc2
49  jump=(ifac-1)*jink
50  ibase=0
51  jbase=0
52  igo=ifac-1
53  IF (igo .gt. 4) RETURN
54  GOTO (10,50,90,130),igo
55 c
56 c coding for factor 2
57 c
58  10 ia=1
59  ja=1
60  ib=ia+iink
61  jb=ja+jink
62  DO 20 l=1,la
63  i=ibase
64  j=jbase
65 cdir$ ivdep
66  DO 15 ijk=1,lot
67  c(ja+j)=a(ia+i)+a(ib+i)
68  d(ja+j)=b(ia+i)+b(ib+i)
69  c(jb+j)=a(ia+i)-a(ib+i)
70  d(jb+j)=b(ia+i)-b(ib+i)
71  i=i+inc3
72  j=j+inc4
73  15 CONTINUE
74  ibase=ibase+inc1
75  jbase=jbase+inc2
76  20 CONTINUE
77  IF (la.eq.m) RETURN
78  la1=la+1
79  jbase=jbase+jump
80  DO 40 k=la1,m,la
81  kb=k+k-2
82  c1=trigs(kb+1)
83  s1=trigs(kb+2)
84  DO 30 l=1,la
85  i=ibase
86  j=jbase
87 cdir$ ivdep
88  DO 25 ijk=1,lot
89  c(ja+j)=a(ia+i)+a(ib+i)
90  d(ja+j)=b(ia+i)+b(ib+i)
91  c(jb+j)=c1*(a(ia+i)-a(ib+i))-s1*(b(ia+i)-b(ib+i))
92  d(jb+j)=s1*(a(ia+i)-a(ib+i))+c1*(b(ia+i)-b(ib+i))
93  i=i+inc3
94  j=j+inc4
95  25 CONTINUE
96  ibase=ibase+inc1
97  jbase=jbase+inc2
98  30 CONTINUE
99  jbase=jbase+jump
100  40 CONTINUE
101  RETURN
102 c
103 c coding for factor 3
104 c
105  50 ia=1
106  ja=1
107  ib=ia+iink
108  jb=ja+jink
109  ic=ib+iink
110  jc=jb+jink
111  DO 60 l=1,la
112  i=ibase
113  j=jbase
114 cdir$ ivdep
115  DO 55 ijk=1,lot
116  c(ja+j)=a(ia+i)+(a(ib+i)+a(ic+i))
117  d(ja+j)=b(ia+i)+(b(ib+i)+b(ic+i))
118  c(jb+j)=(a(ia+i)-p5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)-b(ic+i)))
119  c(jc+j)=(a(ia+i)-p5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)-b(ic+i)))
120  d(jb+j)=(b(ia+i)-p5*(b(ib+i)+b(ic+i)))+(sin60*(a(ib+i)-a(ic+i)))
121  d(jc+j)=(b(ia+i)-p5*(b(ib+i)+b(ic+i)))-(sin60*(a(ib+i)-a(ic+i)))
122  i=i+inc3
123  j=j+inc4
124  55 CONTINUE
125  ibase=ibase+inc1
126  jbase=jbase+inc2
127  60 CONTINUE
128  IF (la.eq.m) RETURN
129  la1=la+1
130  jbase=jbase+jump
131  DO 80 k=la1,m,la
132  kb=k+k-2
133  kc=kb+kb
134  c1=trigs(kb+1)
135  s1=trigs(kb+2)
136  c2=trigs(kc+1)
137  s2=trigs(kc+2)
138  DO 70 l=1,la
139  i=ibase
140  j=jbase
141 cdir$ ivdep
142  DO 65 ijk=1,lot
143  c(ja+j)=a(ia+i)+(a(ib+i)+a(ic+i))
144  d(ja+j)=b(ia+i)+(b(ib+i)+b(ic+i))
145  c(jb+j)=
146  * c1*((a(ia+i)-p5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)-b(ic+i))))
147  * -s1*((b(ia+i)-p5*(b(ib+i)+b(ic+i)))+(sin60*(a(ib+i)-a(ic+i))))
148  d(jb+j)=
149  * s1*((a(ia+i)-p5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)-b(ic+i))))
150  * +c1*((b(ia+i)-p5*(b(ib+i)+b(ic+i)))+(sin60*(a(ib+i)-a(ic+i))))
151  c(jc+j)=
152  * c2*((a(ia+i)-p5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)-b(ic+i))))
153  * -s2*((b(ia+i)-p5*(b(ib+i)+b(ic+i)))-(sin60*(a(ib+i)-a(ic+i))))
154  d(jc+j)=
155  * s2*((a(ia+i)-p5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)-b(ic+i))))
156  * +c2*((b(ia+i)-p5*(b(ib+i)+b(ic+i)))-(sin60*(a(ib+i)-a(ic+i))))
157  i=i+inc3
158  j=j+inc4
159  65 CONTINUE
160  ibase=ibase+inc1
161  jbase=jbase+inc2
162  70 CONTINUE
163  jbase=jbase+jump
164  80 CONTINUE
165  RETURN
166 c
167 c coding for factor 4
168 c
169  90 ia=1
170  ja=1
171  ib=ia+iink
172  jb=ja+jink
173  ic=ib+iink
174  jc=jb+jink
175  id=ic+iink
176  jd=jc+jink
177  DO 100 l=1,la
178  i=ibase
179  j=jbase
180 cdir$ ivdep
181  DO 95 ijk=1,lot
182  c(ja+j)=(a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i))
183  c(jc+j)=(a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))
184  d(ja+j)=(b(ia+i)+b(ic+i))+(b(ib+i)+b(id+i))
185  d(jc+j)=(b(ia+i)+b(ic+i))-(b(ib+i)+b(id+i))
186  c(jb+j)=(a(ia+i)-a(ic+i))-(b(ib+i)-b(id+i))
187  c(jd+j)=(a(ia+i)-a(ic+i))+(b(ib+i)-b(id+i))
188  d(jb+j)=(b(ia+i)-b(ic+i))+(a(ib+i)-a(id+i))
189  d(jd+j)=(b(ia+i)-b(ic+i))-(a(ib+i)-a(id+i))
190  i=i+inc3
191  j=j+inc4
192  95 CONTINUE
193  ibase=ibase+inc1
194  jbase=jbase+inc2
195  100 CONTINUE
196  IF (la.eq.m) RETURN
197  la1=la+1
198  jbase=jbase+jump
199  DO 120 k=la1,m,la
200  kb=k+k-2
201  kc=kb+kb
202  kd=kc+kb
203  c1=trigs(kb+1)
204  s1=trigs(kb+2)
205  c2=trigs(kc+1)
206  s2=trigs(kc+2)
207  c3=trigs(kd+1)
208  s3=trigs(kd+2)
209  DO 110 l=1,la
210  i=ibase
211  j=jbase
212 cdir$ ivdep
213  DO 105 ijk=1,lot
214  c(ja+j)=(a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i))
215  d(ja+j)=(b(ia+i)+b(ic+i))+(b(ib+i)+b(id+i))
216  c(jc+j)=
217  * c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
218  * -s2*((b(ia+i)+b(ic+i))-(b(ib+i)+b(id+i)))
219  d(jc+j)=
220  * s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
221  * +c2*((b(ia+i)+b(ic+i))-(b(ib+i)+b(id+i)))
222  c(jb+j)=
223  * c1*((a(ia+i)-a(ic+i))-(b(ib+i)-b(id+i)))
224  * -s1*((b(ia+i)-b(ic+i))+(a(ib+i)-a(id+i)))
225  d(jb+j)=
226  * s1*((a(ia+i)-a(ic+i))-(b(ib+i)-b(id+i)))
227  * +c1*((b(ia+i)-b(ic+i))+(a(ib+i)-a(id+i)))
228  c(jd+j)=
229  * c3*((a(ia+i)-a(ic+i))+(b(ib+i)-b(id+i)))
230  * -s3*((b(ia+i)-b(ic+i))-(a(ib+i)-a(id+i)))
231  d(jd+j)=
232  * s3*((a(ia+i)-a(ic+i))+(b(ib+i)-b(id+i)))
233  * +c3*((b(ia+i)-b(ic+i))-(a(ib+i)-a(id+i)))
234  i=i+inc3
235  j=j+inc4
236  105 CONTINUE
237  ibase=ibase+inc1
238  jbase=jbase+inc2
239  110 CONTINUE
240  jbase=jbase+jump
241  120 CONTINUE
242  RETURN
243 c
244 c coding for factor 5
245 c
246  130 ia=1
247  ja=1
248  ib=ia+iink
249  jb=ja+jink
250  ic=ib+iink
251  jc=jb+jink
252  id=ic+iink
253  jd=jc+jink
254  ie=id+iink
255  je=jd+jink
256  DO 140 l=1,la
257  i=ibase
258  j=jbase
259 cdir$ ivdep
260  DO 135 ijk=1,lot
261  c(ja+j)=a(ia+i)+(a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))
262  d(ja+j)=b(ia+i)+(b(ib+i)+b(ie+i))+(b(ic+i)+b(id+i))
263  c(jb+j)=(a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i)))
264  * -(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))
265  c(je+j)=(a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i)))
266  * +(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i)))
267  d(jb+j)=(b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i)))
268  * +(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i)))
269  d(je+j)=(b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i)))
270  * -(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i)))
271  c(jc+j)=(a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i)))
272  * -(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))
273  c(jd+j)=(a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i)))
274  * +(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i)))
275  d(jc+j)=(b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))
276  * +(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i)))
277  d(jd+j)=(b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))
278  * -(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i)))
279  i=i+inc3
280  j=j+inc4
281  135 CONTINUE
282  ibase=ibase+inc1
283  jbase=jbase+inc2
284  140 CONTINUE
285  IF (la.eq.m) RETURN
286  la1=la+1
287  jbase=jbase+jump
288  DO 160 k=la1,m,la
289  kb=k+k-2
290  kc=kb+kb
291  kd=kc+kb
292  ke=kd+kb
293  c1=trigs(kb+1)
294  s1=trigs(kb+2)
295  c2=trigs(kc+1)
296  s2=trigs(kc+2)
297  c3=trigs(kd+1)
298  s3=trigs(kd+2)
299  c4=trigs(ke+1)
300  s4=trigs(ke+2)
301  DO 150 l=1,la
302  i=ibase
303  j=jbase
304 cdir$ ivdep
305  DO 145 ijk=1,lot
306  c(ja+j)=a(ia+i)+(a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i))
307  d(ja+j)=b(ia+i)+(b(ib+i)+b(ie+i))+(b(ic+i)+b(id+i))
308  c(jb+j)=
309  * c1*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i)))
310  * -(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i))))
311  * -s1*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i)))
312  * +(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
313  d(jb+j)=
314  * s1*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i)))
315  * -(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i))))
316  * +c1*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i)))
317  * +(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
318  c(je+j)=
319  * c4*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i)))
320  * +(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i))))
321  * -s4*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i)))
322  * -(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
323  d(je+j)=
324  * s4*((a(ia+i)+cos72*(a(ib+i)+a(ie+i))-cos36*(a(ic+i)+a(id+i)))
325  * +(sin72*(b(ib+i)-b(ie+i))+sin36*(b(ic+i)-b(id+i))))
326  * +c4*((b(ia+i)+cos72*(b(ib+i)+b(ie+i))-cos36*(b(ic+i)+b(id+i)))
327  * -(sin72*(a(ib+i)-a(ie+i))+sin36*(a(ic+i)-a(id+i))))
328  c(jc+j)=
329  * c2*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i)))
330  * -(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i))))
331  * -s2*((b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))
332  * +(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
333  d(jc+j)=
334  * s2*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i)))
335  * -(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i))))
336  * +c2*((b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))
337  * +(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
338  c(jd+j)=
339  * c3*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i)))
340  * +(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i))))
341  * -s3*((b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))
342  * -(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
343  d(jd+j)=
344  * s3*((a(ia+i)-cos36*(a(ib+i)+a(ie+i))+cos72*(a(ic+i)+a(id+i)))
345  * +(sin36*(b(ib+i)-b(ie+i))-sin72*(b(ic+i)-b(id+i))))
346  * +c3*((b(ia+i)-cos36*(b(ib+i)+b(ie+i))+cos72*(b(ic+i)+b(id+i)))
347  * -(sin36*(a(ib+i)-a(ie+i))-sin72*(a(ic+i)-a(id+i))))
348  i=i+inc3
349  j=j+inc4
350  145 CONTINUE
351  ibase=ibase+inc1
352  jbase=jbase+inc2
353  150 CONTINUE
354  jbase=jbase+jump
355  160 CONTINUE
356 
357  END SUBROUTINE vpassm
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17