Update version info for release v4.6.1 (#2122)
[WRF.git] / var / da / da_ffts / qpassm.inc
blob56d29b86d8e88039e155a6ecd62d7a0578872bb1
1 !C     SUBROUTINE 'QPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART!C     OF MULTIPLE REAL FFT (FOURIER ANALYSIS) ROUTINE
2 !C
3 !C     A IS FIRST REAL INPUT VECTOR
4 !C         EQUIVALENCE B(1) WITH A(IFAC*LA*INC1+1)
5 !C     C IS FIRST REAL OUTPUT VECTOR
6 !C         EQUIVALENCE D(1) WITH C(LA*INC2+1)
7 !C     TRIGS IS A PRECALCULATED LIST OF SINES & COSINES
8 !C     INC1 IS THE ADDRESSING INCREMENT FOR A
9 !C     INC2 IS THE ADDRESSING INCREMENT FOR C
10 !C     INC3 IS THE INCREMENT BETWEEN INPUT VECTORS A
11 !C     INC4 IS THE INCREMENT BETWEEN OUTPUT VECTORS C
12 !C     LOT IS THE NUMBER OF VECTORS
13 !C     N IS THE LENGTH OF THE VECTORS
14 !C     IFAC IS THE CURRENT FACTOR OF N
15 !C     LA = N/(PRODUCT OF FACTORS USED SO FAR)
16 !C     IERR IS AN ERROR INDICATOR:
17 !C              0 - PASS COMPLETED WITHOUT ERROR
18 !C              1 - LOT GREATER THAN VectorLength
19 !C              2 - IFAC NOT CATERED FOR
20 !C              3 - IFAC ONLY CATERED FOR IF LA=N/IFAC
22 !C-----------------------------------------------------------------------
24      SUBROUTINE QPASSM(A,B,C,D,TRIGS,INC1,INC2,INC3,INC4,LOT,N,IFAC,LA,IERR)
26       INTEGER :: inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr
28       REAL  :: a, b, c, d, trigs
29       DIMENSION A(*),B(*),C(*),D(*),TRIGS(N)
30 ! Local named constants
31  CHARACTER (LEN=*), PARAMETER :: RoutineName = "QPASSM"
33       REAL  :: SIN36, SIN72, QRT5, SIN60
34       DATA SIN36/0.587785252292473/,SIN72/0.951056516295154/, &
35           QRT5/0.559016994374947/,SIN60/0.866025403784437/
37       REAL    :: s1, s2, s3, s4, s5
38       REAL    :: sin45, zsin36, zsin72, zqrt5, zsin60, zsin45, z
39       REAL    :: a0, a1, a2, a3, a4, a5, a6, a10, a11, a20, a21
40       REAL    :: b0, b1, b2, b3, b4, b5, b6, b10, b11, b20, b21
41 !      REAL    :: c0, c1, c2, c3, c4, c5
42       REAL    :: c1, c2, c3, c4, c5
43       INTEGER :: i, ijk, l, k, kb, m, iink, jink, ijump, kstop
44       INTEGER :: ibad, igo, ia, ie, je, ibase, jbase, ja, jb, j, ic
45       INTEGER :: if, jf, kf, ib, jc, kc, id, jd, kd, ke, ig, ih
46       INTEGER :: vectorlength
48 !- End of header ---------------------------------------------------------------
51       M=N/IFAC
52       IINK=LA*INC1
53       JINK=LA*INC2
54       IJUMP=(IFAC-1)*IINK
55       KSTOP=(N-IFAC)/(2*IFAC)
57       IBAD=1
58       VectorLength = lot
59       IF (LOT.GT.VectorLength) GO TO 910
60       IBASE=0
61       JBASE=0
62       IGO=IFAC-1
63       IF (IGO.EQ.7) IGO=6
64       IBAD=2
65       IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910
66       GO TO (200,300,400,500,600,800),IGO
68 !     CODING FOR FACTOR 2
69 !     -------------------
70   200 CONTINUE
71       IA=1
72       IB=IA+IINK
73       JA=1
74       JB=JA+(2*M-LA)*INC2
76       IF (LA.EQ.M) GO TO 290
78       DO 220 L=1,LA
79       I=IBASE
80       J=JBASE
81 !DIR$ IVDEP
82       DO 210 IJK=1,LOT
83       C(JA+J)=A(IA+I)+A(IB+I)
84       C(JB+J)=A(IA+I)-A(IB+I)
85       I=I+INC3
86       J=J+INC4
87   210 CONTINUE
88       IBASE=IBASE+INC1
89       JBASE=JBASE+INC2
90   220 CONTINUE
91       JA=JA+JINK
92       JINK=2*JINK
93       JB=JB-JINK
94       IBASE=IBASE+IJUMP
95       IJUMP=2*IJUMP+IINK
96       IF (JA.EQ.JB) GO TO 260
97       DO 250 K=LA,KSTOP,LA
98       KB=K+K
99       C1=TRIGS(KB+1)
100       S1=TRIGS(KB+2)
101       JBASE=0
102       DO 240 L=1,LA
103       I=IBASE
104       J=JBASE
105 !DIR$ IVDEP
106       DO 230 IJK=1,LOT
107       C(JA+J)=A(IA+I)+(C1*A(IB+I)+S1*B(IB+I))
108       C(JB+J)=A(IA+I)-(C1*A(IB+I)+S1*B(IB+I))
109       D(JA+J)=(C1*B(IB+I)-S1*A(IB+I))+B(IA+I)
110       D(JB+J)=(C1*B(IB+I)-S1*A(IB+I))-B(IA+I)
111       I=I+INC3
112       J=J+INC4
113   230 CONTINUE
114       IBASE=IBASE+INC1
115       JBASE=JBASE+INC2
116   240 CONTINUE
117       IBASE=IBASE+IJUMP
118       JA=JA+JINK
119       JB=JB-JINK
120   250 CONTINUE
121       IF (JA.GT.JB) GO TO 900
122   260 CONTINUE
123       JBASE=0
124       DO 280 L=1,LA
125       I=IBASE
126       J=JBASE
127 !DIR$ IVDEP
128       DO 270 IJK=1,LOT
129       C(JA+J)=A(IA+I)
130       D(JA+J)=-A(IB+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
139   290 CONTINUE
140       Z=1.0/FLOAT(N)
141       DO 294 L=1,LA
142       I=IBASE
143       J=JBASE
144 !DIR$ IVDEP
145       DO 292 IJK=1,LOT
146       C(JA+J)=Z*(A(IA+I)+A(IB+I))
147       C(JB+J)=Z*(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
156 !     CODING FOR FACTOR 3
157 !     -------------------
158   300 CONTINUE
159       IA=1
160       IB=IA+IINK
161       IC=IB+IINK
162       JA=1
163       JB=JA+(2*M-LA)*INC2
164       JC=JB
166       IF (LA.EQ.M) GO TO 390
168       DO 320 L=1,LA
169       I=IBASE
170       J=JBASE
171 !DIR$ IVDEP
172       DO 310 IJK=1,LOT
173       C(JA+J)=A(IA+I)+(A(IB+I)+A(IC+I))
174       C(JB+J)=A(IA+I)-0.5*(A(IB+I)+A(IC+I))
175       D(JB+J)=SIN60*(A(IC+I)-A(IB+I))
176       I=I+INC3
177       J=J+INC4
178   310 CONTINUE
179       IBASE=IBASE+INC1
180       JBASE=JBASE+INC2
181   320 CONTINUE
182       JA=JA+JINK
183       JINK=2*JINK
184       JB=JB+JINK
185       JC=JC-JINK
186       IBASE=IBASE+IJUMP
187       IJUMP=2*IJUMP+IINK
188       IF (JA.EQ.JC) GO TO 360
189       DO 350 K=LA,KSTOP,LA
190       KB=K+K
191       KC=KB+KB
192       C1=TRIGS(KB+1)
193       S1=TRIGS(KB+2)
194       C2=TRIGS(KC+1)
195       S2=TRIGS(KC+2)
196       JBASE=0
197       DO 340 L=1,LA
198       I=IBASE
199       J=JBASE
200 !DIR$ IVDEP
201       DO 330 IJK=1,LOT
202       A1=(C1*A(IB+I)+S1*B(IB+I))+(C2*A(IC+I)+S2*B(IC+I))
203       B1=(C1*B(IB+I)-S1*A(IB+I))+(C2*B(IC+I)-S2*A(IC+I))
204       A2=A(IA+I)-0.5*A1
205       B2=B(IA+I)-0.5*B1
206       A3=SIN60*((C1*A(IB+I)+S1*B(IB+I))-(C2*A(IC+I)+S2*B(IC+I)))
207       B3=SIN60*((C1*B(IB+I)-S1*A(IB+I))-(C2*B(IC+I)-S2*A(IC+I)))
208       C(JA+J)=A(IA+I)+A1
209       D(JA+J)=B(IA+I)+B1
210       C(JB+J)=A2+B3
211       D(JB+J)=B2-A3
212       C(JC+J)=A2-B3
213       D(JC+J)=-(B2+A3)
214       I=I+INC3
215       J=J+INC4
216   330 CONTINUE
217       IBASE=IBASE+INC1
218       JBASE=JBASE+INC2
219   340 CONTINUE
220       IBASE=IBASE+IJUMP
221       JA=JA+JINK
222       JB=JB+JINK
223       JC=JC-JINK
224   350 CONTINUE
225       IF (JA.GT.JC) GO TO 900
226   360 CONTINUE
227       JBASE=0
228       DO 380 L=1,LA
229       I=IBASE
230       J=JBASE
231 !DIR$ IVDEP
232       DO 370 IJK=1,LOT
233       C(JA+J)=A(IA+I)+0.5*(A(IB+I)-A(IC+I))
234       D(JA+J)=-SIN60*(A(IB+I)+A(IC+I))
235       C(JB+J)=A(IA+I)-(A(IB+I)-A(IC+I))
236       I=I+INC3
237       J=J+INC4
238   370 CONTINUE
239       IBASE=IBASE+INC1
240       JBASE=JBASE+INC2
241   380 CONTINUE
242       GO TO 900
244   390 CONTINUE
245       Z=1.0/FLOAT(N)
246       ZSIN60=Z*SIN60
247       DO 394 L=1,LA
248       I=IBASE
249       J=JBASE
250 !DIR$ IVDEP
251       DO 392 IJK=1,LOT
252       C(JA+J)=Z*(A(IA+I)+(A(IB+I)+A(IC+I)))
253       C(JB+J)=Z*(A(IA+I)-0.5*(A(IB+I)+A(IC+I)))
254       D(JB+J)=ZSIN60*(A(IC+I)-A(IB+I))
255       I=I+INC3
256       J=J+INC4
257   392 CONTINUE
258       IBASE=IBASE+INC1
259       JBASE=JBASE+INC2
260   394 CONTINUE
261       GO TO 900
263 !     CODING FOR FACTOR 4
264 !     -------------------
265   400 CONTINUE
266       IA=1
267       IB=IA+IINK
268       IC=IB+IINK
269       ID=IC+IINK
270       JA=1
271       JB=JA+(2*M-LA)*INC2
272       JC=JB+2*M*INC2
273       JD=JB
275       IF (LA.EQ.M) GO TO 490
277       DO 420 L=1,LA
278       I=IBASE
279       J=JBASE
280 !DIR$ IVDEP
281       DO 410 IJK=1,LOT
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=LA,KSTOP,LA
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 L=1,LA
312       I=IBASE
313       J=JBASE
314 !DIR$ IVDEP
315       DO 430 IJK=1,LOT
316       A0=A(IA+I)+(C2*A(IC+I)+S2*B(IC+I))
317       A2=A(IA+I)-(C2*A(IC+I)+S2*B(IC+I))
318       A1=(C1*A(IB+I)+S1*B(IB+I))+(C3*A(ID+I)+S3*B(ID+I))
319       A3=(C1*A(IB+I)+S1*B(IB+I))-(C3*A(ID+I)+S3*B(ID+I))
320       B0=B(IA+I)+(C2*B(IC+I)-S2*A(IC+I))
321       B2=B(IA+I)-(C2*B(IC+I)-S2*A(IC+I))
322       B1=(C1*B(IB+I)-S1*A(IB+I))+(C3*B(ID+I)-S3*A(ID+I))
323       B3=(C1*B(IB+I)-S1*A(IB+I))-(C3*B(ID+I)-S3*A(ID+I))
324       C(JA+J)=A0+A1
325       C(JC+J)=A0-A1
326       D(JA+J)=B0+B1
327       D(JC+J)=B1-B0
328       C(JB+J)=A2+B3
329       C(JD+J)=A2-B3
330       D(JB+J)=B2-A3
331       D(JD+J)=-(B2+A3)
332       I=I+INC3
333       J=J+INC4
334   430 CONTINUE
335       IBASE=IBASE+INC1
336       JBASE=JBASE+INC2
337   440 CONTINUE
338       IBASE=IBASE+IJUMP
339       JA=JA+JINK
340       JB=JB+JINK
341       JC=JC-JINK
342       JD=JD-JINK
343   450 CONTINUE
344       IF (JB.GT.JC) GO TO 900
345   460 CONTINUE
346       SIN45=SQRT(0.5)
347       JBASE=0
348       DO 480 L=1,LA
349       I=IBASE
350       J=JBASE
351 !DIR$ IVDEP
352       DO 470 IJK=1,LOT
353       C(JA+J)=A(IA+I)+SIN45*(A(IB+I)-A(ID+I))
354       C(JB+J)=A(IA+I)-SIN45*(A(IB+I)-A(ID+I))
355       D(JA+J)=-A(IC+I)-SIN45*(A(IB+I)+A(ID+I))
356       D(JB+J)=A(IC+I)-SIN45*(A(IB+I)+A(ID+I))
357       I=I+INC3
358       J=J+INC4
359   470 CONTINUE
360       IBASE=IBASE+INC1
361       JBASE=JBASE+INC2
362   480 CONTINUE
363       GO TO 900
365   490 CONTINUE
366       Z=1.0/FLOAT(N)
367       DO 494 L=1,LA
368       I=IBASE
369       J=JBASE
370 !DIR$ IVDEP
371       DO 492 IJK=1,LOT
372       C(JA+J)=Z*((A(IA+I)+A(IC+I))+(A(IB+I)+A(ID+I)))
373       C(JC+J)=Z*((A(IA+I)+A(IC+I))-(A(IB+I)+A(ID+I)))
374       C(JB+J)=Z*(A(IA+I)-A(IC+I))
375       D(JB+J)=Z*(A(ID+I)-A(IB+I))
376       I=I+INC3
377       J=J+INC4
378   492 CONTINUE
379       IBASE=IBASE+INC1
380       JBASE=JBASE+INC2
381   494 CONTINUE
382       GO TO 900
384 !     CODING FOR FACTOR 5
385 !     -------------------
386   500 CONTINUE
387       IA=1
388       IB=IA+IINK
389       IC=IB+IINK
390       ID=IC+IINK
391       IE=ID+IINK
392       JA=1
393       JB=JA+(2*M-LA)*INC2
394       JC=JB+2*M*INC2
395       JD=JC
396       JE=JB
398       IF (LA.EQ.M) GO TO 590
400       DO 520 L=1,LA
401       I=IBASE
402       J=JBASE
403 !DIR$ IVDEP
404       DO 510 IJK=1,LOT
405       A1=A(IB+I)+A(IE+I)
406       A3=A(IB+I)-A(IE+I)
407       A2=A(IC+I)+A(ID+I)
408       A4=A(IC+I)-A(ID+I)
409       A5=A(IA+I)-0.25*(A1+A2)
410       A6=QRT5*(A1-A2)
411       C(JA+J)=A(IA+I)+(A1+A2)
412       C(JB+J)=A5+A6
413       C(JC+J)=A5-A6
414       D(JB+J)=-SIN72*A3-SIN36*A4
415       D(JC+J)=-SIN36*A3+SIN72*A4
416       I=I+INC3
417       J=J+INC4
418   510 CONTINUE
419       IBASE=IBASE+INC1
420       JBASE=JBASE+INC2
421   520 CONTINUE
422       JA=JA+JINK
423       JINK=2*JINK
424       JB=JB+JINK
425       JC=JC+JINK
426       JD=JD-JINK
427       JE=JE-JINK
428       IBASE=IBASE+IJUMP
429       IJUMP=2*IJUMP+IINK
430       IF (JB.EQ.JD) GO TO 560
431       DO 550 K=LA,KSTOP,LA
432       KB=K+K
433       KC=KB+KB
434       KD=KC+KB
435       KE=KD+KB
436       C1=TRIGS(KB+1)
437       S1=TRIGS(KB+2)
438       C2=TRIGS(KC+1)
439       S2=TRIGS(KC+2)
440       C3=TRIGS(KD+1)
441       S3=TRIGS(KD+2)
442       C4=TRIGS(KE+1)
443       S4=TRIGS(KE+2)
444       JBASE=0
445       DO 540 L=1,LA
446       I=IBASE
447       J=JBASE
448 !DIR$ IVDEP
449       DO 530 IJK=1,LOT
450       A1=(C1*A(IB+I)+S1*B(IB+I))+(C4*A(IE+I)+S4*B(IE+I))
451       A3=(C1*A(IB+I)+S1*B(IB+I))-(C4*A(IE+I)+S4*B(IE+I))
452       A2=(C2*A(IC+I)+S2*B(IC+I))+(C3*A(ID+I)+S3*B(ID+I))
453       A4=(C2*A(IC+I)+S2*B(IC+I))-(C3*A(ID+I)+S3*B(ID+I))
454       B1=(C1*B(IB+I)-S1*A(IB+I))+(C4*B(IE+I)-S4*A(IE+I))
455       B3=(C1*B(IB+I)-S1*A(IB+I))-(C4*B(IE+I)-S4*A(IE+I))
456       B2=(C2*B(IC+I)-S2*A(IC+I))+(C3*B(ID+I)-S3*A(ID+I))
457       B4=(C2*B(IC+I)-S2*A(IC+I))-(C3*B(ID+I)-S3*A(ID+I))
458       A5=A(IA+I)-0.25*(A1+A2)
459       A6=QRT5*(A1-A2)
460       B5=B(IA+I)-0.25*(B1+B2)
461       B6=QRT5*(B1-B2)
462       A10=A5+A6
463       A20=A5-A6
464       B10=B5+B6
465       B20=B5-B6
466       A11=SIN72*B3+SIN36*B4
467       A21=SIN36*B3-SIN72*B4
468       B11=SIN72*A3+SIN36*A4
469       B21=SIN36*A3-SIN72*A4
470       C(JA+J)=A(IA+I)+(A1+A2)
471       C(JB+J)=A10+A11
472       C(JE+J)=A10-A11
473       C(JC+J)=A20+A21
474       C(JD+J)=A20-A21
475       D(JA+J)=B(IA+I)+(B1+B2)
476       D(JB+J)=B10-B11
477       D(JE+J)=-(B10+B11)
478       D(JC+J)=B20-B21
479       D(JD+J)=-(B20+B21)
480       I=I+INC3
481       J=J+INC4
482   530 CONTINUE
483       IBASE=IBASE+INC1
484       JBASE=JBASE+INC2
485   540 CONTINUE
486       IBASE=IBASE+IJUMP
487       JA=JA+JINK
488       JB=JB+JINK
489       JC=JC+JINK
490       JD=JD-JINK
491       JE=JE-JINK
492   550 CONTINUE
493       IF (JB.GT.JD) GO TO 900
494   560 CONTINUE
495       JBASE=0
496       DO 580 L=1,LA
497       I=IBASE
498       J=JBASE
499 !DIR$ IVDEP
500       DO 570 IJK=1,LOT
501       A1=A(IB+I)+A(IE+I)
502       A3=A(IB+I)-A(IE+I)
503       A2=A(IC+I)+A(ID+I)
504       A4=A(IC+I)-A(ID+I)
505       A5=A(IA+I)+0.25*(A3-A4)
506       A6=QRT5*(A3+A4)
507       C(JA+J)=A5+A6
508       C(JB+J)=A5-A6
509       C(JC+J)=A(IA+I)-(A3-A4)
510       D(JA+J)=-SIN36*A1-SIN72*A2
511       D(JB+J)=-SIN72*A1+SIN36*A2
512       I=I+INC3
513       J=J+INC4
514   570 CONTINUE
515       IBASE=IBASE+INC1
516       JBASE=JBASE+INC2
517   580 CONTINUE
518       GO TO 900
520   590 CONTINUE
521       Z=1.0/FLOAT(N)
522       ZQRT5=Z*QRT5
523       ZSIN36=Z*SIN36
524       ZSIN72=Z*SIN72
525       DO 594 L=1,LA
526       I=IBASE
527       J=JBASE
528 !DIR$ IVDEP
529       DO 592 IJK=1,LOT
530       A1=A(IB+I)+A(IE+I)
531       A3=A(IB+I)-A(IE+I)
532       A2=A(IC+I)+A(ID+I)
533       A4=A(IC+I)-A(ID+I)
534       A5=Z*(A(IA+I)-0.25*(A1+A2))
535       A6=ZQRT5*(A1-A2)
536       C(JA+J)=Z*(A(IA+I)+(A1+A2))
537       C(JB+J)=A5+A6
538       C(JC+J)=A5-A6
539       D(JB+J)=-ZSIN72*A3-ZSIN36*A4
540       D(JC+J)=-ZSIN36*A3+ZSIN72*A4
541       I=I+INC3
542       J=J+INC4
543   592 CONTINUE
544       IBASE=IBASE+INC1
545       JBASE=JBASE+INC2
546   594 CONTINUE
547       GO TO 900
549 !     CODING FOR FACTOR 6
550 !     -------------------
551   600 CONTINUE
552       IA=1
553       IB=IA+IINK
554       IC=IB+IINK
555       ID=IC+IINK
556       IE=ID+IINK
557       IF=IE+IINK
558       JA=1
559       JB=JA+(2*M-LA)*INC2
560       JC=JB+2*M*INC2
561       JD=JC+2*M*INC2
562       JE=JC
563       JF=JB
565       IF (LA.EQ.M) GO TO 690
567       DO 620 L=1,LA
568       I=IBASE
569       J=JBASE
570 !DIR$ IVDEP
571       DO 610 IJK=1,LOT
572       A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I))
573       C(JA+J)=(A(IA+I)+A(ID+I))+A11
574       C(JC+J)=(A(IA+I)+A(ID+I)-0.5*A11)
575       D(JC+J)=SIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I)))
576       A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I))
577       C(JB+J)=(A(IA+I)-A(ID+I))-0.5*A11
578       D(JB+J)=SIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I)))
579       C(JD+J)=(A(IA+I)-A(ID+I))+A11
580       I=I+INC3
581       J=J+INC4
582   610 CONTINUE
583       IBASE=IBASE+INC1
584       JBASE=JBASE+INC2
585   620 CONTINUE
586       JA=JA+JINK
587       JINK=2*JINK
588       JB=JB+JINK
589       JC=JC+JINK
590       JD=JD-JINK
591       JE=JE-JINK
592       JF=JF-JINK
593       IBASE=IBASE+IJUMP
594       IJUMP=2*IJUMP+IINK
595       IF (JC.EQ.JD) GO TO 660
596       DO 650 K=LA,KSTOP,LA
597       KB=K+K
598       KC=KB+KB
599       KD=KC+KB
600       KE=KD+KB
601       KF=KE+KB
602       C1=TRIGS(KB+1)
603       S1=TRIGS(KB+2)
604       C2=TRIGS(KC+1)
605       S2=TRIGS(KC+2)
606       C3=TRIGS(KD+1)
607       S3=TRIGS(KD+2)
608       C4=TRIGS(KE+1)
609       S4=TRIGS(KE+2)
610       C5=TRIGS(KF+1)
611       S5=TRIGS(KF+2)
612       JBASE=0
613       DO 640 L=1,LA
614       I=IBASE
615       J=JBASE
616 !DIR$ IVDEP
617       DO 630 IJK=1,LOT
618       A1=C1*A(IB+I)+S1*B(IB+I)
619       B1=C1*B(IB+I)-S1*A(IB+I)
620       A2=C2*A(IC+I)+S2*B(IC+I)
621       B2=C2*B(IC+I)-S2*A(IC+I)
622       A3=C3*A(ID+I)+S3*B(ID+I)
623       B3=C3*B(ID+I)-S3*A(ID+I)
624       A4=C4*A(IE+I)+S4*B(IE+I)
625       B4=C4*B(IE+I)-S4*A(IE+I)
626       A5=C5*A(IF+I)+S5*B(IF+I)
627       B5=C5*B(IF+I)-S5*A(IF+I)
628       A11=(A2+A5)+(A1+A4)
629       A20=(A(IA+I)+A3)-0.5*A11
630       A21=SIN60*((A2+A5)-(A1+A4))
631       B11=(B2+B5)+(B1+B4)
632       B20=(B(IA+I)+B3)-0.5*B11
633       B21=SIN60*((B2+B5)-(B1+B4))
634       C(JA+J)=(A(IA+I)+A3)+A11
635       D(JA+J)=(B(IA+I)+B3)+B11
636       C(JC+J)=A20-B21
637       D(JC+J)=A21+B20
638       C(JE+J)=A20+B21
639       D(JE+J)=A21-B20
640       A11=(A2-A5)+(A4-A1)
641       A20=(A(IA+I)-A3)-0.5*A11
642       A21=SIN60*((A4-A1)-(A2-A5))
643       B11=(B5-B2)-(B4-B1)
644       B20=(B3-B(IA+I))-0.5*B11
645       B21=SIN60*((B5-B2)+(B4-B1))
646       C(JB+J)=A20-B21
647       D(JB+J)=A21-B20
648       C(JD+J)=A11+(A(IA+I)-A3)
649       D(JD+J)=B11+(B3-B(IA+I))
650       C(JF+J)=A20+B21
651       D(JF+J)=A21+B20
652       I=I+INC3
653       J=J+INC4
654   630 CONTINUE
655       IBASE=IBASE+INC1
656       JBASE=JBASE+INC2
657   640 CONTINUE
658       IBASE=IBASE+IJUMP
659       JA=JA+JINK
660       JB=JB+JINK
661       JC=JC+JINK
662       JD=JD-JINK
663       JE=JE-JINK
664       JF=JF-JINK
665   650 CONTINUE
666       IF (JC.GT.JD) GO TO 900
667   660 CONTINUE
668       JBASE=0
669       DO 680 L=1,LA
670       I=IBASE
671       J=JBASE
672 !DIR$ IVDEP
673       DO 670 IJK=1,LOT
674       C(JA+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))+ SIN60*(A(IB+I)-A(IF+I))
675       D(JA+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))-SIN60*(A(IC+I)+A(IE+I))
676       C(JB+J)=A(IA+I)-(A(IC+I)-A(IE+I))
677       D(JB+J)=A(ID+I)-(A(IB+I)+A(IF+I))
678       C(JC+J)=(A(IA+I)+0.5*(A(IC+I)-A(IE+I)))-SIN60*(A(IB+I)-A(IF+I))
679       D(JC+J)=-(A(ID+I)+0.5*(A(IB+I)+A(IF+I)))+SIN60*(A(IC+I)+A(IE+I))
680       I=I+INC3
681       J=J+INC4
682   670 CONTINUE
683       IBASE=IBASE+INC1
684       JBASE=JBASE+INC2
685   680 CONTINUE
686       GO TO 900
688   690 CONTINUE
689       Z=1.0/FLOAT(N)
690       ZSIN60=Z*SIN60
691       DO 694 L=1,LA
692       I=IBASE
693       J=JBASE
694 !DIR$ IVDEP
695       DO 692 IJK=1,LOT
696       A11=(A(IC+I)+A(IF+I))+(A(IB+I)+A(IE+I))
697       C(JA+J)=Z*((A(IA+I)+A(ID+I))+A11)
698       C(JC+J)=Z*((A(IA+I)+A(ID+I))-0.5*A11)
699       D(JC+J)=ZSIN60*((A(IC+I)+A(IF+I))-(A(IB+I)+A(IE+I)))
700       A11=(A(IC+I)-A(IF+I))+(A(IE+I)-A(IB+I))
701       C(JB+J)=Z*((A(IA+I)-A(ID+I))-0.5*A11)
702       D(JB+J)=ZSIN60*((A(IE+I)-A(IB+I))-(A(IC+I)-A(IF+I)))
703       C(JD+J)=Z*((A(IA+I)-A(ID+I))+A11)
704       I=I+INC3
705       J=J+INC4
706   692 CONTINUE
707       IBASE=IBASE+INC1
708       JBASE=JBASE+INC2
709   694 CONTINUE
710       GO TO 900
712 !     CODING FOR FACTOR 8
713 !     -------------------
714   800 CONTINUE
715       IBAD=3
716       IF (LA.NE.M) GO TO 910
717       IA=1
718       IB=IA+IINK
719       IC=IB+IINK
720       ID=IC+IINK
721       IE=ID+IINK
722       IF=IE+IINK
723       IG=IF+IINK
724       IH=IG+IINK
725       JA=1
726       JB=JA+LA*INC2
727       JC=JB+2*M*INC2
728       JD=JC+2*M*INC2
729       JE=JD+2*M*INC2
730       Z=1.0/FLOAT(N)
731       ZSIN45=Z*SQRT(0.5)
733       DO 820 L=1,LA
734       I=IBASE
735       J=JBASE
736 !DIR$ IVDEP
737       DO 810 IJK=1,LOT
738       C(JA+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))+ &
739          ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I))))
740       C(JE+J)=Z*(((A(IA+I)+A(IE+I))+(A(IC+I)+A(IG+I)))- &
741          ((A(ID+I)+A(IH+I))+(A(IB+I)+A(IF+I))))
742       C(JC+J)=Z*((A(IA+I)+A(IE+I))-(A(IC+I)+A(IG+I)))
743       D(JC+J)=Z*((A(ID+I)+A(IH+I))-(A(IB+I)+A(IF+I)))
744       C(JB+J)=Z*(A(IA+I)-A(IE+I))+ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I)))
745       C(JD+J)=Z*(A(IA+I)-A(IE+I))-ZSIN45*((A(IH+I)-A(ID+I))-(A(IF+I)-A(IB+I)))
746       D(JB+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))+Z*(A(IG+I)-A(IC+I))
747       D(JD+J)=ZSIN45*((A(IH+I)-A(ID+I))+(A(IF+I)-A(IB+I)))-Z*(A(IG+I)-A(IC+I))
748       I=I+INC3
749       J=J+INC4
750   810 CONTINUE
751       IBASE=IBASE+INC1
752       JBASE=JBASE+INC2
753   820 CONTINUE
756   900 CONTINUE
757       IBAD=0
758   910 CONTINUE
759       IERR=IBAD
761       RETURN
762       END SUBROUTINE QPASSM