source: flex_extract.git/Source/Fortran/rpassm.f90 @ 4f24798

dev
Last change on this file since 4f24798 was 4f24798, checked in by Anne Tipka <anne.tipka@…>, 22 months ago

elimination of emoslib in fortran code. See ticket #312

  • Property mode set to 100644
File size: 19.8 KB
Line 
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
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG