Mixed Radix Fast Fourier Transforms
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: MXFFT
C AUTHOR: RICHARD C. SINGLETON
C SRI INTERNATIONAL, MENLO PARK, CALIFORNIA 94025
C INPUT: NONE
C-----------------------------------------------------------------------
C
DIMENSION A(4097), B(4097), C(4097), D(4097), NC(17)
DIMENSION IA(64), IB(64), RA(64), RB(64)
COMMON /CSTAK/ DSTAK(2500)
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (IA(1),RA(1))
EQUIVALENCE (IB(1),RB(1))
EQUIVALENCE (ISTAK(1),LOUT)
EQUIVALENCE (ISTAK(3),LUSED)
DATA NC(1), NC(2), NC(3), NC(4), NC(5), NC(6), NC(7), NC(8),
* NC(9), NC(10), NC(11), NC(12), NC(13), NC(14), NC(15),
* NC(16), NC(17) /2,3,5,210,1000,2000,1024,2048,4096,2187,3125,
* 2401,1331,2197,289,361,529/
C
C SET UP MACHINE CONSTANTS, USING PORT SUBPROGRAMS
C
IOUTD = I1MACH(2)
C
N = 16
DO 10 J=1,N
C(J) = UNI(K)
D(J) = UNI(K)
10 CONTINUE
9999 FORMAT (/5X, 1HJ, 8X, 4HA(J), 12X, 4HB(J), 10X, 6HA(J+1), 10X,
* 6HB(J+1))
9998 FORMAT (1X, I5, 4E16.8)
C
DO 20 J=1,N
A(J) = C(J)
B(J) = D(J)
20 CONTINUE
WRITE (IOUTD,9997)
9997 FORMAT (///12X, 44HTEST OF FFT, REALS AND REALT FOR REAL VALUES)
WRITE (IOUTD,9996)
9996 FORMAT (/24H REAL INPUT SEQUENCE ---)
WRITE (IOUTD,9999)
WRITE(IOUTD,9998) (J,A(J),B(J),A(J+1),B(J+1),J=1,N,2)
C
CALL FFT(A, B, 1, N, 1, -1)
CALL REALS(A, B, N, -1)
WRITE (IOUTD,9995)
9995 FORMAT (/41H FOURIER COSINE AND SINE COEFFICIENTS ---)
WRITE (IOUTD,9999)
WRITE(IOUTD,9998) (J,A(J),B(J),A(J+1),B(J+1),J=1,N,2)
J = N + 1
WRITE (IOUTD,9998) J, A(J), B(J)
C
C THE NEXT CALL ON REALT DOES THE SAME THING AS: CALL REALS(A,B,N,1)
C
CALL REALT(A, B, 1, N, 1, 1)
CALL FFT(A, B, 1, N, 1, 1)
WRITE (IOUTD,9994)
9994 FORMAT (/29H TRANSFORM-INVERSE RESULT ---)
WRITE (IOUTD,9999)
WRITE(IOUTD,9998) (J,A(J),B(J),A(J+1),B(J+1),J=1,N,2)
C
DO 30 J=1,N
IA(J) = J
IB(J) = J
A(J) = A(J) - C(J)
B(J) = B(J) - D(J)
30 CONTINUE
CALL SORTG(A, N, RA)
CALL SORTG(B, N, RB)
WRITE (IOUTD,9993)
9993 FORMAT (/24H SORTED ERROR VALUES ---//9X, 1HJ, 11X, 4HA(J), 14X,
* 1HJ, 11X, 4HB(J))
WRITE (IOUTD,9992) (IA(J),A(J),IB(J),B(J),J=1,N)
9992 FORMAT (1X, I9, E20.8, I10, E20.8)
C
DO 40 J=1,N
A(2*J-1) = C(J)
A(2*J) = D(J)
B(2*J-1) = 0.0
B(2*J) = 0.0
40 CONTINUE
WRITE (IOUTD,9991)
9991 FORMAT (///16X, 27HTEST OF FFT FOR REAL VALUES)
C
N = N + N
CALL FFT(A, B, 1, N, 1, -1)
WRITE (IOUTD,9995)
WRITE (IOUTD,9999)
WRITE(IOUTD,9998) (J,A(J),B(J),A(J+1),B(J+1),J=1,N,2)
C
CALL FFT(A, B, 1, N, 1, 1)
WRITE (IOUTD,9994)
WRITE (IOUTD,9999)
WRITE(IOUTD,9998) (J,A(J),B(J),A(J+1),B(J+1),J=1,N,2)
C
C GENERATE TEST DATA FOR TRANSFORM OF MAXIMUM SIZE 4096
C
DO 50 J=1,4096
CALL NORMAL(C(J), D(J))
50 CONTINUE
C
WRITE (IOUTD,9990)
9990 FORMAT (///15X, 35HSINGLE-VARIATE TESTS OF SUBROUTINES)
WRITE (IOUTD,9989)
9989 FORMAT (/40X, 9HRMS ERROR/25X, 1HN, 9X, 7HA ARRAY, 5X, 7HB ARRAY)
C
DO 90 I=1,17
N = NC(I)
WRITE (IOUTD,9998)
C
DO 60 J=1,N
A(J) = C(J)
B(J) = D(J)
60 CONTINUE
CALL FFT(A, B, 1, N, 1, -1)
CALL FFT(A, B, 1, N, 1, 1)
CALL RMS(A, B, C, D, N, SS1, SS2)
WRITE (IOUTD,9988) N, SS1, SS2
9988 FORMAT (10X, 3HFFT, I13, 4X, 2E12.3)
C
DO 70 J=1,N
A(J) = C(J)
B(J) = D(J)
70 CONTINUE
CALL REALS(A, B, N, -1)
CALL REALS(A, B, N, 1)
CALL RMS(A, B, C, D, N, SS1, SS2)
WRITE (IOUTD,9987) N, SS1, SS2
9987 FORMAT (10X, 5HREALS, I11, 4X, 2E12.3)
C
DO 80 J=1,N
A(J) = C(J)
B(J) = D(J)
80 CONTINUE
CALL REALT(A, B, 1, N, 1, -1)
CALL REALT(A, B, 1, N, 1, 1)
CALL RMS(A, B, C, D, N, SS1, SS2)
WRITE (IOUTD,9986) N, SS1, SS2
9986 FORMAT (10X, 5HREALT, I11, 4X, 2E12.3)
C
90 CONTINUE
C
WRITE (IOUTD,9985)
9985 FORMAT (///15X, 34HMULTI-VARIATE TESTS OF SUBROUTINES)
WRITE (IOUTD,9989)
C
N1 = 4
N2 = 30
N3 = 9
N = N1*N2*N3
C
DO 100 J=1,N
A(J) = C(J)
B(J) = D(J)
100 CONTINUE
CALL FFT(A, B, N2*N3, N1, 1, -1)
CALL FFT(A, B, N3, N2, N1, -1)
CALL FFT(A, B, 1, N3, N1*N2, -1)
CALL FFT(A, B, 1, N3, N1*N2, 1)
CALL FFT(A, B, N3, N2, N1, 1)
CALL FFT(A, B, N2*N3, N1, 1, 1)
CALL RMS(A, B, C, D, N, SS1, SS2)
WRITE (IOUTD,9988) N, SS1, SS2
C
M = N + N
DO 110 J=1,M
A(J) = C(J)
B(J) = D(J)
110 CONTINUE
CALL FFT(A, A(2), N3, N2, N1, -2)
CALL FFT(A, A(2), N3, N2, N1, 2)
CALL RMS(A, B, C, D, M, SS1, SS2)
WRITE (IOUTD,9988) N, SS1, SS2
C
DO 120 J=1,N
A(J) = C(J)
B(J) = D(J)
120 CONTINUE
CALL REALT(A, B, N3, N2, N1, -1)
CALL REALT(A, B, N3, N2, N1, 1)
CALL RMS(A, B, C, D, N, SS1, SS2)
WRITE (IOUTD,9986) N, SS1, SS2
C
DO 130 J=1,M
A(J) = C(J)
B(J) = D(J)
130 CONTINUE
CALL REALT(A, A(2), N3, N2, N1, -2)
CALL REALT(A, A(2), N3, N2, N1, 2)
CALL RMS(A, B, C, D, M, SS1, SS2)
WRITE (IOUTD,9986) N, SS1, SS2
C
WRITE (IOUTD,9984)
9984 FORMAT (///15X, 46H8 BY 9 REAL TRANSFORM --- FFT VS. FFT AND REAL,
* 1HT)
WRITE (IOUTD,9983)
9983 FORMAT (/5X, 1HJ, 8X, 4HA(J), 12X, 4HB(J), 12X, 4HC(J), 12X,
* 4HD(J))
C
N1 = 4
N2 = 9
N = N1*N2
C
DO 140 J=1,N
A(2*J-1) = C(J)
B(2*J-1) = 0.0
A(2*J) = D(J)
B(2*J) = 0.0
140 CONTINUE
CALL FFT(C, D, N2, N1, 1, -1)
CALL REALT(C, D, N2, N1, 1, -1)
CALL FFT(C, D, 1, N2, N1, -1)
CALL FFT(C(N+1), D(N+1), 1, N2, 1, -1)
N = N + N
N1 = N1 + N1
CALL FFT(A, B, N2, N1, 1, -1)
CALL FFT(A, B, 1, N2, N1, -1)
C
N1 = (N1/2+1)*N2
N2 = N1 + 1
WRITE (IOUTD,9998) (K,A(K),B(K),C(K),D(K),K=1,N1)
WRITE (IOUTD,9982) (K,A(K),B(K),K=N2,N)
9982 FORMAT (1X, I5, 2E16.8)
C
WRITE (IOUTD,9981) LUSED
9981 FORMAT (/21H MAXIMUM STACK SIZE =, I6)
C
C DE-ALLOCATE WORKING STORAGE...AS A FINAL CHECK ON ARRAY BOUNDS
C
J = LOUT
IF (J.NE.0) CALL ISTKRL(J)
C
STOP
END
C
C-----------------------------------------------------------------------
C BLOCK DATA: INITIALIZES LABELED COMMON
C-----------------------------------------------------------------------
C
BLOCK DATA
C
COMMON /CSTAK/ DSTAK(2500)
C
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
INTEGER ISIZE(5)
C
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (ISTAK(1),LOUT)
EQUIVALENCE (ISTAK(2),LNOW)
EQUIVALENCE (ISTAK(3),LUSED)
EQUIVALENCE (ISTAK(4),LMAX)
EQUIVALENCE (ISTAK(5),LBOOK)
EQUIVALENCE (ISTAK(6),ISIZE(1))
C
DATA ISIZE(1), ISIZE(2), ISIZE(3), ISIZE(4), ISIZE(5) /1,1,1,2,2/
DATA LOUT, LNOW, LUSED, LMAX, LBOOK /0,10,10,5000,10/
C
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFT
C MULTIVARIATE COMPLEX FOURIER TRANSFORM, COMPUTED IN PLACE
C USING MIXED-RADIX FAST FOURIER TRANSFORM ALGORITHM.
C-----------------------------------------------------------------------
C
SUBROUTINE FFT(A, B, NSEG, N, NSPN, ISN)
C
C ARRAYS A AND B ORIGINALLY HOLD THE REAL AND IMAGINARY
C COMPONENTS OF THE DATA, AND RETURN THE REAL AND
C IMAGINARY COMPONENTS OF THE RESULTING FOURIER COEFFICIENTS.
C MULTIVARIATE DATA IS INDEXED ACCORDING TO THE FORTRAN
C ARRAY ELEMENT SUCCESSOR FUNCTION, WITHOUT LIMIT
C ON THE NUMBER OF IMPLIED MULTIPLE SUBSCRIPTS.
C THE SUBROUTINE IS CALLED ONCE FOR EACH VARIATE.
C THE CALLS FOR A MULTIVARIATE TRANSFORM MAY BE IN ANY ORDER.
C
C N IS THE DIMENSION OF THE CURRENT VARIABLE.
C NSPN IS THE SPACING OF CONSECUTIVE DATA VALUES
C WHILE INDEXING THE CURRENT VARIABLE.
C NSEG*N*NSPN IS THE TOTAL NUMBER OF COMPLEX DATA VALUES.
C THE SIGN OF ISN DETERMINES THE SIGN OF THE COMPLEX
C EXPONENTIAL, AND THE MAGNITUDE OF ISN IS NORMALLY ONE.
C THE MAGNITUDE OF ISN DETERMINES THE INDEXING INCREMENT FOR A&B.
C
C IF FFT IS CALLED TWICE, WITH OPPOSITE SIGNS ON ISN, AN
C IDENTITY TRANSFORMATION IS DONE...CALLS CAN BE IN EITHER ORDER.
C THE RESULTS ARE SCALED BY 1/N WHEN THE SIGN OF ISN IS POSITIVE.
C
C A TRI-VARIATE TRANSFORM WITH A(N1,N2,N3), B(N1,N2,N3)
C IS COMPUTED BY
C CALL FFT(A,B,N2*N3,N1,1,-1)
C CALL FFT(A,B,N3,N2,N1,-1)
C CALL FFT(A,B,1,N3,N1*N2,-1)
C
C A SINGLE-VARIATE TRANSFORM OF N COMPLEX DATA VALUES IS COMPUTED BY
C CALL FFT(A,B,1,N,1,-1)
C
C THE DATA MAY ALTERNATIVELY BE STORED IN A SINGLE COMPLEX
C ARRAY A, THEN THE MAGNITUDE OF ISN CHANGED TO TWO TO
C GIVE THE CORRECT INDEXING INCREMENT AND A(2) USED TO
C PASS THE INITIAL ADDRESS FOR THE SEQUENCE OF IMAGINARY
C VALUES, E.G.
C CALL FFT(A,A(2),NSEG,N,NSPN,-2)
C
C ARRAY NFAC IS WORKING STORAGE FOR FACTORING N. THE SMALLEST
C NUMBER EXCEEDING THE 15 LOCATIONS PROVIDED IS 12,754,584.
C
DIMENSION A(1), B(1), NFAC(15)
C
COMMON /CSTAK/ DSTAK(2500)
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
REAL RSTAK(5000)
C
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (DSTAK(1),RSTAK(1))
C
C DETERMINE THE FACTORS OF N
C
M = 0
NF = IABS(N)
K = NF
IF (NF.EQ.1) RETURN
NSPAN = IABS(NF*NSPN)
NTOT = IABS(NSPAN*NSEG)
IF (ISN*NTOT.NE.0) GO TO 20
IERR = I1MACH(4)
WRITE (IERR,9999) NSEG, N, NSPN, ISN
9999 FORMAT (31H ERROR - ZERO IN FFT PARAMETERS, 4I10)
RETURN
C
10 M = M + 1
NFAC(M) = 4
K = K/16
20 IF (K-(K/16)*16.EQ.0) GO TO 10
J = 3
JJ = 9
GO TO 40
30 M = M + 1
NFAC(M) = J
K = K/JJ
40 IF (MOD(K,JJ).EQ.0) GO TO 30
J = J + 2
JJ = J**2
IF (JJ.LE.K) GO TO 40
IF (K.GT.4) GO TO 50
KT = M
NFAC(M+1) = K
IF (K.NE.1) M = M + 1
GO TO 90
50 IF (K-(K/4)*4.NE.0) GO TO 60
M = M + 1
NFAC(M) = 2
K = K/4
C ALL SQUARE FACTORS OUT NOW, BUT K .GE. 5 STILL
60 KT = M
MAXP = MAX0(KT+KT+2,K-1)
J = 2
70 IF (MOD(K,J).NE.0) GO TO 80
M = M + 1
NFAC(M) = J
K = K/J
80 J = ((J+1)/2)*2 + 1
IF (J.LE.K) GO TO 70
90 IF (M.LE.KT+1) MAXP = M + KT + 1
IF (M+KT.GT.15) GO TO 120
IF (KT.EQ.0) GO TO 110
J = KT
100 M = M + 1
NFAC(M) = NFAC(J)
J = J - 1
IF (J.NE.0) GO TO 100
C
110 MAXF = M - KT
MAXF = NFAC(MAXF)
IF (KT.GT.0) MAXF = MAX0(NFAC(KT),MAXF)
J = ISTKGT(MAXF*4,3)
JJ = J + MAXF
J2 = JJ + MAXF
J3 = J2 + MAXF
K = ISTKGT(MAXP,2)
CALL FFTMX(A, B, NTOT, NF, NSPAN, ISN, M, KT, RSTAK(J),
* RSTAK(JJ), RSTAK(J2), RSTAK(J3), ISTAK(K), NFAC)
CALL ISTKRL(2)
RETURN
C
120 IERR = I1MACH(4)
WRITE (IERR,9998) N
9998 FORMAT (50H ERROR - FFT PARAMETER N HAS MORE THAN 15 FACTORS-,
* I20)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FFTMX
C CALLED BY SUBROUTINE 'FFT' TO COMPUTE MIXED-RADIX FOURIER TRANSFORM
C-----------------------------------------------------------------------
C
SUBROUTINE FFTMX(A, B, NTOT, N, NSPAN, ISN, M, KT, AT, CK, BT,
* SK, NP, NFAC)
C
DIMENSION A(1), B(1), AT(1), CK(1), BT(1), SK(1), NP(1), NFAC(1)
C
INC = IABS(ISN)
NT = INC*NTOT
KS = INC*NSPAN
RAD = ATAN(1.0)
S72 = RAD/0.625
C72 = COS(S72)
S72 = SIN(S72)
S120 = SQRT(0.75)
IF (ISN.GT.0) GO TO 10
S72 = -S72
S120 = -S120
RAD = -RAD
GO TO 30
C
C SCALE BY 1/N FOR ISN .GT. 0
C
10 AK = 1.0/FLOAT(N)
DO 20 J=1,NT,INC
A(J) = A(J)*AK
B(J) = B(J)*AK
20 CONTINUE
C
30 KSPAN = KS
NN = NT - INC
JC = KS/N
C
C SIN, COS VALUES ARE RE-INITIALIZED EACH LIM STEPS
C
LIM = 32
KLIM = LIM*JC
I = 0
JF = 0
MAXF = M - KT
MAXF = NFAC(MAXF)
IF (KT.GT.0) MAXF = MAX0(NFAC(KT),MAXF)
C
C COMPUTE FOURIER TRANSFORM
C
40 DR = 8.0*FLOAT(JC)/FLOAT(KSPAN)
CD = 2.0*SIN(0.5*DR*RAD)**2
SD = SIN(DR*RAD)
KK = 1
I = I + 1
IF (NFAC(I).NE.2) GO TO 110
C
C TRANSFORM FOR FACTOR OF 2 (INCLUDING ROTATION FACTOR)
C
KSPAN = KSPAN/2
K1 = KSPAN + 2
50 K2 = KK + KSPAN
AK = A(K2)
BK = B(K2)
A(K2) = A(KK) - AK
B(K2) = B(KK) - BK
A(KK) = A(KK) + AK
B(KK) = B(KK) + BK
KK = K2 + KSPAN
IF (KK.LE.NN) GO TO 50
KK = KK - NN
IF (KK.LE.JC) GO TO 50
IF (KK.GT.KSPAN) GO TO 350
60 C1 = 1.0 - CD
S1 = SD
MM = MIN0(K1/2,KLIM)
GO TO 80
70 AK = C1 - (CD*C1+SD*S1)
S1 = (SD*C1-CD*S1) + S1
C
C THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C ERROR. IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE
C C1=AK
C
C1 = 0.5/(AK**2+S1**2) + 0.5
S1 = C1*S1
C1 = C1*AK
80 K2 = KK + KSPAN
AK = A(KK) - A(K2)
BK = B(KK) - B(K2)
A(KK) = A(KK) + A(K2)
B(KK) = B(KK) + B(K2)
A(K2) = C1*AK - S1*BK
B(K2) = S1*AK + C1*BK
KK = K2 + KSPAN
IF (KK.LT.NT) GO TO 80
K2 = KK - NT
C1 = -C1
KK = K1 - K2
IF (KK.GT.K2) GO TO 80
KK = KK + JC
IF (KK.LE.MM) GO TO 70
IF (KK.LT.K2) GO TO 90
K1 = K1 + INC + INC
KK = (K1-KSPAN)/2 + JC
IF (KK.LE.JC+JC) GO TO 60
GO TO 40
90 S1 = FLOAT((KK-1)/JC)*DR*RAD
C1 = COS(S1)
S1 = SIN(S1)
MM = MIN0(K1/2,MM+KLIM)
GO TO 80
C
C TRANSFORM FOR FACTOR OF 3 (OPTIONAL CODE)
C
100 K1 = KK + KSPAN
K2 = K1 + KSPAN
AK = A(KK)
BK = B(KK)
AJ = A(K1) + A(K2)
BJ = B(K1) + B(K2)
A(KK) = AK + AJ
B(KK) = BK + BJ
AK = -0.5*AJ + AK
BK = -0.5*BJ + BK
AJ = (A(K1)-A(K2))*S120
BJ = (B(K1)-B(K2))*S120
A(K1) = AK - BJ
B(K1) = BK + AJ
A(K2) = AK + BJ
B(K2) = BK - AJ
KK = K2 + KSPAN
IF (KK.LT.NN) GO TO 100
KK = KK - NN
IF (KK.LE.KSPAN) GO TO 100
GO TO 290
C
C TRANSFORM FOR FACTOR OF 4
C
110 IF (NFAC(I).NE.4) GO TO 230
KSPNN = KSPAN
KSPAN = KSPAN/4
120 C1 = 1.0
S1 = 0
MM = MIN0(KSPAN,KLIM)
GO TO 150
130 C2 = C1 - (CD*C1+SD*S1)
S1 = (SD*C1-CD*S1) + S1
C
C THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C ERROR. IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE
C C1=C2
C
C1 = 0.5/(C2**2+S1**2) + 0.5
S1 = C1*S1
C1 = C1*C2
140 C2 = C1**2 - S1**2
S2 = C1*S1*2.0
C3 = C2*C1 - S2*S1
S3 = C2*S1 + S2*C1
150 K1 = KK + KSPAN
K2 = K1 + KSPAN
K3 = K2 + KSPAN
AKP = A(KK) + A(K2)
AKM = A(KK) - A(K2)
AJP = A(K1) + A(K3)
AJM = A(K1) - A(K3)
A(KK) = AKP + AJP
AJP = AKP - AJP
BKP = B(KK) + B(K2)
BKM = B(KK) - B(K2)
BJP = B(K1) + B(K3)
BJM = B(K1) - B(K3)
B(KK) = BKP + BJP
BJP = BKP - BJP
IF (ISN.LT.0) GO TO 180
AKP = AKM - BJM
AKM = AKM + BJM
BKP = BKM + AJM
BKM = BKM - AJM
IF (S1.EQ.0.0) GO TO 190
160 A(K1) = AKP*C1 - BKP*S1
B(K1) = AKP*S1 + BKP*C1
A(K2) = AJP*C2 - BJP*S2
B(K2) = AJP*S2 + BJP*C2
A(K3) = AKM*C3 - BKM*S3
B(K3) = AKM*S3 + BKM*C3
KK = K3 + KSPAN
IF (KK.LE.NT) GO TO 150
170 KK = KK - NT + JC
IF (KK.LE.MM) GO TO 130
IF (KK.LT.KSPAN) GO TO 200
KK = KK - KSPAN + INC
IF (KK.LE.JC) GO TO 120
IF (KSPAN.EQ.JC) GO TO 350
GO TO 40
180 AKP = AKM + BJM
AKM = AKM - BJM
BKP = BKM - AJM
BKM = BKM + AJM
IF (S1.NE.0.0) GO TO 160
190 A(K1) = AKP
B(K1) = BKP
A(K2) = AJP
B(K2) = BJP
A(K3) = AKM
B(K3) = BKM
KK = K3 + KSPAN
IF (KK.LE.NT) GO TO 150
GO TO 170
200 S1 = FLOAT((KK-1)/JC)*DR*RAD
C1 = COS(S1)
S1 = SIN(S1)
MM = MIN0(KSPAN,MM+KLIM)
GO TO 140
C
C TRANSFORM FOR FACTOR OF 5 (OPTIONAL CODE)
C
210 C2 = C72**2 - S72**2
S2 = 2.0*C72*S72
220 K1 = KK + KSPAN
K2 = K1 + KSPAN
K3 = K2 + KSPAN
K4 = K3 + KSPAN
AKP = A(K1) + A(K4)
AKM = A(K1) - A(K4)
BKP = B(K1) + B(K4)
BKM = B(K1) - B(K4)
AJP = A(K2) + A(K3)
AJM = A(K2) - A(K3)
BJP = B(K2) + B(K3)
BJM = B(K2) - B(K3)
AA = A(KK)
BB = B(KK)
A(KK) = AA + AKP + AJP
B(KK) = BB + BKP + BJP
AK = AKP*C72 + AJP*C2 + AA
BK = BKP*C72 + BJP*C2 + BB
AJ = AKM*S72 + AJM*S2
BJ = BKM*S72 + BJM*S2
A(K1) = AK - BJ
A(K4) = AK + BJ
B(K1) = BK + AJ
B(K4) = BK - AJ
AK = AKP*C2 + AJP*C72 + AA
BK = BKP*C2 + BJP*C72 + BB
AJ = AKM*S2 - AJM*S72
BJ = BKM*S2 - BJM*S72
A(K2) = AK - BJ
A(K3) = AK + BJ
B(K2) = BK + AJ
B(K3) = BK - AJ
KK = K4 + KSPAN
IF (KK.LT.NN) GO TO 220
KK = KK - NN
IF (KK.LE.KSPAN) GO TO 220
GO TO 290
C
C TRANSFORM FOR ODD FACTORS
C
230 K = NFAC(I)
KSPNN = KSPAN
KSPAN = KSPAN/K
IF (K.EQ.3) GO TO 100
IF (K.EQ.5) GO TO 210
IF (K.EQ.JF) GO TO 250
JF = K
S1 = RAD/(FLOAT(K)/8.0)
C1 = COS(S1)
S1 = SIN(S1)
CK(JF) = 1.0
SK(JF) = 0.0
J = 1
240 CK(J) = CK(K)*C1 + SK(K)*S1
SK(J) = CK(K)*S1 - SK(K)*C1
K = K - 1
CK(K) = CK(J)
SK(K) = -SK(J)
J = J + 1
IF (J.LT.K) GO TO 240
250 K1 = KK
K2 = KK + KSPNN
AA = A(KK)
BB = B(KK)
AK = AA
BK = BB
J = 1
K1 = K1 + KSPAN
260 K2 = K2 - KSPAN
J = J + 1
AT(J) = A(K1) + A(K2)
AK = AT(J) + AK
BT(J) = B(K1) + B(K2)
BK = BT(J) + BK
J = J + 1
AT(J) = A(K1) - A(K2)
BT(J) = B(K1) - B(K2)
K1 = K1 + KSPAN
IF (K1.LT.K2) GO TO 260
A(KK) = AK
B(KK) = BK
K1 = KK
K2 = KK + KSPNN
J = 1
270 K1 = K1 + KSPAN
K2 = K2 - KSPAN
JJ = J
AK = AA
BK = BB
AJ = 0.0
BJ = 0.0
K = 1
280 K = K + 1
AK = AT(K)*CK(JJ) + AK
BK = BT(K)*CK(JJ) + BK
K = K + 1
AJ = AT(K)*SK(JJ) + AJ
BJ = BT(K)*SK(JJ) + BJ
JJ = JJ + J
IF (JJ.GT.JF) JJ = JJ - JF
IF (K.LT.JF) GO TO 280
K = JF - J
A(K1) = AK - BJ
B(K1) = BK + AJ
A(K2) = AK + BJ
B(K2) = BK - AJ
J = J + 1
IF (J.LT.K) GO TO 270
KK = KK + KSPNN
IF (KK.LE.NN) GO TO 250
KK = KK - NN
IF (KK.LE.KSPAN) GO TO 250
C
C MULTIPLY BY ROTATION FACTOR (EXCEPT FOR FACTORS OF 2 AND 4)
C
290 IF (I.EQ.M) GO TO 350
KK = JC + 1
300 C2 = 1.0 - CD
S1 = SD
MM = MIN0(KSPAN,KLIM)
GO TO 320
310 C2 = C1 - (CD*C1+SD*S1)
S1 = S1 + (SD*C1-CD*S1)
C
C THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C ERROR. IF ROUNDED ARITHMETIC IS USED, THEY MAY
C BE DELETED.
C
C1 = 0.5/(C2**2+S1**2) + 0.5
S1 = C1*S1
C2 = C1*C2
320 C1 = C2
S2 = S1
KK = KK + KSPAN
330 AK = A(KK)
A(KK) = C2*AK - S2*B(KK)
B(KK) = S2*AK + C2*B(KK)
KK = KK + KSPNN
IF (KK.LE.NT) GO TO 330
AK = S1*S2
S2 = S1*C2 + C1*S2
C2 = C1*C2 - AK
KK = KK - NT + KSPAN
IF (KK.LE.KSPNN) GO TO 330
KK = KK - KSPNN + JC
IF (KK.LE.MM) GO TO 310
IF (KK.LT.KSPAN) GO TO 340
KK = KK - KSPAN + JC + INC
IF (KK.LE.JC+JC) GO TO 300
GO TO 40
340 S1 = FLOAT((KK-1)/JC)*DR*RAD
C2 = COS(S1)
S1 = SIN(S1)
MM = MIN0(KSPAN,MM+KLIM)
GO TO 320
C
C PERMUTE THE RESULTS TO NORMAL ORDER---DONE IN TWO STAGES
C PERMUTATION FOR SQUARE FACTORS OF N
C
350 NP(1) = KS
IF (KT.EQ.0) GO TO 440
K = KT + KT + 1
IF (M.LT.K) K = K - 1
J = 1
NP(K+1) = JC
360 NP(J+1) = NP(J)/NFAC(J)
NP(K) = NP(K+1)*NFAC(J)
J = J + 1
K = K - 1
IF (J.LT.K) GO TO 360
K3 = NP(K+1)
KSPAN = NP(2)
KK = JC + 1
K2 = KSPAN + 1
J = 1
IF (N.NE.NTOT) GO TO 400
C
C PERMUTATION FOR SINGLE-VARIATE TRANSFORM (OPTIONAL CODE)
C
370 AK = A(KK)
A(KK) = A(K2)
A(K2) = AK
BK = B(KK)
B(KK) = B(K2)
B(K2) = BK
KK = KK + INC
K2 = KSPAN + K2
IF (K2.LT.KS) GO TO 370
380 K2 = K2 - NP(J)
J = J + 1
K2 = NP(J+1) + K2
IF (K2.GT.NP(J)) GO TO 380
J = 1
390 IF (KK.LT.K2) GO TO 370
KK = KK + INC
K2 = KSPAN + K2
IF (K2.LT.KS) GO TO 390
IF (KK.LT.KS) GO TO 380
JC = K3
GO TO 440
C
C PERMUTATION FOR MULTIVARIATE TRANSFORM
C
400 K = KK + JC
410 AK = A(KK)
A(KK) = A(K2)
A(K2) = AK
BK = B(KK)
B(KK) = B(K2)
B(K2) = BK
KK = KK + INC
K2 = K2 + INC
IF (KK.LT.K) GO TO 410
KK = KK + KS - JC
K2 = K2 + KS - JC
IF (KK.LT.NT) GO TO 400
K2 = K2 - NT + KSPAN
KK = KK - NT + JC
IF (K2.LT.KS) GO TO 400
420 K2 = K2 - NP(J)
J = J + 1
K2 = NP(J+1) + K2
IF (K2.GT.NP(J)) GO TO 420
J = 1
430 IF (KK.LT.K2) GO TO 400
KK = KK + JC
K2 = KSPAN + K2
IF (K2.LT.KS) GO TO 430
IF (KK.LT.KS) GO TO 420
JC = K3
440 IF (2*KT+1.GE.M) RETURN
KSPNN = NP(KT+1)
C
C PERMUTATION FOR SQUARE-FREE FACTORS OF N
C
J = M - KT
NFAC(J+1) = 1
450 NFAC(J) = NFAC(J)*NFAC(J+1)
J = J - 1
IF (J.NE.KT) GO TO 450
KT = KT + 1
NN = NFAC(KT) - 1
JJ = 0
J = 0
GO TO 480
460 JJ = JJ - K2
K2 = KK
K = K + 1
KK = NFAC(K)
470 JJ = KK + JJ
IF (JJ.GE.K2) GO TO 460
NP(J) = JJ
480 K2 = NFAC(KT)
K = KT + 1
KK = NFAC(K)
J = J + 1
IF (J.LE.NN) GO TO 470
C
C DETERMINE THE PERMUTATION CYCLES OF LENGTH GREATER THAN 1
C
J = 0
GO TO 500
490 K = KK
KK = NP(K)
NP(K) = -KK
IF (KK.NE.J) GO TO 490
K3 = KK
500 J = J + 1
KK = NP(J)
IF (KK.LT.0) GO TO 500
IF (KK.NE.J) GO TO 490
NP(J) = -J
IF (J.NE.NN) GO TO 500
MAXF = INC*MAXF
C
C REORDER A AND B, FOLLOWING THE PERMUTATION CYCLES
C
GO TO 570
510 J = J - 1
IF (NP(J).LT.0) GO TO 510
JJ = JC
520 KSPAN = JJ
IF (JJ.GT.MAXF) KSPAN = MAXF
JJ = JJ - KSPAN
K = NP(J)
KK = JC*K + I + JJ
K1 = KK + KSPAN
K2 = 0
530 K2 = K2 + 1
AT(K2) = A(K1)
BT(K2) = B(K1)
K1 = K1 - INC
IF (K1.NE.KK) GO TO 530
540 K1 = KK + KSPAN
K2 = K1 - JC*(K+NP(K))
K = -NP(K)
550 A(K1) = A(K2)
B(K1) = B(K2)
K1 = K1 - INC
K2 = K2 - INC
IF (K1.NE.KK) GO TO 550
KK = K2
IF (K.NE.J) GO TO 540
K1 = KK + KSPAN
K2 = 0
560 K2 = K2 + 1
A(K1) = AT(K2)
B(K1) = BT(K2)
K1 = K1 - INC
IF (K1.NE.KK) GO TO 560
IF (JJ.NE.0) GO TO 520
IF (J.NE.1) GO TO 510
570 J = K3 + 1
NT = NT - KSPNN
I = NT - INC + 1
IF (NT.GE.0) GO TO 510
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: REALS
C USED WITH 'FFT' TO COMPUTE FOURIER TRANSFORM OR INVERSE FOR REAL DATA
C-----------------------------------------------------------------------
C
SUBROUTINE REALS(A, B, N, ISN)
C
C IF ISN=-1, THIS SUBROUTINE COMPLETES THE FOURIER TRANSFORM
C OF 2*N REAL DATA VALUES, WHERE THE ORIGINAL DATA VALUES ARE
C STORED ALTERNATELY IN ARRAYS A AND B, AND ARE FIRST
C TRANSFORMED BY A COMPLEX FOURIER TRANSFORM OF DIMENSION N.
C THE COSINE COEFFICIENTS ARE IN A(1),A(2),...A(N),A(N+1)
C AND THE SINE COEFFICIENTS ARE IN B(1),B(2),...B(N),B(N+1).
C NOTE THAT THE ARRAYS A AND B MUST HAVE DIMENSION N+1.
C A TYPICAL CALLING SEQUENCE IS
C CALL FFT(A,B,N,N,N,-1)
C CALL REALS(A,B,N,-1)
C
C IF ISN=1, THE INVERSE TRANSFORMATION IS DONE, THE FIRST
C STEP IN EVALUATING A REAL FOURIER SERIES.
C A TYPICAL CALLING SEQUENCE IS
C CALL REALS(A,B,N,1)
C CALL FFT(A,B,N,N,N,1)
C THE TIME DOMAIN RESULTS ALTERNATE IN ARRAYS A AND B,
C I.E. A(1),B(1),A(2),B(2),...A(N),B(N).
C
C THE DATA MAY ALTERNATIVELY BE STORED IN A SINGLE COMPLEX
C ARRAY A, THEN THE MAGNITUDE OF ISN CHANGED TO TWO TO
C GIVE THE CORRECT INDEXING INCREMENT AND A(2) USED TO
C PASS THE INITIAL ADDRESS FOR THE SEQUENCE OF IMAGINARY
C VALUES, E.G.
C CALL FFT(A,A(2),N,N,N,-2)
C CALL REALS(A,A(2),N,-2)
C IN THIS CASE, THE COSINE AND SINE COEFFICIENTS ALTERNATE IN A.
C
DIMENSION A(1), B(1)
INC = IABS(ISN)
NF = IABS(N)
IF (NF*ISN.NE.0) GO TO 10
IERR = I1MACH(4)
WRITE (IERR,9999) N, ISN
9999 FORMAT (33H ERROR - ZERO IN REALS PARAMETERS, 2I10)
RETURN
C
10 NK = NF*INC + 2
NH = NK/2
RAD = ATAN(1.0)
DR = -4.0/FLOAT(NF)
CD = 2.0*SIN(0.5*DR*RAD)**2
SD = SIN(DR*RAD)
C
C SIN,COS VALUES ARE RE-INITIALIZED EACH LIM STEPS
C
LIM = 32
MM = LIM
ML = 0
SN = 0.0
IF (ISN.GT.0) GO TO 40
CN = 1.0
A(NK-1) = A(1)
B(NK-1) = B(1)
20 DO 30 J=1,NH,INC
K = NK - J
AA = A(J) + A(K)
AB = A(J) - A(K)
BA = B(J) + B(K)
BB = B(J) - B(K)
RE = CN*BA + SN*AB
EM = SN*BA - CN*AB
B(K) = (EM-BB)*0.5
B(J) = (EM+BB)*0.5
A(K) = (AA-RE)*0.5
A(J) = (AA+RE)*0.5
ML = ML + 1
IF (ML.EQ.MM) GO TO 50
AA = CN - (CD*CN+SD*SN)
SN = (SD*CN-CD*SN) + SN
C
C THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C ERROR. IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE
C CN=AA
C
CN = 0.5/(AA**2+SN**2) + 0.5
SN = CN*SN
CN = CN*AA
30 CONTINUE
RETURN
C
40 CN = -1.0
SD = -SD
GO TO 20
C
50 MM = MM + LIM
SN = FLOAT(ML)*DR*RAD
CN = COS(SN)
IF (ISN.GT.0) CN = -CN
SN = SIN(SN)
GO TO 30
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: REALT
C USED WITH 'FFT' OR ANY OTHER COMPLEX FOURIER TRANSFORM TO COMPUTE
C TRANSFORM OR INVERSE FOR REAL DATA
C THE DATA MAY BE EITHER SINGLE-VARIATE OR MULTI-VARIATE
C-----------------------------------------------------------------------
C
SUBROUTINE REALT(A, B, NSEG, N, NSPN, ISN)
C
C IF ISN=-1, THIS SUBROUTINE COMPLETES THE FOURIER TRANSFORM
C OF 2*N REAL DATA VALUES, WHERE THE ORIGINAL DATA VALUES ARE
C STORED ALTERNATELY IN ARRAYS A AND B, AND ARE FIRST
C TRANSFORMED BY A COMPLEX FOURIER TRANSFORM OF DIMENSION N.
C THE COSINE COEFFICIENTS ARE IN A(1),A(2),...A(N),A(N+1)
C AND THE SINE COEFFICIENTS ARE IN B(1),B(2),...B(N),B(N+1).
C NOTE THAT THE ARRAYS A AND B MUST HAVE DIMENSION N+1.
C A TYPICAL CALLING SEQUENCE IS
C CALL FFT(A,B,1,N,1,-1)
C CALL REALT(A,B,1,N,1,-1)
C
C IF ISN=1, THE INVERSE TRANSFORMATION IS DONE, THE FIRST
C STEP IN EVALUATING A REAL FOURIER SERIES.
C A TYPICAL CALLING SEQUENCE IS
C CALL REALT(A,B,1,N,1,1)
C CALL FFT(A,B,1,N,1,1)
C THE TIME DOMAIN RESULTS ALTERNATE IN ARRAYS A AND B,
C I.E. A(1),B(1),A(2),B(2),...A(N),B(N).
C
C THE DATA MAY ALTERNATIVELY BE STORED IN A SINGLE COMPLEX
C ARRAY A, THEN THE MAGNITUDE OF ISN CHANGED TO TWO TO
C GIVE THE CORRECT INDEXING INCREMENT AND A(2) USED TO
C PASS THE INITIAL ADDRESS FOR THE SEQUENCE OF IMAGINARY
C VALUES, E.G.
C CALL FFT(A,A(2),1,N,1,-2)
C CALL REALT(A,A(2),1,N,1,-2)
C IN THIS CASE, THE COSINE AND SINE COEFFICIENTS ALTERNATE IN A.
C
C THIS SUBROUTINE IS SET UP TO DO THE ABOVE-DESCRIBED OPERATION ON
C ALL SUB-VECTORS WITHIN ANY DIMENSION OF A MULTI-DIMENSIONAL
C FOURIER TRANSFORM. THE PARAMETERS NSEG, N, NSPN, AND INC
C SHOULD AGREE WITH THOSE USED IN THE ASSOCIATED CALL OF 'FFT'.
C THE FOLDING FREQUENCY COSINE COEFFICIENTS ARE STORED AT THE END
C OF ARRAY A (WITH ZEROS IN CORRESPONDING LOCATIONS IN ARRAY B),
C IN A SUB-MATRIX OF DIMENSION ONE LESS THAN THE MAIN ARRAY. THE
C DELETED DIMENSION IS THAT CORRESPONDING TO THE PARAMETER N IN
C THE CALL OF REALT. THUS ARRAYS A AND B MUST HAVE DIMENSION
C NSEG*NSPN*(N+1).
C
DIMENSION A(1), B(1)
INC = IABS(ISN)
KS = IABS(NSPN)*INC
NF = IABS(N)
NS = KS*NF
NT = IABS(NS*NSEG)
IF (ISN*NT.NE.0) GO TO 10
IERR = I1MACH(4)
WRITE (IERR,9999) NSEG, N, NSPN, ISN
9999 FORMAT (33H ERROR - ZERO IN REALT PARAMETERS, 3I10, I9)
RETURN
C
10 JC = KS
K2 = IABS(KS*NSEG) - INC
KD = NS
NH = NS/2 + 1
NN = NT - INC
NT = NT + 1
KK = 1
RAD = ATAN(1.0)
DR = -4.0/FLOAT(NF)
CD = 2.0*SIN(0.5*DR*RAD)**2
SD = SIN(DR*RAD)
C
C SIN,COS VALUES ARE RE-INITIALIZED EACH LIM STEPS
C
LIM = 32
KLIM = LIM*KS
MM = MIN0(NH,KLIM)
SN = 0.0
IF (ISN.GT.0) GO TO 70
C
20 AA = A(KK)
BA = B(KK)
B(KK) = 0
A(KK) = AA + BA
A(NT) = AA - BA
B(NT) = 0
NT = NT + JC
KK = KK + NS
IF (KK.LE.NN) GO TO 20
NT = NT - K2
KK = KK - NN
IF (KK.LE.JC) GO TO 20
CN = 1.0
30 IF (NF.EQ.1) RETURN
C
40 AA = CN - (CD*CN+SD*SN)
SN = (SD*CN-CD*SN) + SN
C
C THE FOLLOWING THREE STATEMENTS COMPENSATE FOR TRUNCATION
C ERROR. IF ROUNDED ARITHMETIC IS USED, SUBSTITUTE
C CN=AA
C
CN = 0.5/(AA**2+SN**2) + 0.5
SN = CN*SN
CN = CN*AA
50 JC = JC + KS
KD = KD - KS - KS
60 K2 = KK + KD
AA = A(KK) + A(K2)
AB = A(KK) - A(K2)
BA = B(KK) + B(K2)
BB = B(KK) - B(K2)
RE = CN*BA + SN*AB
EM = SN*BA - CN*AB
B(K2) = (EM-BB)*0.5
B(KK) = (EM+BB)*0.5
A(K2) = (AA-RE)*0.5
A(KK) = (AA+RE)*0.5
KK = KK + NS
IF (KK.LE.NN) GO TO 60
KK = KK - NN
IF (KK.LE.JC) GO TO 60
IF (KK.LE.MM) GO TO 40
IF (KK.GT.NH) RETURN
SN = FLOAT(JC/KS)*DR*RAD
CN = COS(SN)
IF (ISN.GT.0) CN = -CN
SN = SIN(SN)
MM = MIN0(NH,MM+KLIM)
GO TO 50
C
70 AA = A(KK)
BA = A(NT)
A(KK) = (AA+BA)*0.5
B(KK) = (AA-BA)*0.5
NT = NT + JC
KK = KK + NS
IF (KK.LE.NN) GO TO 70
NT = NT - K2
KK = KK - NN
IF (KK.LE.JC) GO TO 70
CN = -1.0
SD = -SD
GO TO 30
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: SORTG
C SORTS ARRAY A INTO INCREASING ORDER, FROM A(1) TO A(N)
C THE ARRAY TAG IS PERMUTED THE SAME AS ARRAY A
C-----------------------------------------------------------------------
C
SUBROUTINE SORTG(A, N, TAG)
C
C TO SORT N ELEMENTS STARTING WITH A(K), CALL WITH A(K) AND TAG(K).
C AN EARLIER VERSION OF THIS ALGORITHM, WITHOUT THE TAG ARRAY, WAS
C PUBLISHED BY R.C. SINGLETON AS ACM ALGORITHM 347,
C COMM. ACM 12 (MARCH 1969), 1865-1866. THE CURRENT VERSION
C SOLVES A MACHINE-DEPENDENT PROBLEM PRESENT IN THE EARLIER
C VERSION AND ALMOST ALL OTHER SORT SUBROUTINES. ON MANY
C COMPUTERS, COMPARING A VERY LARGE NEGATIVE NUMBER WITH A
C VERY LARGE POSITIVE NUMBER GIVES A WRONG RESULT AND A BAD SORT.
C THIS PROBLEM WAS NOTED BY R. GRIFFIN AND K.A. REDISH, "REMARK
C ON ALGORITHM 347,", COMM. ACM 13 (JANUARY 1970), 54.
C THE PROBLEM IS AVOIDED HERE BY AN INITIAL SPLIT ON ZERO.
C TIME IS PROPORTIONAL TO N*LOG(N)
C AS FAR AS THE AUTHOR IS AWARE, NO FASTER IN-PLACE SORT METHOD HAS
C BEEN PUBLISHED SINCE THE ORIGINAL APPEARANCE OF THIS ALGORITHM.
C
C WORKING STORAGE ARRAYS IL AND IU SHOULD HAVE DIMENSION
C INT(ALOG(FLOAT(N))/ALOG(2.0))
C A DIMENSION OF 20 ALLOWS VALUES OF N UP TO 2**21-1
C
DIMENSION A(1), IU(20), IL(20), TAG(1)
M = 1
I = 1
J = N
K = I
L = J
IF (I.GE.J) RETURN
T = 0
IF (A(I)) 30, 30, 10
10 IF (A(L)) 90, 90, 20
20 L = L - 1
IF (L-I) 70, 70, 10
30 IF (A(J)) 40, 110, 110
40 IF (A(K)) 50, 90, 90
50 K = K + 1
IF (J-K) 70, 70, 40
60 IF (I.GE.J) GO TO 140
70 K = I
IJ = (J+I)/2
T = A(IJ)
IF (A(I).LE.T) GO TO 80
A(IJ) = A(I)
A(I) = T
T = A(IJ)
TG = TAG(IJ)
TAG(IJ) = TAG(I)
TAG(I) = TG
80 L = J
IF (A(J).GE.T) GO TO 110
A(IJ) = A(J)
A(J) = T
T = A(IJ)
TG = TAG(IJ)
TAG(IJ) = TAG(J)
TAG(J) = TG
IF (A(I).LE.T) GO TO 110
A(IJ) = A(I)
A(I) = T
T = A(IJ)
TG = TAG(IJ)
TAG(IJ) = TAG(I)
TAG(I) = TG
GO TO 110
90 TT = A(L)
100 A(L) = A(K)
A(K) = TT
TG = TAG(L)
TAG(L) = TAG(K)
TAG(K) = TG
110 L = L - 1
IF (A(L).GT.T) GO TO 110
TT = A(L)
120 K = K + 1
IF (A(K).LT.T) GO TO 120
IF (K.LE.L) GO TO 100
IF (L-I.LE.J-K) GO TO 130
IL(M) = I
IU(M) = L
I = K
M = M + 1
GO TO 150
130 IL(M) = K
IU(M) = J
J = L
M = M + 1
GO TO 150
140 M = M - 1
IF (M.EQ.0) RETURN
I = IL(M)
J = IU(M)
150 IF (J-I.GT.10) GO TO 70
IF (I.EQ.1) GO TO 60
I = I - 1
160 I = I + 1
IF (I.EQ.J) GO TO 140
T = A(I+1)
IF (A(I).LE.T) GO TO 160
TG = TAG(I+1)
K = I
170 A(K+1) = A(K)
TAG(K+1) = TAG(K)
K = K - 1
IF (T.LT.A(K)) GO TO 170
A(K+1) = T
TAG(K+1) = TG
GO TO 160
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: NORMAL
C GENERATES AN INDEPENDENT PAIR OF RANDOM NORMAL DEVIATES
C METHOD DUE TO G. MARSAGLIA AND T.A. BRAY,
C SIAM REVIEW, VOL. 6, NO. 3, JULY 1964. 260-264
C-----------------------------------------------------------------------
C
SUBROUTINE NORMAL(X, Y)
C
C OUTPUT: X,Y = INDEPENDENT PAIR OF RANDOM NORMAL DEVIATES
C FUNCTION UNI GENERATES PSEUDO-RANDOM NUMBER BETWEEN 0.0 AND 1.0
C
10 RX = UNI(NRM)*2.0 - 1.0
RY = UNI(NRM)*2.0 - 1.0
R = RX**2 + RY**2
IF (R.GE.1.0) GO TO 10
R = SQRT(-2.0*ALOG(R)/R)
X = RX*R
Y = RY*R
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: RMS
C COMPUTES RMS ERROR FOR TRANSFORM-INVERSE PAIR
C-----------------------------------------------------------------------
C
SUBROUTINE RMS(A, B, C, D, N, EA, EB)
C
C ARRAYS: A,B = TRANSFORM, INVERSE RESULTS
C C,D = ORIGINAL DATA
C INPUT: N = DIMENSION OF ARRAYS A, B, C AND D
C OUTPUT: EA,EB = RMS ERRORS FOR A AND B ARRAYS
C
DIMENSION A(1), B(1), C(1), D(1)
SSA = 0.0
SSB = 0.0
DO 10 J=1,N
SSA = (A(J)-C(J))**2 + SSA
SSB = (B(J)-D(J))**2 + SSB
10 CONTINUE
EA = SQRT(SSA/FLOAT(N))
EB = SQRT(SSB/FLOAT(N))
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: ISTKGT(NITEMS,ITYPE)
C ALLOCATES WORKING STORAGE FOR NITEMS OF ITYPE, AS FOLLOWS
C
C 1 - LOGICAL
C 2 - INTEGER
C 3 - REAL
C 4 - DOUBLE PRECISION
C 5 - COMPLEX
C
C-----------------------------------------------------------------------
C
INTEGER FUNCTION ISTKGT(NITEMS, ITYPE)
C
COMMON /CSTAK/ DSTAK(2500)
C
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
INTEGER ISIZE(5)
C
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (ISTAK(1),LOUT)
EQUIVALENCE (ISTAK(2),LNOW)
EQUIVALENCE (ISTAK(3),LUSED)
EQUIVALENCE (ISTAK(4),LMAX)
EQUIVALENCE (ISTAK(5),LBOOK)
EQUIVALENCE (ISTAK(6),ISIZE(1))
C
ISTKGT = (LNOW*ISIZE(2)-1)/ISIZE(ITYPE) + 2
I = ((ISTKGT-1+NITEMS)*ISIZE(ITYPE)-1)/ISIZE(2) + 3
IF (I.GT.LMAX) GO TO 10
ISTAK(I-1) = ITYPE
ISTAK(I) = LNOW
LOUT = LOUT + 1
LNOW = I
LUSED = MAX0(LUSED,LNOW)
RETURN
C
10 IERR = I1MACH(4)
WRITE (IERR,9999) I
9999 FORMAT (1H , 39HOVERFLOW OF COMMON ARRAY ISTAK --- NEED, I10)
WRITE (IERR,9998) (ISTAK(J),J=1,10), ISTAK(LNOW-1), ISTAK(LNOW)
9998 FORMAT (1X, I5, 11I6)
STOP
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: ISTKRL(K)
C DE-ALLOCATES THE LAST K WORKING STORAGE AREAS
C-----------------------------------------------------------------------
C
SUBROUTINE ISTKRL(K)
C
COMMON /CSTAK/ DSTAK(2500)
C
DOUBLE PRECISION DSTAK
INTEGER ISTAK(5000)
C
EQUIVALENCE (DSTAK(1),ISTAK(1))
EQUIVALENCE (ISTAK(1),LOUT)
EQUIVALENCE (ISTAK(2),LNOW)
EQUIVALENCE (ISTAK(3),LUSED)
EQUIVALENCE (ISTAK(4),LMAX)
EQUIVALENCE (ISTAK(5),LBOOK)
C
IN = K
C
IF (LBOOK.LE.LNOW .AND. LNOW.LE.LUSED .AND. LUSED.LE.LMAX) GO TO
* 10
IERR = I1MACH(4)
WRITE (IERR,9999)
9999 FORMAT (53H WARNING...ISTAK(2),ISTAK(3),ISTAK(4) OR ISTAK(5) HIT)
WRITE (IERR,9997) (ISTAK(J),J=1,10), ISTAK(LNOW-1), ISTAK(LNOW)
C
10 IF (IN.LE.0) RETURN
IF (LBOOK.GT.ISTAK(LNOW) .OR. ISTAK(LNOW).GE.LNOW-1) GO TO 20
LOUT = LOUT - 1
LNOW = ISTAK(LNOW)
IN = IN - 1
GO TO 10
C
20 IERR = I1MACH(4)
WRITE (IERR,9998)
9998 FORMAT (45H WARNING...POINTER AT ISTAK(LNOW) OVERWRITTEN/11X,
* 27HDE-ALLOCATION NOT COMPLETED)
WRITE (IERR,9997) (ISTAK(J),J=1,10), ISTAK(LNOW-1), ISTAK(LNOW)
9997 FORMAT (1X, I5, 11I6)
RETURN
C
END
Return to Main Software Page