1 !C SUBROUTINE 'QPASSM' - PERFORMS ONE PASS THROUGH DATA AS PART!C OF MULTIPLE REAL FFT (FOURIER ANALYSIS) ROUTINE
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 ---------------------------------------------------------------
55 KSTOP=(N-IFAC)/(2*IFAC)
59 IF (LOT.GT.VectorLength) GO TO 910
65 IF (IGO.LT.1.OR.IGO.GT.6) GO TO 910
66 GO TO (200,300,400,500,600,800),IGO
76 IF (LA.EQ.M) GO TO 290
83 C(JA+J)=A(IA+I)+A(IB+I)
84 C(JB+J)=A(IA+I)-A(IB+I)
96 IF (JA.EQ.JB) GO TO 260
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)
121 IF (JA.GT.JB) GO TO 900
146 C(JA+J)=Z*(A(IA+I)+A(IB+I))
147 C(JB+J)=Z*(A(IA+I)-A(IB+I))
156 ! CODING FOR FACTOR 3
157 ! -------------------
166 IF (LA.EQ.M) GO TO 390
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))
188 IF (JA.EQ.JC) GO TO 360
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))
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)))
225 IF (JA.GT.JC) GO TO 900
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))
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))
263 ! CODING FOR FACTOR 4
264 ! -------------------
275 IF (LA.EQ.M) GO TO 490
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)
299 IF (JB.EQ.JC) GO TO 460
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))
344 IF (JB.GT.JC) GO TO 900
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))
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))
384 ! CODING FOR FACTOR 5
385 ! -------------------
398 IF (LA.EQ.M) GO TO 590
409 A5=A(IA+I)-0.25*(A1+A2)
411 C(JA+J)=A(IA+I)+(A1+A2)
414 D(JB+J)=-SIN72*A3-SIN36*A4
415 D(JC+J)=-SIN36*A3+SIN72*A4
430 IF (JB.EQ.JD) GO TO 560
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)
460 B5=B(IA+I)-0.25*(B1+B2)
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)
475 D(JA+J)=B(IA+I)+(B1+B2)
493 IF (JB.GT.JD) GO TO 900
505 A5=A(IA+I)+0.25*(A3-A4)
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
534 A5=Z*(A(IA+I)-0.25*(A1+A2))
536 C(JA+J)=Z*(A(IA+I)+(A1+A2))
539 D(JB+J)=-ZSIN72*A3-ZSIN36*A4
540 D(JC+J)=-ZSIN36*A3+ZSIN72*A4
549 ! CODING FOR FACTOR 6
550 ! -------------------
565 IF (LA.EQ.M) GO TO 690
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
595 IF (JC.EQ.JD) GO TO 660
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)
629 A20=(A(IA+I)+A3)-0.5*A11
630 A21=SIN60*((A2+A5)-(A1+A4))
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
641 A20=(A(IA+I)-A3)-0.5*A11
642 A21=SIN60*((A4-A1)-(A2-A5))
644 B20=(B3-B(IA+I))-0.5*B11
645 B21=SIN60*((B5-B2)+(B4-B1))
648 C(JD+J)=A11+(A(IA+I)-A3)
649 D(JD+J)=B11+(B3-B(IA+I))
666 IF (JC.GT.JD) GO TO 900
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))
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)
712 ! CODING FOR FACTOR 8
713 ! -------------------
716 IF (LA.NE.M) GO TO 910
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))
762 END SUBROUTINE QPASSM