1 | !c Copyright 1981-2016 ECMWF. |
---|
2 | !c |
---|
3 | !c This software is licensed under the terms of the Apache Licence |
---|
4 | !c Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
5 | !c |
---|
6 | !c In applying this licence, ECMWF does not waive the privileges and immunities |
---|
7 | !c granted to it by virtue of its status as an intergovernmental organisation |
---|
8 | !c nor does it submit to any jurisdiction. |
---|
9 | !c |
---|
10 | |
---|
11 | SUBROUTINE QPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,ILOT,N,IFAC,ILA,IERR) |
---|
12 | REAL A(*),B(*),C(*),D(*),TRIGS(*) |
---|
13 | !C |
---|
14 | !C SUBROUTINE 'QPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART |
---|
15 | !C OF MULTIPLE REAL FFT (FOURIER ANALYSIS) ROUTINE |
---|
16 | !C |
---|
17 | !C A IS FIRST REAL INPUT VECTOR |
---|
18 | !C EQUIVALENCE B(1) WITH A(IFAC*ILA*INC1+1) |
---|
19 | !C C IS FIRST REAL OUTPUT VECTOR |
---|
20 | !C EQUIVALENCE D(1) WITH C(ILA*INC2+1) |
---|
21 | !C TRIGS IS A PRECALCULATED LIST OF SINES & COSINES |
---|
22 | !C INC1 IS THE ADDRESSING INCREMENT FOR A |
---|
23 | !C INC2 IS THE ADDRESSING INCREMENT FOR C |
---|
24 | !C INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A |
---|
25 | !C INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C |
---|
26 | !C ILOT IS THE NUMBER OF VECTORS |
---|
27 | !C N IS THE LENGTH OF THE VECTORS |
---|
28 | !C IFAC IS THE CURRENT FACTOR OF N |
---|
29 | !C ILA = N/(PRODUCT OF FACTORS USED SO FAR) |
---|
30 | !C IERR IS AN ERROR INDICATOR: |
---|
31 | !C 0 - PASS COMPLETED WITHOUT ERROR |
---|
32 | !C 1 - ILOT GREATER THAN 64 |
---|
33 | !C 2 - IFAC NOT CATERED FOR |
---|
34 | !C 3 - IFAC ONLY CATERED FOR IF ILA=N/IFAC |
---|
35 | !C |
---|
36 | !C----------------------------------------------------------------------- |
---|
37 | !C |
---|
38 | SAVE SIN36, SIN72, QRT5, SIN60 |
---|
39 | |
---|
40 | DATA SIN36/0.587785252292473/,SIN72/0.951056516295154/,& |
---|
41 | QRT5/0.559016994374947/,SIN60/0.866025403784437/ |
---|
42 | |
---|
43 | M=N/IFAC |
---|
44 | IINK=ILA*INC1 |
---|
45 | JINK=ILA*INC2 |
---|
46 | IJUMP=(IFAC-1)*IINK |
---|
47 | KSTOP=(N-IFAC)/(2*IFAC) |
---|
48 | |
---|
49 | IBAD=1 |
---|
50 | IF (ILOT.GT.512) GO TO 910 |
---|
51 | IBASE=0 |
---|
52 | JBASE=0 |
---|
53 | IGO=IFAC-1 |
---|
54 | IF (IGO.EQ.7) IGO=6 |
---|
55 | IBAD=2 |
---|
56 | IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910 |
---|
57 | GO TO (200,300,400,500,600,800),IGO |
---|
58 | |
---|
59 | ! CODING FOR FACTOR 2 |
---|
60 | ! ------------------- |
---|
61 | 200 CONTINUE |
---|
62 | IA=1 |
---|
63 | IB=IA+IINK |
---|
64 | JA=1 |
---|
65 | JB=JA+(2*M-ILA)*INC2 |
---|
66 | |
---|
67 | IF (ILA.EQ.M) GO TO 290 |
---|
68 | |
---|
69 | DO 220 JL=1,ILA |
---|
70 | I=IBASE |
---|
71 | J=JBASE |
---|
72 | !DIR$ IVDEP |
---|
73 | !VOCL LOOP,NOVREC |
---|
74 | DO 210 IJK=1,ILOT |
---|
75 | C(JA+J)=A(IA+I)+A(IB+I) |
---|
76 | C(JB+J)=A(IA+I)-A(IB+I) |
---|
77 | I=I+INC3 |
---|
78 | J=J+INC4 |
---|
79 | 210 CONTINUE |
---|
80 | IBASE=IBASE+INC1 |
---|
81 | JBASE=JBASE+INC2 |
---|
82 | 220 CONTINUE |
---|
83 | JA=JA+JINK |
---|
84 | JINK=2*JINK |
---|
85 | JB=JB-JINK |
---|
86 | IBASE=IBASE+IJUMP |
---|
87 | IJUMP=2*IJUMP+IINK |
---|
88 | IF (JA.EQ.JB) GO TO 260 |
---|
89 | DO 250 K=ILA,KSTOP,ILA |
---|
90 | KB=K+K |
---|
91 | C1=TRIGS(KB+1) |
---|
92 | S1=TRIGS(KB+2) |
---|
93 | JBASE=0 |
---|
94 | DO 240 JL=1,ILA |
---|
95 | I=IBASE |
---|
96 | J=JBASE |
---|
97 | !DIR$ IVDEP |
---|
98 | !VOCL LOOP,NOVREC |
---|
99 | DO 230 IJK=1,ILOT |
---|
100 | C(JA+J)=A(IA+I)+(C1*A(IB+I)+S1*B(IB+I)) |
---|
101 | C(JB+J)=A(IA+I)-(C1*A(IB+I)+S1*B(IB+I)) |
---|
102 | D(JA+J)=(C1*B(IB+I)-S1*A(IB+I))+B(IA+I) |
---|
103 | D(JB+J)=(C1*B(IB+I)-S1*A(IB+I))-B(IA+I) |
---|
104 | I=I+INC3 |
---|
105 | J=J+INC4 |
---|
106 | 230 CONTINUE |
---|
107 | IBASE=IBASE+INC1 |
---|
108 | JBASE=JBASE+INC2 |
---|
109 | 240 CONTINUE |
---|
110 | IBASE=IBASE+IJUMP |
---|
111 | JA=JA+JINK |
---|
112 | JB=JB-JINK |
---|
113 | 250 CONTINUE |
---|
114 | IF (JA.GT.JB) GO TO 900 |
---|
115 | 260 CONTINUE |
---|
116 | JBASE=0 |
---|
117 | DO 280 JL=1,ILA |
---|
118 | I=IBASE |
---|
119 | J=JBASE |
---|
120 | !DIR$ IVDEP |
---|
121 | !VOCL LOOP,NOVREC |
---|
122 | DO 270 IJK=1,ILOT |
---|
123 | C(JA+J)=A(IA+I) |
---|
124 | D(JA+J)=-A(IB+I) |
---|
125 | I=I+INC3 |
---|
126 | J=J+INC4 |
---|
127 | 270 CONTINUE |
---|
128 | IBASE=IBASE+INC1 |
---|
129 | JBASE=JBASE+INC2 |
---|
130 | 280 CONTINUE |
---|
131 | GO TO 900 |
---|
132 | |
---|
133 | 290 CONTINUE |
---|
134 | Z=1.0/FLOAT(N) |
---|
135 | DO 294 JL=1,ILA |
---|
136 | I=IBASE |
---|
137 | J=JBASE |
---|
138 | !DIR$ IVDEP |
---|
139 | !VOCL LOOP,NOVREC |
---|
140 | DO 292 IJK=1,ILOT |
---|
141 | C(JA+J)=Z*(A(IA+I)+A(IB+I)) |
---|
142 | C(JB+J)=Z*(A(IA+I)-A(IB+I)) |
---|
143 | I=I+INC3 |
---|
144 | J=J+INC4 |
---|
145 | 292 CONTINUE |
---|
146 | IBASE=IBASE+INC1 |
---|
147 | JBASE=JBASE+INC2 |
---|
148 | 294 CONTINUE |
---|
149 | GO TO 900 |
---|
150 | |
---|
151 | ! CODING FOR FACTOR 3 |
---|
152 | ! ------------------- |
---|
153 | 300 CONTINUE |
---|
154 | IA=1 |
---|
155 | IB=IA+IINK |
---|
156 | IC=IB+IINK |
---|
157 | JA=1 |
---|
158 | JB=JA+(2*M-ILA)*INC2 |
---|
159 | JC=JB |
---|
160 | |
---|
161 | IF (ILA.EQ.M) GO TO 390 |
---|
162 | |
---|
163 | DO 320 JL=1,ILA |
---|
164 | I=IBASE |
---|
165 | J=JBASE |
---|
166 | !DIR$ IVDEP |
---|
167 | !VOCL LOOP,NOVREC |
---|
168 | DO 310 IJK=1,ILOT |
---|
169 | C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I)) |
---|
170 | C(JB+J)=A(IA+I)-0.5*(A(IB+I)+A(IC+I)) |
---|
171 | D(JB+J)=SIN60*(A(IC+I)-A(IB+I)) |
---|
172 | I=I+INC3 |
---|
173 | J=J+INC4 |
---|
174 | 310 CONTINUE |
---|
175 | IBASE=IBASE+INC1 |
---|
176 | JBASE=JBASE+INC2 |
---|
177 | 320 CONTINUE |
---|
178 | JA=JA+JINK |
---|
179 | JINK=2*JINK |
---|
180 | JB=JB+JINK |
---|
181 | JC=JC-JINK |
---|
182 | IBASE=IBASE+IJUMP |
---|
183 | IJUMP=2*IJUMP+IINK |
---|
184 | IF (JA.EQ.JC) GO TO 360 |
---|
185 | DO 350 K=ILA,KSTOP,ILA |
---|
186 | KB=K+K |
---|
187 | KC=KB+KB |
---|
188 | C1=TRIGS(KB+1) |
---|
189 | S1=TRIGS(KB+2) |
---|
190 | C2=TRIGS(KC+1) |
---|
191 | S2=TRIGS(KC+2) |
---|
192 | JBASE=0 |
---|
193 | DO 340 JL=1,ILA |
---|
194 | I=IBASE |
---|
195 | J=JBASE |
---|
196 | !DIR$ IVDEP |
---|
197 | !VOCL LOOP,NOVREC |
---|
198 | DO 330 IJK=1,ILOT |
---|
199 | A1=(C1*A(IB+I)+S1*B(IB+I))+(C2*A(IC+I)+S2*B(IC+I)) |
---|
200 | B1=(C1*B(IB+I)-S1*A(IB+I))+(C2*B(IC+I)-S2*A(IC+I)) |
---|
201 | A2=A(IA+I)-0.5*A1 |
---|
202 | B2=B(IA+I)-0.5*B1 |
---|
203 | A3=SIN60*((C1*A(IB+I)+S1*B(IB+I))-(C2*A(IC+I)+S2*B(IC+I))) |
---|
204 | B3=SIN60*((C1*B(IB+I)-S1*A(IB+I))-(C2*B(IC+I)-S2*A(IC+I))) |
---|
205 | C(JA+J)=A(IA+I)+A1 |
---|
206 | D(JA+J)=B(IA+I)+B1 |
---|
207 | C(JB+J)=A2+B3 |
---|
208 | D(JB+J)=B2-A3 |
---|
209 | C(JC+J)=A2-B3 |
---|
210 | D(JC+J)=-(B2+A3) |
---|
211 | I=I+INC3 |
---|
212 | J=J+INC4 |
---|
213 | 330 CONTINUE |
---|
214 | IBASE=IBASE+INC1 |
---|
215 | JBASE=JBASE+INC2 |
---|
216 | 340 CONTINUE |
---|
217 | IBASE=IBASE+IJUMP |
---|
218 | JA=JA+JINK |
---|
219 | JB=JB+JINK |
---|
220 | JC=JC-JINK |
---|
221 | 350 CONTINUE |
---|
222 | IF (JA.GT.JC) GO TO 900 |
---|
223 | 360 CONTINUE |
---|
224 | JBASE=0 |
---|
225 | DO 380 JL=1,ILA |
---|
226 | I=IBASE |
---|
227 | J=JBASE |
---|
228 | !DIR$ IVDEP |
---|
229 | !VOCL LOOP,NOVREC |
---|
230 | DO 370 IJK=1,ILOT |
---|
231 | C(JA+J)=A(IA+I)+0.5*(A(IB+I)-A(IC+I)) |
---|
232 | D(JA+J)=-SIN60*(A(IB+I)+A(IC+I)) |
---|
233 | C(JB+J)=A(IA+I)-(A(IB+I)-A(IC+I)) |
---|
234 | I=I+INC3 |
---|
235 | J=J+INC4 |
---|
236 | 370 CONTINUE |
---|
237 | IBASE=IBASE+INC1 |
---|
238 | JBASE=JBASE+INC2 |
---|
239 | 380 CONTINUE |
---|
240 | GO TO 900 |
---|
241 | |
---|
242 | 390 CONTINUE |
---|
243 | Z=1.0/FLOAT(N) |
---|
244 | ZSIN60=Z*SIN60 |
---|
245 | DO 394 JL=1,ILA |
---|
246 | I=IBASE |
---|
247 | J=JBASE |
---|
248 | !DIR$ IVDEP |
---|
249 | !VOCL LOOP,NOVREC |
---|
250 | DO 392 IJK=1,ILOT |
---|
251 | C(JA+J)=Z*(A(IA+I)+(A(IB+I)+A(IC+I))) |
---|
252 | C(JB+J)=Z*(A(IA+I)-0.5*(A(IB+I)+A(IC+I))) |
---|
253 | D(JB+J)=ZSIN60*(A(IC+I)-A(IB+I)) |
---|
254 | I=I+INC3 |
---|
255 | J=J+INC4 |
---|
256 | 392 CONTINUE |
---|
257 | IBASE=IBASE+INC1 |
---|
258 | JBASE=JBASE+INC2 |
---|
259 | 394 CONTINUE |
---|
260 | GO TO 900 |
---|
261 | |
---|
262 | ! CODING FOR FACTOR 4 |
---|
263 | ! ------------------- |
---|
264 | 400 CONTINUE |
---|
265 | IA=1 |
---|
266 | IB=IA+IINK |
---|
267 | IC=IB+IINK |
---|
268 | ID=IC+IINK |
---|
269 | JA=1 |
---|
270 | JB=JA+(2*M-ILA)*INC2 |
---|
271 | JC=JB+2*M*INC2 |
---|
272 | JD=JB |
---|
273 | |
---|
274 | IF (ILA.EQ.M) GO TO 490 |
---|
275 | |
---|
276 | DO 420 JL=1,ILA |
---|
277 | I=IBASE |
---|
278 | J=JBASE |
---|
279 | !DIR$ IVDEP |
---|
280 | !VOCL LOOP,NOVREC |
---|
281 | DO 410 IJK=1,ILOT |
---|
282 | C(JA+J)=(A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)) |
---|
283 | C(JC+J)=(A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)) |
---|
284 | C(JB+J)=A(IA+I)-A(IC+I) |
---|
285 | D(JB+J)=A(ID+I)-A(IB+I) |
---|
286 | I=I+INC3 |
---|
287 | J=J+INC4 |
---|
288 | 410 CONTINUE |
---|
289 | IBASE=IBASE+INC1 |
---|
290 | JBASE=JBASE+INC2 |
---|
291 | 420 CONTINUE |
---|
292 | JA=JA+JINK |
---|
293 | JINK=2*JINK |
---|
294 | JB=JB+JINK |
---|
295 | JC=JC-JINK |
---|
296 | JD=JD-JINK |
---|
297 | IBASE=IBASE+IJUMP |
---|
298 | IJUMP=2*IJUMP+IINK |
---|
299 | IF (JB.EQ.JC) GO TO 460 |
---|
300 | DO 450 K=ILA,KSTOP,ILA |
---|
301 | KB=K+K |
---|
302 | KC=KB+KB |
---|
303 | KD=KC+KB |
---|
304 | C1=TRIGS(KB+1) |
---|
305 | S1=TRIGS(KB+2) |
---|
306 | C2=TRIGS(KC+1) |
---|
307 | S2=TRIGS(KC+2) |
---|
308 | C3=TRIGS(KD+1) |
---|
309 | S3=TRIGS(KD+2) |
---|
310 | JBASE=0 |
---|
311 | DO 440 JL=1,ILA |
---|
312 | I=IBASE |
---|
313 | J=JBASE |
---|
314 | !DIR$ IVDEP |
---|
315 | !VOCL LOOP,NOVREC |
---|
316 | DO 430 IJK=1,ILOT |
---|
317 | A0=A(IA+I)+(C2*A(IC+I)+S2*B(IC+I)) |
---|
318 | A2=A(IA+I)-(C2*A(IC+I)+S2*B(IC+I)) |
---|
319 | A1=(C1*A(IB+I)+S1*B(IB+I))+(C3*A(ID+I)+S3*B(ID+I)) |
---|
320 | A3=(C1*A(IB+I)+S1*B(IB+I))-(C3*A(ID+I)+S3*B(ID+I)) |
---|
321 | B0=B(IA+I)+(C2*B(IC+I)-S2*A(IC+I)) |
---|
322 | B2=B(IA+I)-(C2*B(IC+I)-S2*A(IC+I)) |
---|
323 | B1=(C1*B(IB+I)-S1*A(IB+I))+(C3*B(ID+I)-S3*A(ID+I)) |
---|
324 | B3=(C1*B(IB+I)-S1*A(IB+I))-(C3*B(ID+I)-S3*A(ID+I)) |
---|
325 | C(JA+J)=A0+A1 |
---|
326 | C(JC+J)=A0-A1 |
---|
327 | D(JA+J)=B0+B1 |
---|
328 | D(JC+J)=B1-B0 |
---|
329 | C(JB+J)=A2+B3 |
---|
330 | C(JD+J)=A2-B3 |
---|
331 | D(JB+J)=B2-A3 |
---|
332 | D(JD+J)=-(B2+A3) |
---|
333 | I=I+INC3 |
---|
334 | J=J+INC4 |
---|
335 | 430 CONTINUE |
---|
336 | IBASE=IBASE+INC1 |
---|
337 | JBASE=JBASE+INC2 |
---|
338 | 440 CONTINUE |
---|
339 | IBASE=IBASE+IJUMP |
---|
340 | JA=JA+JINK |
---|
341 | JB=JB+JINK |
---|
342 | JC=JC-JINK |
---|
343 | JD=JD-JINK |
---|
344 | 450 CONTINUE |
---|
345 | IF (JB.GT.JC) GO TO 900 |
---|
346 | 460 CONTINUE |
---|
347 | SIN45=SQRT(0.5) |
---|
348 | JBASE=0 |
---|
349 | DO 480 JL=1,ILA |
---|
350 | I=IBASE |
---|
351 | J=JBASE |
---|
352 | !DIR$ IVDEP |
---|
353 | !VOCL LOOP,NOVREC |
---|
354 | DO 470 IJK=1,ILOT |
---|
355 | C(JA+J)=A(IA+I)+SIN45*(A(IB+I)-A(ID+I)) |
---|
356 | C(JB+J)=A(IA+I)-SIN45*(A(IB+I)-A(ID+I)) |
---|
357 | D(JA+J)=-A(IC+I)-SIN45*(A(IB+I)+A(ID+I)) |
---|
358 | D(JB+J)=A(IC+I)-SIN45*(A(IB+I)+A(ID+I)) |
---|
359 | I=I+INC3 |
---|
360 | J=J+INC4 |
---|
361 | 470 CONTINUE |
---|
362 | IBASE=IBASE+INC1 |
---|
363 | JBASE=JBASE+INC2 |
---|
364 | 480 CONTINUE |
---|
365 | GO TO 900 |
---|
366 | |
---|
367 | 490 CONTINUE |
---|
368 | Z=1.0/FLOAT(N) |
---|
369 | DO 494 JL=1,ILA |
---|
370 | I=IBASE |
---|
371 | J=JBASE |
---|
372 | !DIR$ IVDEP |
---|
373 | !VOCL LOOP,NOVREC |
---|
374 | DO 492 IJK=1,ILOT |
---|
375 | C(JA+J)=Z*((A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I))) |
---|
376 | C(JC+J)=Z*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I))) |
---|
377 | C(JB+J)=Z*(A(IA+I)-A(IC+I)) |
---|
378 | D(JB+J)=Z*(A(ID+I)-A(IB+I)) |
---|
379 | I=I+INC3 |
---|
380 | J=J+INC4 |
---|
381 | 492 CONTINUE |
---|
382 | IBASE=IBASE+INC1 |
---|
383 | JBASE=JBASE+INC2 |
---|
384 | 494 CONTINUE |
---|
385 | GO TO 900 |
---|
386 | |
---|
387 | ! CODING FOR FACTOR 5 |
---|
388 | ! ------------------- |
---|
389 | 500 CONTINUE |
---|
390 | IA=1 |
---|
391 | IB=IA+IINK |
---|
392 | IC=IB+IINK |
---|
393 | ID=IC+IINK |
---|
394 | IE=ID+IINK |
---|
395 | JA=1 |
---|
396 | JB=JA+(2*M-ILA)*INC2 |
---|
397 | JC=JB+2*M*INC2 |
---|
398 | JD=JC |
---|
399 | JE=JB |
---|
400 | |
---|
401 | IF (ILA.EQ.M) GO TO 590 |
---|
402 | |
---|
403 | DO 520 JL=1,ILA |
---|
404 | I=IBASE |
---|
405 | J=JBASE |
---|
406 | !DIR$ IVDEP |
---|
407 | !VOCL LOOP,NOVREC |
---|
408 | DO 510 IJK=1,ILOT |
---|
409 | A1=A(IB+I)+A(IE+I) |
---|
410 | A3=A(IB+I)-A(IE+I) |
---|
411 | A2=A(IC+I)+A(ID+I) |
---|
412 | A4=A(IC+I)-A(ID+I) |
---|
413 | A5=A(IA+I)-0.25*(A1+A2) |
---|
414 | A6=QRT5*(A1-A2) |
---|
415 | C(JA+J)=A(IA+I)+(A1+A2) |
---|
416 | C(JB+J)=A5+A6 |
---|
417 | C(JC+J)=A5-A6 |
---|
418 | D(JB+J)=-SIN72*A3-SIN36*A4 |
---|
419 | D(JC+J)=-SIN36*A3+SIN72*A4 |
---|
420 | I=I+INC3 |
---|
421 | J=J+INC4 |
---|
422 | 510 CONTINUE |
---|
423 | IBASE=IBASE+INC1 |
---|
424 | JBASE=JBASE+INC2 |
---|
425 | 520 CONTINUE |
---|
426 | JA=JA+JINK |
---|
427 | JINK=2*JINK |
---|
428 | JB=JB+JINK |
---|
429 | JC=JC+JINK |
---|
430 | JD=JD-JINK |
---|
431 | JE=JE-JINK |
---|
432 | IBASE=IBASE+IJUMP |
---|
433 | IJUMP=2*IJUMP+IINK |
---|
434 | IF (JB.EQ.JD) GO TO 560 |
---|
435 | DO 550 K=ILA,KSTOP,ILA |
---|
436 | KB=K+K |
---|
437 | KC=KB+KB |
---|
438 | KD=KC+KB |
---|
439 | KE=KD+KB |
---|
440 | C1=TRIGS(KB+1) |
---|
441 | S1=TRIGS(KB+2) |
---|
442 | C2=TRIGS(KC+1) |
---|
443 | S2=TRIGS(KC+2) |
---|
444 | C3=TRIGS(KD+1) |
---|
445 | S3=TRIGS(KD+2) |
---|
446 | C4=TRIGS(KE+1) |
---|
447 | S4=TRIGS(KE+2) |
---|
448 | JBASE=0 |
---|
449 | DO 540 JL=1,ILA |
---|
450 | I=IBASE |
---|
451 | J=JBASE |
---|
452 | !DIR$ IVDEP |
---|
453 | !VOCL LOOP,NOVREC |
---|
454 | DO 530 IJK=1,ILOT |
---|
455 | A1=(C1*A(IB+I)+S1*B(IB+I))+(C4*A(IE+I)+S4*B(IE+I)) |
---|
456 | A3=(C1*A(IB+I)+S1*B(IB+I))-(C4*A(IE+I)+S4*B(IE+I)) |
---|
457 | A2=(C2*A(IC+I)+S2*B(IC+I))+(C3*A(ID+I)+S3*B(ID+I)) |
---|
458 | A4=(C2*A(IC+I)+S2*B(IC+I))-(C3*A(ID+I)+S3*B(ID+I)) |
---|
459 | B1=(C1*B(IB+I)-S1*A(IB+I))+(C4*B(IE+I)-S4*A(IE+I)) |
---|
460 | B3=(C1*B(IB+I)-S1*A(IB+I))-(C4*B(IE+I)-S4*A(IE+I)) |
---|
461 | B2=(C2*B(IC+I)-S2*A(IC+I))+(C3*B(ID+I)-S3*A(ID+I)) |
---|
462 | B4=(C2*B(IC+I)-S2*A(IC+I))-(C3*B(ID+I)-S3*A(ID+I)) |
---|
463 | A5=A(IA+I)-0.25*(A1+A2) |
---|
464 | A6=QRT5*(A1-A2) |
---|
465 | B5=B(IA+I)-0.25*(B1+B2) |
---|
466 | B6=QRT5*(B1-B2) |
---|
467 | A10=A5+A6 |
---|
468 | A20=A5-A6 |
---|
469 | B10=B5+B6 |
---|
470 | B20=B5-B6 |
---|
471 | A11=SIN72*B3+SIN36*B4 |
---|
472 | A21=SIN36*B3-SIN72*B4 |
---|
473 | B11=SIN72*A3+SIN36*A4 |
---|
474 | B21=SIN36*A3-SIN72*A4 |
---|
475 | C(JA+J)=A(IA+I)+(A1+A2) |
---|
476 | C(JB+J)=A10+A11 |
---|
477 | C(JE+J)=A10-A11 |
---|
478 | C(JC+J)=A20+A21 |
---|
479 | C(JD+J)=A20-A21 |
---|
480 | D(JA+J)=B(IA+I)+(B1+B2) |
---|
481 | D(JB+J)=B10-B11 |
---|
482 | D(JE+J)=-(B10+B11) |
---|
483 | D(JC+J)=B20-B21 |
---|
484 | D(JD+J)=-(B20+B21) |
---|
485 | I=I+INC3 |
---|
486 | J=J+INC4 |
---|
487 | 530 CONTINUE |
---|
488 | IBASE=IBASE+INC1 |
---|
489 | JBASE=JBASE+INC2 |
---|
490 | 540 CONTINUE |
---|
491 | IBASE=IBASE+IJUMP |
---|
492 | JA=JA+JINK |
---|
493 | JB=JB+JINK |
---|
494 | JC=JC+JINK |
---|
495 | JD=JD-JINK |
---|
496 | JE=JE-JINK |
---|
497 | 550 CONTINUE |
---|
498 | IF (JB.GT.JD) GO TO 900 |
---|
499 | 560 CONTINUE |
---|
500 | JBASE=0 |
---|
501 | DO 580 JL=1,ILA |
---|
502 | I=IBASE |
---|
503 | J=JBASE |
---|
504 | !DIR$ IVDEP |
---|
505 | !VOCL LOOP,NOVREC |
---|
506 | DO 570 IJK=1,ILOT |
---|
507 | A1=A(IB+I)+A(IE+I) |
---|
508 | A3=A(IB+I)-A(IE+I) |
---|
509 | A2=A(IC+I)+A(ID+I) |
---|
510 | A4=A(IC+I)-A(ID+I) |
---|
511 | A5=A(IA+I)+0.25*(A3-A4) |
---|
512 | A6=QRT5*(A3+A4) |
---|
513 | C(JA+J)=A5+A6 |
---|
514 | C(JB+J)=A5-A6 |
---|
515 | C(JC+J)=A(IA+I)-(A3-A4) |
---|
516 | D(JA+J)=-SIN36*A1-SIN72*A2 |
---|
517 | D(JB+J)=-SIN72*A1+SIN36*A2 |
---|
518 | I=I+INC3 |
---|
519 | J=J+INC4 |
---|
520 | 570 CONTINUE |
---|
521 | IBASE=IBASE+INC1 |
---|
522 | JBASE=JBASE+INC2 |
---|
523 | 580 CONTINUE |
---|
524 | GO TO 900 |
---|
525 | |
---|
526 | 590 CONTINUE |
---|
527 | Z=1.0/FLOAT(N) |
---|
528 | ZQRT5=Z*QRT5 |
---|
529 | ZSIN36=Z*SIN36 |
---|
530 | ZSIN72=Z*SIN72 |
---|
531 | DO 594 JL=1,ILA |
---|
532 | I=IBASE |
---|
533 | J=JBASE |
---|
534 | !DIR$ IVDEP |
---|
535 | !VOCL LOOP,NOVREC |
---|
536 | DO 592 IJK=1,ILOT |
---|
537 | A1=A(IB+I)+A(IE+I) |
---|
538 | A3=A(IB+I)-A(IE+I) |
---|
539 | A2=A(IC+I)+A(ID+I) |
---|
540 | A4=A(IC+I)-A(ID+I) |
---|
541 | A5=Z*(A(IA+I)-0.25*(A1+A2)) |
---|
542 | A6=ZQRT5*(A1-A2) |
---|
543 | C(JA+J)=Z*(A(IA+I)+(A1+A2)) |
---|
544 | C(JB+J)=A5+A6 |
---|
545 | C(JC+J)=A5-A6 |
---|
546 | D(JB+J)=-ZSIN72*A3-ZSIN36*A4 |
---|
547 | D(JC+J)=-ZSIN36*A3+ZSIN72*A4 |
---|
548 | I=I+INC3 |
---|
549 | J=J+INC4 |
---|
550 | 592 CONTINUE |
---|
551 | IBASE=IBASE+INC1 |
---|
552 | JBASE=JBASE+INC2 |
---|
553 | 594 CONTINUE |
---|
554 | GO TO 900 |
---|
555 | |
---|
556 | ! CODING FOR FACTOR 6 |
---|
557 | ! ------------------- |
---|
558 | 600 CONTINUE |
---|
559 | IA=1 |
---|
560 | IB=IA+IINK |
---|
561 | IC=IB+IINK |
---|
562 | ID=IC+IINK |
---|
563 | IE=ID+IINK |
---|
564 | IF=IE+IINK |
---|
565 | JA=1 |
---|
566 | JB=JA+(2*M-ILA)*INC2 |
---|
567 | JC=JB+2*M*INC2 |
---|
568 | JD=JC+2*M*INC2 |
---|
569 | JE=JC |
---|
570 | JF=JB |
---|
571 | |
---|
572 | IF (ILA.EQ.M) GO TO 690 |
---|
573 | |
---|
574 | DO 620 JL=1,ILA |
---|
575 | I=IBASE |
---|
576 | J=JBASE |
---|
577 | !DIR$ IVDEP |
---|
578 | !VOCL LOOP,NOVREC |
---|
579 | DO 610 IJK=1,ILOT |
---|
580 | A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I)) |
---|
581 | C(JA+J)=(A(IA+I)+A(ID+I))+A11 |
---|
582 | C(JC+J)=(A(IA+I)+A(ID+I)-0.5*A11) |
---|
583 | D(JC+J)=SIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I))) |
---|
584 | A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I)) |
---|
585 | C(JB+J)=(A(IA+I)-A(ID+I))-0.5*A11 |
---|
586 | D(JB+J)=SIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) |
---|
587 | C(JD+J)=(A(IA+I)-A(ID+I))+A11 |
---|
588 | I=I+INC3 |
---|
589 | J=J+INC4 |
---|
590 | 610 CONTINUE |
---|
591 | IBASE=IBASE+INC1 |
---|
592 | JBASE=JBASE+INC2 |
---|
593 | 620 CONTINUE |
---|
594 | JA=JA+JINK |
---|
595 | JINK=2*JINK |
---|
596 | JB=JB+JINK |
---|
597 | JC=JC+JINK |
---|
598 | JD=JD-JINK |
---|
599 | JE=JE-JINK |
---|
600 | JF=JF-JINK |
---|
601 | IBASE=IBASE+IJUMP |
---|
602 | IJUMP=2*IJUMP+IINK |
---|
603 | IF (JC.EQ.JD) GO TO 660 |
---|
604 | DO 650 K=ILA,KSTOP,ILA |
---|
605 | KB=K+K |
---|
606 | KC=KB+KB |
---|
607 | KD=KC+KB |
---|
608 | KE=KD+KB |
---|
609 | KF=KE+KB |
---|
610 | C1=TRIGS(KB+1) |
---|
611 | S1=TRIGS(KB+2) |
---|
612 | C2=TRIGS(KC+1) |
---|
613 | S2=TRIGS(KC+2) |
---|
614 | C3=TRIGS(KD+1) |
---|
615 | S3=TRIGS(KD+2) |
---|
616 | C4=TRIGS(KE+1) |
---|
617 | S4=TRIGS(KE+2) |
---|
618 | C5=TRIGS(KF+1) |
---|
619 | S5=TRIGS(KF+2) |
---|
620 | JBASE=0 |
---|
621 | DO 640 JL=1,ILA |
---|
622 | I=IBASE |
---|
623 | J=JBASE |
---|
624 | !DIR$ IVDEP |
---|
625 | !VOCL LOOP,NOVREC |
---|
626 | DO 630 IJK=1,ILOT |
---|
627 | A1=C1*A(IB+I)+S1*B(IB+I) |
---|
628 | B1=C1*B(IB+I)-S1*A(IB+I) |
---|
629 | A2=C2*A(IC+I)+S2*B(IC+I) |
---|
630 | B2=C2*B(IC+I)-S2*A(IC+I) |
---|
631 | A3=C3*A(ID+I)+S3*B(ID+I) |
---|
632 | B3=C3*B(ID+I)-S3*A(ID+I) |
---|
633 | A4=C4*A(IE+I)+S4*B(IE+I) |
---|
634 | B4=C4*B(IE+I)-S4*A(IE+I) |
---|
635 | A5=C5*A(IF+I)+S5*B(IF+I) |
---|
636 | B5=C5*B(IF+I)-S5*A(IF+I) |
---|
637 | A11=(A2+A5)+(A1+A4) |
---|
638 | A20=(A(IA+I)+A3)-0.5*A11 |
---|
639 | A21=SIN60*((A2+A5)-(A1+A4)) |
---|
640 | B11=(B2+B5)+(B1+B4) |
---|
641 | B20=(B(IA+I)+B3)-0.5*B11 |
---|
642 | B21=SIN60*((B2+B5)-(B1+B4)) |
---|
643 | C(JA+J)=(A(IA+I)+A3)+A11 |
---|
644 | D(JA+J)=(B(IA+I)+B3)+B11 |
---|
645 | C(JC+J)=A20-B21 |
---|
646 | D(JC+J)=A21+B20 |
---|
647 | C(JE+J)=A20+B21 |
---|
648 | D(JE+J)=A21-B20 |
---|
649 | A11=(A2-A5)+(A4-A1) |
---|
650 | A20=(A(IA+I)-A3)-0.5*A11 |
---|
651 | A21=SIN60*((A4-A1)-(A2-A5)) |
---|
652 | B11=(B5-B2)-(B4-B1) |
---|
653 | B20=(B3-B(IA+I))-0.5*B11 |
---|
654 | B21=SIN60*((B5-B2)+(B4-B1)) |
---|
655 | C(JB+J)=A20-B21 |
---|
656 | D(JB+J)=A21-B20 |
---|
657 | C(JD+J)=A11+(A(IA+I)-A3) |
---|
658 | D(JD+J)=B11+(B3-B(IA+I)) |
---|
659 | C(JF+J)=A20+B21 |
---|
660 | D(JF+J)=A21+B20 |
---|
661 | I=I+INC3 |
---|
662 | J=J+INC4 |
---|
663 | 630 CONTINUE |
---|
664 | IBASE=IBASE+INC1 |
---|
665 | JBASE=JBASE+INC2 |
---|
666 | 640 CONTINUE |
---|
667 | IBASE=IBASE+IJUMP |
---|
668 | JA=JA+JINK |
---|
669 | JB=JB+JINK |
---|
670 | JC=JC+JINK |
---|
671 | JD=JD-JINK |
---|
672 | JE=JE-JINK |
---|
673 | JF=JF-JINK |
---|
674 | 650 CONTINUE |
---|
675 | IF (JC.GT.JD) GO TO 900 |
---|
676 | 660 CONTINUE |
---|
677 | JBASE=0 |
---|
678 | DO 680 JL=1,ILA |
---|
679 | I=IBASE |
---|
680 | J=JBASE |
---|
681 | !DIR$ IVDEP |
---|
682 | !VOCL LOOP,NOVREC |
---|
683 | DO 670 IJK=1,ILOT |
---|
684 | C(JA+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))+ SIN60*(A(IB+I)-A(IF+I)) |
---|
685 | D(JA+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))-SIN60*(A(IC+I)+A(IE+I)) |
---|
686 | C(JB+J)=A(IA+I)-(A(IC+I)-A(IE+I)) |
---|
687 | D(JB+J)=A(ID+I)-(A(IB+I)+A(IF+I)) |
---|
688 | C(JC+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))-SIN60*(A(IB+I)-A(IF+I)) |
---|
689 | D(JC+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))+SIN60*(A(IC+I)+A(IE+I)) |
---|
690 | I=I+INC3 |
---|
691 | J=J+INC4 |
---|
692 | 670 CONTINUE |
---|
693 | IBASE=IBASE+INC1 |
---|
694 | JBASE=JBASE+INC2 |
---|
695 | 680 CONTINUE |
---|
696 | GO TO 900 |
---|
697 | |
---|
698 | 690 CONTINUE |
---|
699 | Z=1.0/FLOAT(N) |
---|
700 | ZSIN60=Z*SIN60 |
---|
701 | DO 694 JL=1,ILA |
---|
702 | I=IBASE |
---|
703 | J=JBASE |
---|
704 | !DIR$ IVDEP |
---|
705 | !VOCL LOOP,NOVREC |
---|
706 | DO 692 IJK=1,ILOT |
---|
707 | A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I)) |
---|
708 | C(JA+J)=Z*((A(IA+I)+A(ID+I))+A11) |
---|
709 | C(JC+J)=Z*((A(IA+I)+A(ID+I))-0.5*A11) |
---|
710 | D(JC+J)=ZSIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I))) |
---|
711 | A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I)) |
---|
712 | C(JB+J)=Z*((A(IA+I)-A(ID+I))-0.5*A11) |
---|
713 | D(JB+J)=ZSIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I))) |
---|
714 | C(JD+J)=Z*((A(IA+I)-A(ID+I))+A11) |
---|
715 | I=I+INC3 |
---|
716 | J=J+INC4 |
---|
717 | 692 CONTINUE |
---|
718 | IBASE=IBASE+INC1 |
---|
719 | JBASE=JBASE+INC2 |
---|
720 | 694 CONTINUE |
---|
721 | GO TO 900 |
---|
722 | |
---|
723 | ! CODING FOR FACTOR 8 |
---|
724 | ! ------------------- |
---|
725 | 800 CONTINUE |
---|
726 | IBAD=3 |
---|
727 | IF (ILA.NE.M) GO TO 910 |
---|
728 | IA=1 |
---|
729 | IB=IA+IINK |
---|
730 | IC=IB+IINK |
---|
731 | ID=IC+IINK |
---|
732 | IE=ID+IINK |
---|
733 | IF=IE+IINK |
---|
734 | IG=IF+IINK |
---|
735 | IH=IG+IINK |
---|
736 | JA=1 |
---|
737 | JB=JA+ILA*INC2 |
---|
738 | JC=JB+2*M*INC2 |
---|
739 | JD=JC+2*M*INC2 |
---|
740 | JE=JD+2*M*INC2 |
---|
741 | Z=1.0/FLOAT(N) |
---|
742 | ZSIN45=Z*SQRT(0.5) |
---|
743 | |
---|
744 | DO 820 JL=1,ILA |
---|
745 | I=IBASE |
---|
746 | J=JBASE |
---|
747 | !DIR$ IVDEP |
---|
748 | !VOCL LOOP,NOVREC |
---|
749 | DO 810 IJK=1,ILOT |
---|
750 | C(JA+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))+& |
---|
751 | ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I)))) |
---|
752 | C(JE+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))-& |
---|
753 | ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I)))) |
---|
754 | C(JC+J)=Z*((A(IA+I)+A(IE+I))-(A(IC+I)+A(IG+I))) |
---|
755 | D(JC+J)=Z*((A(ID+I)+A(IH+I))-(A(IB+I)+A(IF+I))) |
---|
756 | C(JB+J)=Z*(A(IA+I)-A(IE+I))& |
---|
757 | +ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I))) |
---|
758 | C(JD+J)=Z*(A(IA+I)-A(IE+I))& |
---|
759 | -ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I))) |
---|
760 | D(JB+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))& |
---|
761 | +Z*(A(IG+I)-A(IC+I)) |
---|
762 | D(JD+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))& |
---|
763 | -Z*(A(IG+I)-A(IC+I)) |
---|
764 | I=I+INC3 |
---|
765 | J=J+INC4 |
---|
766 | 810 CONTINUE |
---|
767 | IBASE=IBASE+INC1 |
---|
768 | JBASE=JBASE+INC2 |
---|
769 | 820 CONTINUE |
---|
770 | |
---|
771 | ! RETURN |
---|
772 | ! ------ |
---|
773 | 900 CONTINUE |
---|
774 | IBAD=0 |
---|
775 | 910 CONTINUE |
---|
776 | IERR=IBAD |
---|
777 | RETURN |
---|
778 | END |
---|