source: flex_extract.git/Source/Fortran/qpassm.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: 17.7 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 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
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG