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