Program for Minimum-p Synthesis of Recursive Digital Filters
Return to Main Software Page
C
C-----------------------------------------------------------------------
C MAIN PROGRAM: COMBINED AMPLITUDE AND GROUP DELAY SYNTHESIS USING A
C MINIMUM P ERROR CRITERION AND THE FLETCHER POWELL
C SUBROUTINE FOR FUNCTION MINIMIZATION
C
C AUTHOR: A. G. DECZKY
C INSTITUT DE MICROTECHNIQUE, UNIVERSITE DE NEUCHATEL
C RUE DE LA MALADIERE 71, 2000 NEUCHATEL, SWITZERLAND
C
C REFERENCE: 'SYNTHESIS OF RECURSIVE DIGITAL FILTERS USING THE
C MINIMUM P ERROR CRITERION', A. G. DECZKY, IEEE TRANS.
C AUDIO ELECTROACOUST., VOL. AU-20, NO. 5, PP. 257-263,
C OCTOBER 1972.
C-----------------------------------------------------------------------
C
LOGICAL DEB, DEBFP, LONG
INTEGER WEIGHT
COMMON /CM05/ AS(201), DS(201), WF(201), WG(201)
COMMON /CM10/ NP, NPD
COMMON /CM15/ DEB
COMMON /CM20/ C0(20), C1(20), C2(20), D1(20), D2(20)
COMMON /CM25/ DEBFP, LONG
COMMON /CM30/ FV(100), XIT(100), KOUNT
COMMON /CM35/ FD(201), FDN(201)
COMMON /CM40/ FIK(10), NPK(10), NK
COMMON /CM45/ WEIGHT(10)
COMMON /CM50/ N1, N2, N3, N4, N5, MODE
COMMON /CM55/ ALFA, BETA, NIT
COMMON /CM60/ WFK(10), WGK(10)
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
DIMENSION X(40), G(40), IP(10), IQ(10)
C
C INITIALIZE MACHINE DEPENDENT CONSTANTS
C
IND = I1MACH(1)
INP = I1MACH(1)
IOUTD = I1MACH(2)
IOUTP = I1MACH(3)
PI = 4.0*ATAN(1.0)
R1 = R1MACH(1)
SPM = SQRT(R1)
C
C INPUT DATA
C
C MODE = TYPE OF APPROXIMATION
C = 1 MAGNITUDE ONLY
C = 2 GROUP DELAY ONLY
C = 3 GROUP DELAY EQUALISATION
C = 4 MAGNITUDE AND GROUP DELAY
C
READ (IND,9999) MODE
C
C NZR = NO OF REAL ZEROS
C NZU = NO OF ZEROS ON THE UNIT CIRCLE
C NZC = NO OF OF COMPLEX ZEROS
C NPR = NO OF REAL POLES
C NPC = NO OF OF COMPLEX POLES
C
READ (IND,9999) NZR, NZU, NZC, NPR, NPC
N1 = NZR
N2 = N1 + NZU
N3 = N2 + NZC*2
N4 = N3 + NPR
N5 = N4 + NPC*2 - 1
N = N5 + 2
IF (MODE.EQ.4) N = N + 1
C
C X = INITIAL PARAMETER VECTOR
C ZEROES AND POLES SPECIFIED IN POLAR COORDINATES AS FOLLOWS
C REAL ZEROS COME FIRST, THEN ZEROS ON THE UNIT CIRCLE, THEN COMPLEX
C ZEROS, THEN REAL POLES AND FINALLY COMPLEX POLES
C REAL ZEROES (POLES) ARE SPECIFIED AS RADIUS1,RADIUS2,...
C ZEROES ON THE UNIT CIRCLE ARE SPECIFIED AS ANGLE1,ANGLE2,...
C COMPLEX ZEROES (POLES) ARE SPECIFIED AS ANGLE1,RADIUS1,ANGLE2,
C RADIUS2,...
C WHERE THE ANGLES ARE SPECIFIED IN HZ AS PHI/(2 PI)*FSAMPLE
C AND ONLY ZEROES (POLES) WITH POSITIVE ANGLES ARE SPECIFIED
C THEIR COMPLEX CONJUGATE PAIRS BEING UNDERSTOOD
C THE LAST PARAMETER(S) IS (ARE): K0 IF MODE=1, TAU0 IF MODE=2,3,
C K0,TAU0 IF MODE=4
C
READ (IND,9998) (X(I),I=1,N)
C
C NK = NO OF INTERVALS
C
READ (IND,9999) NK
C
C NPK(I) = NO. OF POINTS ON THE INTERVAL FIK(I) - FIK(I+1)
C
READ (IND,9999) (NPK(I),I=1,NK)
C
C WEIGHT(I) = SPACING OF THE POINTS NPK(I) ON THE I TH INTERVAL
C = 1 POINTS EQUISPACED
C = 2 POINTS SPACED ON A SINE ABSCISSA, HALF CYCLE
C = 3 POINTS SPACED ON A COS ABSCISSA, HALF CYCLE
C = 4 POINTS SPACED ON A COS ABSCISSA, FULL CYCLE
C
READ (IND,9999) (WEIGHT(I),I=1,NK)
C
C WFK(I) = RELATIVE WEIGHTING OF THE I TH INTERVAL FOR THE MAGNITUDE
C
READ (IND,9998) (WFK(I),I=1,NK)
C WGK(I) = RELATIVE WEIGHTING OF THE I TH INTERVAL FOR THE GROUP DELAY
C
READ (IND,9998) (WGK(I),I=1,NK)
C
NK = NK + 1
C
C FIK(I) = STARTING POINT OF THE I TH APPROX. INTERVAL
C
READ (IND,9998) (FIK(I),I=1,NK)
C
C IPK = NO OF SUCCESSIVE INDICES OF APPROXIMATION
C
READ (IND,9999) IPK
C
C IP(II) = SUCCESSIVE INDICES OF APPROXIMATION FOR THE MAGNITUDE
C
READ (IND,9999) (IP(II),II=1,IPK)
C
C IQ(II) = SUCCESSIVE INDICES OF APPROXIMATION FOR THE GROUP DELAY
C
READ (IND,9999) (IQ(II),II=1,IPK)
C
C INV = NUMBER OF TOTAL CYCLES IF NO CONVERGENCE
C
READ (IND,9999) INV
C
C F1 = LOWER FREQUENCY LIMIT FOR THE FREQUENCY RESPONSE
C F2 = UPPER FREQUENCY LIMIT FOR THE FREQUENCY RESPONSE
C FSAMPL = SAMPLING FREQUENCY
C IZ = NUMBER OF POINTS FOR THE FREQUENCY RESPONSE
C (IF IZ LE 0 NO RESPONSE CALCULATED )
C
READ (IND,9998) F1, F2, FSAMPL
READ (IND,9999) IZ
WS = FSAMPL/2.0
C
C ALFA= RELATIVE WEIGHTING OF THE MAGNITUDE
C RELATIVE WEIGHTING OF THE GROUP DELAY = 1.0 - ALFA
C THIS PARAMETER IS ONLY USED IF MODE = 4
C
READ (IND,9998) ALFA
BETA = 1.0 - ALFA
C
C LONG = LOGICAL PARAMETER FOR TYPE OF OUTPUT
C = .TRUE. FOR LONG OUTPUT
C = .FALSE. FOR SHORT OUTPUT
C
C
READ (IND,9997) LONG
C
C LIMIT = NUMBER OF ITERATIONS PER CYCLE BEFORE RESTARTING WITH
C STEEPEST DESCENT
C
LIMIT = 4*N
C
C EPS = ERROR PARAMETER FOR FMFP SUBROUTINE
C
EPS = 1.E-3
C
C DEBUG = .TRUE . FOR DEBUGGING PRINTOUT
C
DEB = .FALSE.
C
C DEBFP = .TRUE . FOR DEBUGGING PRINTOUT FROM FMFP SUBROUTINE
C
DEBFP = .FALSE.
READ (IND,9997) DEB, DEBFP
C
9999 FORMAT (I4)
9998 FORMAT (F16.9)
9997 FORMAT (L1)
CALL COSYFP(X, G, IZ, F1, F2, WS, N, EPS, IP, IQ, LIMIT, INV, IPK)
STOP
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: COSYFP
C MAIN SUBROUTINE FOR AMPLITUDE AND/OR G.D. SYNTHESIS USING FP
C-----------------------------------------------------------------------
C
SUBROUTINE COSYFP(X, G, IZ, F1, F2, WS, N, EPS, IP, IQ, LIMIT,
* INV, IPK)
LOGICAL DEB, DEBFP, LONG
INTEGER WEIGHT
COMMON /CM05/ AS(201), DS(201), WF(201), WG(201)
COMMON /CM10/ NP, NPD
COMMON /CM15/ DEB
COMMON /CM20/ C0(20), C1(20), C2(20), D1(20), D2(20)
COMMON /CM25/ DEBFP, LONG
COMMON /CM30/ FV(100), XIT(100), KOUNT
COMMON /CM35/ FD(201), FDN(201)
COMMON /CM40/ FIK(10), NPK(10), NK
COMMON /CM45/ WEIGHT(10)
COMMON /CM50/ N1, N2, N3, N4, N5, MODE
COMMON /CM55/ ALFA, BETA, NIT
COMMON /CM60/ WFK(10), WGK(10)
COMMON /CM65/ IPX, IQX
COMMON /CM70/ A(201), B(201)
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
DIMENSION GDH(201)
DIMENSION X(1), G(1), H(201), U0(20), V0(20), U(20), V(20),
* A0(20), A1(20), A2(20), B0(20), B1(20), B2(20), AMPL(201),
* PHI(201), DB(201), GD(201), IP(10), IQ(10)
C
C INITIALIZE CONSTANTS
C
NPD = 0
NPA = 0
NP = 0
KP = 1
KOUNT = 0
NIT = 0
IER = 0
NIT = 0
NKOUNT = 0
IQX = IQ(1)
IPX = IP(1)
NN = N/4
CONST = 1.0
DO 10 I=1,20
A0(I) = 1.0
B0(I) = 1.0
A1(I) = 0.0
A2(I) = 0.0
B1(I) = 0.0
B2(I) = 0.0
10 CONTINUE
C
C INITIALIZE X(N) AND X(N-1) (IF NECESSARY)
C
GO TO (20, 30, 30, 40), MODE
20 X(N) = 1.0
GO TO 50
30 X(N) = 0.0
GO TO 50
40 X(N-1) = 1.0
X(N) = 0.0
50 WSPI = PI/WS
C
C COMPUTE FREQUENCY POINTS ON THE VARIOUS INTERVALS
C
DO 110 K=2,NK
NPX = NPK(K-1)
NPX1 = NPX - 1
FNPX = FLOAT(NPX)
FI2 = FIK(K)
FI1 = FIK(K-1)
IWG = WEIGHT(K-1)
WGK1 = WGK(K-1)
WFK1 = WFK(K-1)
C
IF (NPX.EQ.0) GO TO 110
C
DO 100 J=1,NPX
JK = J + NP
WG(JK) = WGK1
WF(JK) = WFK1
C
FLJNP = FLOAT(J-1)/FNPX
GO TO (60, 70, 80, 90), IWG
60 FDN(JK) = FLJNP*(FI2-FI1) + FI1
GO TO 100
70 FDN(JK) = SIN(FLJNP*PI/2.0)*(FI2-FI1) + FI1
GO TO 100
80 FDN(JK) = COS(FLJNP*PI/2.0)*(FI1-FI2) + FI2
GO TO 100
90 FDN(JK) = 0.5*(FI1+FI2) + COS(FLJNP*PI)*(FI1-FI2)*0.5
100 CONTINUE
C
IF (WFK1.NE.0.) NPA = NPA + NPX
IF (WGK1.NE.0.) NPD = NPD + NPX
NP = NP + NPX
110 CONTINUE
C
DO 120 J=1,NP
FD(J) = FDN(J)*WSPI
120 CONTINUE
C
N1P = N1 + 1
N2P = N2 + 1
N4P = N4 + 1
C
C NORMALISE THE ANGLE OF THE POLES AND ZEROES
C
IF (N1P.GT.N2) GO TO 140
DO 130 J=N1P,N2
X(J) = X(J)/WS
130 CONTINUE
140 IF (N2P.GT.N3) GO TO 160
DO 150 J=N2P,N3,2
X(J) = X(J)/WS
150 CONTINUE
160 IF (N4P.GT.N5) GO TO 180
DO 170 J=N4P,N5,2
X(J) = X(J)/WS
170 CONTINUE
C
C READ IN FILTER COEFFICIENTS
C
180 IF (MODE.NE.3) GO TO 190
READ (INP,9999) CONST0, M
9999 FORMAT (E14.7, I4)
CALL INCOCA(M, CONST0, C0, C1, C2, D1, D2)
CALL DGSPEC(GDH, M)
C
C SET UP ARRAYS OF DESIRED FUNCTIONS FOR MAGN. AND GD
C
190 DO 240 J=1,NP
FJ = FDN(J)
GO TO (200, 210, 220, 230), MODE
200 AS(J) = FS(FJ)
GO TO 240
210 DS(J) = DGS(FJ)
GO TO 240
220 DS(J) = -GDH(J)
GO TO 240
230 AS(J) = FS(FJ)
DS(J) = DGS(FJ)
240 CONTINUE
C
CALL FUNCT(N, X, F, G)
C
C COMPUTE INITIAL APPROX. TO X(N) AND X(N-1)
C
CALL XNXM1(X, N, NP)
C
GO TO (250, 260, 260, 270), MODE
250 IF (LONG) WRITE (IOUTD,9983) X(N)
GO TO 280
260 IF (LONG) WRITE (IOUTD,9982) X(N)
GO TO 280
270 IF (LONG) WRITE (IOUTD,9981) X(N-1), X(N)
280 CALL FUNCT(N, X, F, G)
IF (LONG) WRITE (IOUTD,9998) F, IER, KOUNT, NIT
IF (LONG) WRITE (IOUTD,9994)
IF (LONG) WRITE (IOUTD,9995) (J,X(J),J,G(J),J=1,N)
C
C MAIN LOOP CALLING THE F.P. SUBROUTINE
C
290 CALL FMFP(N, X, F, G, EPS, LIMIT, IER, H)
C
NKOUNT = NKOUNT + KOUNT
IF (LONG) WRITE (IOUTD,9998) F, IER, KOUNT, NIT
IF (LONG) WRITE (IOUTD,9995) (J,X(J),J,G(J),J=1,N)
C
IF (IER.EQ.0) GO TO 310
IF (INV-1) 310, 310, 300
300 INV = INV - 1
GO TO 290
310 IF (IER.NE.0 .OR. KP.GE.IPK) GO TO 320
KP = KP + 1
IQX = IQ(KP)
IPX = IP(KP)
GO TO 290
C
C COMPUTE QUADRATIC FACTORS, POLES AND ZEROS FROM ARG. VECTOR
C
320 CALL FACTOR(X, A1, A2, B1, B2, U0, V0, U, V, CONST, JA, JB, JZ,
* JP, WS)
C
NN = MAX0(JA-1,JB-1)
JZ1 = JZ - 1
JP1 = JP - 1
C
IF (MODE.NE.3) GO TO 350
JA = JA + M
JB = JB + M
NNP = NN + 1
NM = NN + M
C
IF (NNP.GT.NM) GO TO 340
DO 330 J=NNP,NM
JM = J - NN
A1(J) = C1(JM)
A2(J) = C2(JM)
B1(J) = D1(JM)
B2(J) = D2(JM)
330 CONTINUE
340 NN = NM
C
350 IF (DEB) WRITE (IOUTD,9980) JA, JB
IF (DEB) WRITE (IOUTD,9979) (A1(J),A2(J),B1(J),B2(J),J=1,NN)
C
GO TO (360, 390, 370, 380), MODE
360 CONST = X(N)
GO TO 390
370 CONST = CONST*CONST0
GO TO 390
380 CONST = X(N-1)
C
C OUTPUT SECTION
C
390 WRITE (IOUTD,9996)
WRITE (IOUTD,9986) IPX, IQX
WRITE (IOUTD,9998) F, IER, NKOUNT, NIT
WRITE (IOUTD,9993)
WRITE (IOUTD,9995) (J,X(J),J,G(J),J=1,N)
IF (JZ1.GT.0) WRITE (IOUTD,9989) (U0(J),V0(J),J=1,JZ1)
IF (JP1.GT.0) WRITE (IOUTD,9988) (U(J),V(J),J=1,JP1)
WRITE (IOUTD,9987) CONST
WRITE (IOUTD,9985) (A0(J),A1(J),A2(J),J=1,NN)
WRITE (IOUTD,9984) (B0(J),B1(J),B2(J),J=1,NN)
IF (IZ.LE.0) GO TO 400
CALL FREDIC(NN, CONST, A0, A1, A2, B1, B2, F1, F2, WS, IZ, AMPL,
* PHI, DB, GD, FD)
WRITE (IOUTD,9997) (FD(J),AMPL(J),DB(J),PHI(J),GD(J),J=1,IZ)
C
400 WRITE (IOUTP,9992) NN
WRITE (IOUTP,9991) CONST, NN
WRITE (IOUTP,9990) (A0(J),A1(J),A2(J),B1(J),B2(J),J=1,NN)
RETURN
C
C FORMATS
C
9998 FORMAT (/3H F=, E12.5, 5X, 4HIER=, I2, 5X, 17HNO OF ITERATIONS=,
* I3//1X, 24HNO OF FUNCT EVALUATIONS=, I4)
9997 FORMAT (1H1//19H FREQUENCY RESPONSE//1X, 9HFREQUENCY, 6X,
* 9HMAGNITUDE, 5X, 10HLOSS IN DB, 8X, 5HPHASE, 7X, 9HGROUP DEL,
* 2HAY//(E11.4, 4E15.5))
9996 FORMAT (1H1//49H SYNTHESIS OF DIGITAL FILTER WITH SPEC. MAGNITUDE,
* 20H AND/OR GROUP DELAY /)
9995 FORMAT (/(10X, 2HX(, I2, 2H)=, E15.8, 10X, 2HG(, I2, 2H)=, E15.8))
9994 FORMAT (//33H INITIAL ARG. VECTOR AND GRADIENT)
9993 FORMAT (//34H COMPUTED ARG. VECTOR AND GRADIENT)
9992 FORMAT (52H COEFFICIENTS OF DIGITAL FILTER WITH SPEC. AMPLITUDE,
* 19H AND/OR GROUP DELAY, / 3H N=, I2)
9991 FORMAT (1X, E14.7, I4)
9990 FORMAT (1X, 5E14.7)
9989 FORMAT (/24X, 13HZ PLANE ZEROS//12X, 12HANGLE/2PI*FS, 16X,
* 6HRADIUS//(2E25.8))
9988 FORMAT (/24X, 13HZ PLANE POLES//12X, 12HANGLE/2PI*FS, 16X,
* 6HRADIUS//(2E25.8))
9987 FORMAT (/1X, 6HCONST=, E15.8)
9986 FORMAT (/36H THE ERROR CRITERION USED IS MINIMUM, I2, 9H FOR THE ,
* 13HMAGNITUDE AND//29X, 7HMINIMUM, I2, 14H FOR THE GROUP,
* 6H DELAY)
9985 FORMAT (//24X, 22HNUMERATOR COEFFICIENTS///14X, 2HA0, 20X, 2HA1,
* 20X, 2HA2//(3E22.8))
9984 FORMAT (//23X, 24HDENOMINATOR COEFFICIENTS///14X, 2HB0, 20X, 2HB1,
* 20X, 2HB2//(3E22.8))
9983 FORMAT (1H1//5H K0 =, E15.8)
9982 FORMAT (1H1//5H T0 =, E15.8)
9981 FORMAT (1H1//5H K0 =, E15.8, 10X, 5H T0 =, E15.8)
9980 FORMAT (1H1//2I10//)
9979 FORMAT (//(4E18.5))
C
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: XNXM1
C SUBROUTINE TO CALCULATE INITIAL APPROX. TO X(N),X(N-1)
C-----------------------------------------------------------------------
C
SUBROUTINE XNXM1(X, N, NP)
DIMENSION X(N)
COMMON /CM05/ AS(201), DS(201), WF(201), WG(201)
COMMON /CM50/ N1, N2, N3, N4, N5, MODE
COMMON /CM70/ A(201), B(201)
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
C
AGD = 0.0
AGN = 0.0
AD = 0.0
AN = 0.0
DO 40 J=1,NP
GO TO (10, 20, 20, 30), MODE
10 WFJ = WF(J)
IF (WFJ.LT.SPM) WFJ = SPM
HJ = A(J)/WFJ + AS(J)
WHJ = HJ*WFJ
AN = AN + WHJ*AS(J)
AD = AD + HJ*WHJ
GO TO 40
20 AGN = AGN + B(J)
AGD = AGD + WG(J)
GO TO 40
30 WFJ = WF(J)
IF (WFJ.LT.SPM) WFJ = SPM
HJ = A(J)/WFJ + AS(J)
WHJ = HJ*WFJ
AN = AN + WHJ*AS(J)
AD = AD + HJ*WHJ
AGN = AGN + B(J)
AGD = AGD + WG(J)
40 CONTINUE
C
GO TO (50, 60, 60, 70), MODE
50 X(N) = AN/AD
GO TO 80
60 X(N) = AGN/AGD
GO TO 80
70 X(N) = AGN/AGD
X(N-1) = AN/AD
80 RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: RANGE
C FUNCTION TO LIMIT THE RANGE OF ANGLES (X) TO 0 TO PI
C-----------------------------------------------------------------------
C
FUNCTION RANGE(X)
C
X = ABS(X)
IX = IFIX(X)
MX = 2*(IX/2)
RANGE = X - FLOAT(MX)
IF (RANGE.GT.1.0) RANGE = 2.0 - RANGE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: DGSPEC
C TO CALCULATE THE GROUP DELAY OF THE NON-EQUALISED FILTER
C-----------------------------------------------------------------------
C
SUBROUTINE DGSPEC(GDH, M)
DIMENSION GDH(201)
COMMON /CM10/ NP, NPD
COMMON /CM20/ C0(20), C1(20), C2(20), D1(20), D2(20)
COMMON /CM35/ FD(201), FDN(201)
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
COMPLEX NUMC, NUMD, DENC, DEND
C
DO 20 J=1,NP
FJ = FD(J)
GDH(J) = 0.
CS = COS(FJ)
SN = SIN(FJ)
CS2 = CS + CS
SN2 = SN + SN
C
DO 10 I=1,M
C11 = C0(I) + C2(I)
C12 = C1(I)
C22 = C0(I) - C2(I)
D11 = 1. + D2(I)
D12 = D1(I)
D22 = 1. - D2(I)
NUMC = CMPLX(CS2+C12,SN2)
NUMD = CMPLX(CS2+D12,SN2)
DENC = CMPLX(C11*CS+C12,C22*SN)
DEND = CMPLX(D11*CS+D12,D22*SN)
GDH(J) = GDH(J) + REAL(NUMD/DEND) - REAL(NUMC/DENC)
10 CONTINUE
C
20 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: INCOCA
C READ COEFFICIENTS FOR ONE DIGITAL CASCADE FILTER
C E14.7 FORMAT, FILE ON UNIT INP
C-----------------------------------------------------------------------
C
SUBROUTINE INCOCA(NTERM, CONST, A0, A1, A2, B1, B2)
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
REAL A0(NTERM), A1(NTERM), A2(NTERM), B1(NTERM), B2(NTERM)
DO 10 N=1,NTERM
READ (INP,9999) A0(N), A1(N), A2(N), B1(N), B2(N)
10 CONTINUE
9999 FORMAT (5E14.7)
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FREDIC
C SUBROUTINE TO COMPUTE THE FREQUENCY RESPONSE OF A DIGITAL FILTER IN
C CASCADE FORM
C-----------------------------------------------------------------------
C
SUBROUTINE FREDIC(NTERM, CONST, A0, A1, A2, B1, B2, FU, FO, WS,
* IZ, AMPL, PHI, DB, GD, FD)
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
COMPLEX N, D, GWS, NN, DD, N1, N2
DIMENSION A0(NTERM), A1(NTERM), A2(NTERM), B1(NTERM), B2(NTERM)
DIMENSION AMPL(IZ), PHI(IZ), DB(IZ), FD(IZ), GD(IZ)
DIMENSION AA(40), AB(40), BA(40), BB(40)
C
WSXPI = PI/WS
WU = FU
Z = (FO-FU)/FLOAT(IZ-1)
C
DO 10 I=1,NTERM
AA(I) = A0(I) + A2(I)
AB(I) = A0(I) - A2(I)
BA(I) = 1.0 + B2(I)
BB(I) = 1.0 - B2(I)
10 CONTINUE
C
DO 40 J=1,IZ
WT = WU + Z*FLOAT(J-1)
FD(J) = WT
WT = WT*WSXPI
CWT = COS(WT)
SWT = SIN(WT)
D = CMPLX(1.0,0.)
N = D
QL = 0.0
QQ = QL
C
DO 20 I=1,NTERM
NN = CMPLX(AA(I)*CWT+A1(I),AB(I)*SWT)
DD = CMPLX(BA(I)*CWT+B1(I),BB(I)*SWT)
N2 = CMPLX(2.0*CWT+B1(I),2.0*SWT)
N1 = CMPLX(2.0*A0(I)*CWT+A1(I),2.0*A0(I)*SWT)
N = N*NN
D = D*DD
C
IF (CABS(DD).LT.SPM) DD = CMPLX(SPM,SPM)
Q2 = REAL(N2/DD)
C
IF (CABS(NN).LT.SPM) NN = CMPLX(SPM,SPM)
Q1 = REAL(N1/NN)
IF (ABS(Q1-QL).GT.50.) Q1 = 1.0
QL = Q1
QQ = QQ + Q2 - Q1
20 CONTINUE
C
GD(J) = QQ
C
IF (CABS(D).LT.SPM) D = CMPLX(SPM,SPM)
GWS = N/D
C
AMPL(J) = CONST*CABS(GWS)
GABS = AMPL(J)
IF (GABS.NE.0.) GO TO 30
C
PHI(J) = 0.0
DB(J) = 1000.
GO TO 40
C
30 GWSI = AIMAG(GWS)
GWSR = REAL(GWS)
PHI(J) = ATAN(GWSI/GWSR)
DB(J) = -ALOG10(GABS)*20.0
40 CONTINUE
C
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FACTOR
C SUBROUTINE TO RECALCULATE QUADRATIC FACTORS
C AND POLES AND ZEROS FROM PARAMETER VECTOR X
C-----------------------------------------------------------------------
C
SUBROUTINE FACTOR(X, A1, A2, B1, B2, U0, V0, U, V, CONST, JA, JB,
* JZ, JP, WS)
COMMON /CM50/ N1, N2, N3, N4, N5, MODE
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
DIMENSION X(1), A1(1), A2(1), B1(1), B2(1), U0(1), V0(1)
DIMENSION U(1), V(1)
C
J = 0
JA = 1
JB = 1
JZ = 1
CONST = 1.0
10 IF (J.GE.N1) GO TO 40
J = J + 1
R = X(J)
GO TO (20, 30, 30, 30), MODE
20 A2(JA) = 0.0
A1(JA) = -R
U0(JZ) = 0.
V0(JZ) = R
JA = JA + 1
JZ = JZ + 1
GO TO 80
30 A1(JA) = -R - 1.0/R
A2(JA) = 1.0
U0(JZ) = 0.0
V0(JZ) = R
U0(JZ+1) = 0.0
V0(JZ+1) = 1.0/R
JA = JA + 1
JZ = JZ + 2
GO TO 80
C
40 IF (J.GE.N2) GO TO 50
J = J + 1
XJ = X(J)
XJ = RANGE(XJ)
A1(JA) = -2.0*COS(XJ*PI)
A2(JA) = 1.0
V0(JZ) = 1.0
U0(JZ) = XJ*WS
JA = JA + 1
JZ = JZ + 1
GO TO 80
C
50 IF (J.GE.N3) GO TO 90
J = J + 2
XJM = X(J-1)
XJM = RANGE(XJM)
CSP = COS(XJM*PI)
R = X(J)
GO TO (60, 70, 70, 70), MODE
60 A2(JA) = R*R
A1(JA) = -(R+R)*CSP
U0(JZ) = XJM*WS
V0(JZ) = R
JA = JA + 1
JZ = JZ + 1
GO TO 80
70 A2(JA) = R*R
A1(JA) = -(R+R)*CSP
A2(JA+1) = 1.0/A2(JA)
A1(JA+1) = -2.0/R*CSP
V0(JZ) = R
V0(JZ+1) = 1.0/R
U0(JZ) = XJM*WS
U0(JZ+1) = XJM*WS
JA = JA + 2
JZ = JZ + 2
80 IF (J.LT.N3) GO TO 10
C
90 NN = MAX0(N4,N5)
JP = 1
100 IF (J.GE.N4) GO TO 140
J = J + 1
R = X(J)
GO TO (110, 120, 120, 110), MODE
110 B2(JB) = 0.0
B1(JB) = -R
U(JP) = 0.0
V(JP) = R
GO TO 130
120 B2(JB) = 0.0
B1(JB) = -R
A1(JA) = -1.0/R
A2(JA) = 0.0
CONST = CONST*R
U(JP) = 0.0
V(JP) = R
U0(JZ) = 0.0
V0(JZ) = R
JA = JA + 1
JZ = JZ + 1
130 JB = JB + 1
JP = JP + 1
GO TO 180
C
140 IF (J.GT.N5) GO TO 190
J = J + 2
R = X(J)
XJM = X(J-1)
XJM = RANGE(XJM)
GO TO (150, 160, 160, 150), MODE
150 B2(JB) = R*R
B1(JB) = -(R+R)*COS(XJM*PI)
U(JP) = XJM*WS
V(JP) = R
GO TO 170
160 B2(JB) = R*R
B1(JB) = -(R+R)*COS(XJM*PI)
A1(JA) = B1(JB)/B2(JB)
A2(JA) = 1.0/B2(JB)
CONST = CONST*B2(JB)
U(JP) = XJM*WS
V(JP) = R
U0(JZ) = XJM*WS
V0(JZ) = 1.0/R
JA = JA + 1
JZ = JZ + 1
170 JB = JB + 1
JP = JP + 1
180 IF (J.LT.NN) GO TO 100
C
190 RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FUNCT
C SUBROUTINE TO CALCULATE THE FUNCTION TO BE MINIMIZED
C FROM THE PARAMETER VECTOR X AS A FUNCTION OF FREQUENCY
C-----------------------------------------------------------------------
C
SUBROUTINE FUNCT(N, X, F, G)
REAL NUM, NUMM, NUMP
LOGICAL DEB
COMMON /CM05/ AS(201), DS(201), WF(201), WG(201)
COMMON /CM10/ NP, NPD
COMMON /CM15/ DEB
COMMON /CM35/ FD(201), FDN(201)
COMMON /CM50/ N1, N2, N3, N4, N5, MODE
COMMON /CM55/ ALFA, BETA, NIT
COMMON /CM65/ IP, IQ
COMMON /CM70/ A(201), B(201)
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
DIMENSION X(N), G(N), GZR(20), GZU(20), GPR(20), GZCP(20),
* GZCR(20), GPCP(20), GPCR(20), GPRD(20), GPCPD(20), GPCRD(20)
C
RSPM = SQRT(SPM)
F = 0.0
C
DO 10 J=1,N
G(J) = 0.0
10 CONTINUE
XN1 = X(N-1)
XN = X(N)
NIT = NIT + 1
IF (DEB) WRITE (IOUTD,9999) NIT
9999 FORMAT (//27H NO. OF FUNCT. EVALUATIONS=, I4//)
C
DO 400 I=1,NP
C
FI = FD(I)
CS = COS(FI)
CS2 = 2.0*CS
H = XN
IF (MODE.EQ.4) H = XN1
HD = 1.0
HN = 1.0
DG = 0.0
C
NN = MAX0(N4,N5)
C
J = 0
20 J = J + 1
IF (J.GT.N1) GO TO 50
R = X(J)
GO TO (30, 190, 190, 40), MODE
30 NUM = 1. + R*R - R*CS2
IF (NUM.LT.SPM) NUM = SPM
HN = HN*NUM
GZR(J) = (R-CS)/NUM
GO TO 190
40 NUM = 1.0/R + R - CS2
IF (NUM.LT.SPM) NUM = SPM
HN = HN*NUM**2
GZR(J) = (1.0-1.0/R**2)/NUM
GO TO 190
C
50 IF (J.GT.N2) GO TO 70
GO TO (60, 190, 190, 60), MODE
60 P = X(J)*PI
NUMP = 2.0*(1.0-COS(FI+P))
NUMM = 2.0*(1.0-COS(FI-P))
IF (NUMM.LT.SPM) NUMM = SPM
IF (NUMP.LT.SPM) NUMP = SPM
HN = HN*NUMM*NUMP
GZU(J) = -SIN(FI-P)/NUMM + SIN(FI+P)/NUMP
GO TO 190
C
70 IF (J.GE.N3) GO TO 100
R = X(J+1)
P = X(J)*PI
GO TO (80, 190, 190, 90), MODE
80 RP = 1.0 + R*R
R2 = R + R
CSNM = COS(FI-P)
CSNP = COS(FI+P)
NUMM = RP - R2*CSNM
NUMP = RP - R2*CSNP
IF (NUMP.LT.SPM) NUMP = SPM
IF (NUMM.LT.SPM) NUMM = SPM
HN = HN*NUMM*NUMP
GZCP(J) = R*(-SIN(FI-P)/NUMM+SIN(FI+P)/NUMP)
GZCR(J) = (R-CSNM)/NUMM + (R-CSNP)/NUMP
J = J + 1
GO TO 190
90 CSNP = 2.0*COS(FI+P)
CSNM = 2.0*COS(FI-P)
R2 = R + 1.0/R
NUMM = R2 - CSNM
NUMP = R2 - CSNP
IF (NUMP.LT.RSPM) NUMP = RSPM
IF (NUMM.LT.RSPM) NUMM = RSPM
HN = HN*(NUMM*NUMP)**2
GZCP(J) = 2.0*(-SIN(FI-P)/NUMM+SIN(FI+P)/NUMP)
GZCR(J) = (1.0-1.0/R**2)*(1.0/NUMP+1.0/NUMM)
J = J + 1
GO TO 190
C
100 IF (J.GT.N4) GO TO 140
R = X(J)
RSQ = R**2
R2 = 2.0*R
DEN = 1.0 + RSQ - R2*CS
IF (DEN.LT.RSPM) DEN = RSPM
GO TO (110, 120, 120, 130), MODE
110 HD = HD*DEN
GPR(J) = (R-CS)/DEN
GO TO 190
120 DG = DG + 2.0*(1.0-R*CS)/DEN
GPRD(J) = 2.0*((1.0+RSQ)*CS-R2)/DEN**2
GO TO 190
130 HD = HD*DEN
GPR(J) = (R-CS)/DEN
DG = DG + (1.0-R*CS)/DEN
GPRD(J) = ((1.0+RSQ)*CS-R2)/DEN**2
GO TO 190
C
140 P = X(J)*PI
R = X(J+1)
R2 = R + R
RP = 1.0 + R*R
CSNP = COS(FI+P)
CSNM = COS(FI-P)
SNDP = SIN(FI+P)
SNDM = SIN(FI-P)
DENP = RP - R2*CSNP
DENM = RP - R2*CSNM
IF (ABS(DENP).LT.RSPM) DENP = RSPM
IF (ABS(DENM).LT.RSPM) DENM = RSPM
GO TO (150, 160, 160, 170), MODE
150 HD = HD*DENM*DENP
GPCR(J) = (R-CSNM)/DENM + (R-CSNP)/DENP
GPCP(J) = R*(-SNDM/DENM+SNDP/DENP)
GO TO 180
160 DENP2 = DENP**2
DENM2 = DENM**2
DG = DG + 2.0*((1.0-R*CSNM)/DENM+(1.0-R*CSNP)/DENP)
GPCPD(J) = R2*(2.0-RP)*(SNDM/DENM2-SNDP/DENP2)
GPCRD(J) = 2.0*((RP*CSNM-R2)/DENM2+(RP*CSNP-R2)/DENP2)
GO TO 180
170 HD = HD*DENM*DENP
GPCR(J) = (R-CSNM)/DENM + (R-CSNP)/DENP
GPCP(J) = R*(-SNDM/DENM+SNDP/DENP)
DENP2 = DENP**2
DENM2 = DENM**2
DG = DG + (1.0-R*CSNM)/DENM + (1.0-R*CSNP)/DENP
GPCPD(J) = R*(2.0-RP)*(SNDM/DENM2-SNDP/DENP2)
GPCRD(J) = (RP*CSNM-R2)/DENM2 + (RP*CSNP-R2)/DENP2
180 J = J + 1
C
190 IF (J.LT.NN) GO TO 20
XNP = FLOAT(NP)
XIP = FLOAT(IP)
XNPD = FLOAT(NPD)
XIQ = FLOAT(IQ)
GO TO (200, 210, 210, 220), MODE
200 H = H*SQRT(HN/HD)
E = H - AS(I)
IF (ABS(E).LT.SPM) E = SPM
A(I) = E*WF(I)
EP = E**IP/XNP
E2H = XIP*(EP/E)*H*WF(I)
E2P = E2H*PI
GO TO 230
C
210 ED = DG - DS(I) - XN
IF (ABS(ED).LT.SPM) ED = SPM
B(I) = ED*WG(I)
EDP = ED**IQ/XNPD
ED2H = XIQ*(EDP/ED)*WG(I)
ED2P = ED2H*PI
GO TO 230
C
220 H = H*SQRT(HN/HD)
E = H - AS(I)
IF (ABS(E).LT.SPM) E = SPM
A(I) = E*WF(I)
EP = E**IP/XNP
E2H = XIP*(EP/E)*H*WF(I)*ALFA
E2P = E2H*PI
ED = DG - DS(I) - XN
IF (ABS(ED).LT.SPM) ED = SPM
B(I) = ED*WG(I)
EDP = ED**IQ/XNPD
ED2H = XIQ*(EDP/ED)*WG(I)*BETA
ED2P = ED2H*PI
C
230 J = 0
240 J = J + 1
IF (J.GT.N1) GO TO 250
G(J) = G(J) + E2H*GZR(J)
GO TO 360
C
250 IF (J.GT.N2) GO TO 260
G(J) = G(J) + E2P*GZU(J)
GO TO 360
C
260 IF (J.GE.N3) GO TO 270
JJ = J + 1
G(J) = G(J) + E2P*GZCP(J)
G(JJ) = G(JJ) + E2H*GZCR(J)
J = JJ
GO TO 360
C
270 IF (J.GT.N4) GO TO 310
GO TO (280, 290, 290, 300), MODE
280 G(J) = G(J) - E2H*GPR(J)
GO TO 360
290 G(J) = G(J) + ED2H*GPRD(J)
GO TO 360
300 G(J) = G(J) - E2H*GPR(J) + ED2H*GPRD(J)
GO TO 360
C
310 JJ = J + 1
GO TO (320, 330, 330, 340), MODE
320 G(J) = G(J) - E2P*GPCP(J)
G(JJ) = G(JJ) - E2H*GPCR(J)
GO TO 350
330 G(J) = G(J) + ED2P*GPCPD(J)
G(JJ) = G(JJ) + ED2H*GPCRD(J)
GO TO 350
340 G(J) = G(J) - E2P*GPCP(J) + ED2P*GPCPD(J)
G(JJ) = G(JJ) - E2H*GPCR(J) + ED2H*GPCRD(J)
350 J = JJ
360 IF (J.LT.NN) GO TO 240
C
GO TO (370, 380, 380, 390), MODE
370 F = F + EP*WF(I)
G(N) = G(N) + E2H/XN
GO TO 400
380 F = F + EDP*WG(I)
G(N) = G(N) - ED2H
GO TO 400
390 F = F + EP*WF(I)*ALFA + EDP*WG(I)*BETA
G(N-1) = G(N-1) + E2H/XN1
G(N) = G(N) - ED2H
C
400 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C SUBROUTINE: FMFP
C TO FIND THE LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES
C BY THE METHOD OF FLETCHER AND POWELL
C-----------------------------------------------------------------------
C
SUBROUTINE FMFP(N, X, F, G, EPS, LIMIT, IER, H)
COMMON /CM25/ DEB, LONG
COMMON /CM30/ FV(100), XIT(100), KOUNT
COMMON /CM50/ M1, M2, M3, M4, M5, MODE
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
LOGICAL DEB, LONG
DIMENSION H(201), X(N), G(N), Y(40)
C
C COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUME
C
CALL FUNCT(N, X, F, G)
C
C RESET ITERATION COUNTER AND GENERATE IDENTITY MATRIX
C
KOUNT = 0
IER = 0
N2 = N + N
N3 = N2 + N
N31 = N3 + 1
10 K = N31
IF (DEB) WRITE (IOUTD,9982)
DO 30 J=1,N
NJ = N - J
H(K) = 1.0
IF (NJ.LE.0) GO TO 40
DO 20 L=1,NJ
KL = K + L
H(KL) = 0.0
20 CONTINUE
K = KL + 1
30 CONTINUE
C
C START ITERATION LOOP, UPDATE COUNTER AND SAVE FUNCTION VALUE
C
40 KOUNT = KOUNT + 1
FV(KOUNT) = F
XIT(KOUNT) = KOUNT
C
C SAVE FUNCTION VALUE, ARGUMENT VECTOR AND GRADIENT VECTOR
C
IF (DEB .OR. LONG) WRITE (IOUTD,9999) F, KOUNT
IF (DEB) WRITE (IOUTD,9997)
IF (DEB) WRITE (IOUTD,9998) (X(J),G(J),J=1,N)
OLDF = F
C
C DETERMINE DIRECTION VECTOR H
C
DO 70 J=1,N
K = N + J
H(K) = G(J)
K = K + N
H(K) = X(J)
T = 0.0
K = J + N3
DO 60 L=1,N
T = T - G(L)*H(K)
IF (L.GE.J) GO TO 50
K = K + N - L
GO TO 60
50 K = K + 1
60 CONTINUE
H(J) = T
70 CONTINUE
IF (DEB) WRITE (IOUTD,9981)
IF (DEB) WRITE (IOUTD,9996) (J,H(J),J=1,N)
C
C CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H
C
GNRM = 0.0
HNRM = 0.0
DY = 0.0
C
C CALCULATE DIRECTIONAL DERIVATIVE AND TEST VALUES FOR DIRECTION
C VECTOR H AND GRADIENT VECTOR G.
C
DO 80 J=1,N
HNRM = HNRM + ABS(H(J))
GNRM = GNRM + ABS(G(J))
DY = DY + H(J)*G(J)
80 CONTINUE
IF (DEB) WRITE (IOUTD,9995) HNRM, GNRM
IF (DEB) WRITE (IOUTD,9992) DY, F
C
C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTIONAL
C DERIVATIVE APPEARS TO BE POSITIVE OR ZERO OR IF H IS
C
C
IF (DY.GE.0.0 .OR. (HNRM/GNRM).LE.EPS) GO TO 500
C
C SEARCH MINIMUM ALONG DIRECTION H
C SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE
C
FY = F
HN = 1.0
IF (HNRM.GT.1.0) HN = 1.0/HNRM
EPSC = 2.*R1MACH(3)
UC = 1.0 - EPSC
C
C FIND MAX. STEP SIZE TO ENSURE STABILITY
C
IP1 = M3 + 1
IP2 = M4
IP3 = M4 + 2
IP4 = M5 + 1
C
C CONSIDER REAL POLES
C
IF (M3.EQ.M4) GO TO 100
DO 90 J=IP1,IP2
IF (ABS(X(J)+H(J)*HN).LT.UC) GO TO 90
IF (H(J).LT.0.0) HN = (-UC-X(J))/H(J)
IF (H(J).GT.0.0) HN = (UC-X(J))/H(J)
90 CONTINUE
C
C CONSIDER COMPLEX POLES
C
100 IF (M4.EQ.M5+1) GO TO 120
DO 110 J=IP3,IP4,2
R = X(J) + H(J)*HN
IF (R.GT.0. .AND. R.LT.UC) GO TO 110
IF (H(J).GT.0.0) HN = (UC-X(J))/H(J)
IF (H(J).LT.0.0) HN = -X(J)/H(J)
110 CONTINUE
120 IF (DEB) WRITE (IOUTD,9994) HN
DO 130 J=1,N
Y(J) = X(J) + HN*H(J)
130 CONTINUE
AMBDA = HN
CALL FUNCT(N, Y, F, G)
DYN = DY*HN
DYY = 0.0
DO 140 J=1,N
DYY = DYY + H(J)*G(J)
140 CONTINUE
IF (DEB) WRITE (IOUTD,9991) DYY, F
IF (DYY.NE.0.0) GO TO 160
DO 150 J=1,N
X(J) = Y(J)
150 CONTINUE
GO TO 390
160 CONTINUE
DEN = F - FY - DYN
IF (DEN.EQ.0.0) GO TO 170
ALFA = -0.5*DYN/DEN*HN
IF (DEB) WRITE (IOUTD,9990) ALFA, HN
C
C USE ESTIMATED STEP SIZE IF POSITIVE OTHERWISE USE HN
C
IF (ALFA.LE.0.0) GO TO 170
IF (F.GT.FY .AND. DYY.LT.0.) GO TO 170
IF (DYY*(ALFA-HN).GE.0.0) GO TO 170
AMBDA = ALFA
170 ALFA = 0.0
C
C SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT
C
180 DX = DY
FX = FY
C
C STEP ARGUMENT ALONG H
C FIND MAX. STEP SIZE TO ENSURE STABILITY
C CONSIDER REAL POLES FIRST
C
IF (M3.EQ.M4) GO TO 200
DO 190 J=IP1,IP2
IF (ABS(X(J)+H(J)*AMBDA).LT.UC) GO TO 190
IF (H(J).LT.0.0) AMBDA = (-UC-X(J))/H(J)
IF (H(J).GT.0.0) AMBDA = (UC-X(J))/H(J)
190 CONTINUE
C
C CONSIDER COMPLEX POLES
C
200 IF (M4.EQ.M5+1) GO TO 220
DO 210 J=IP3,IP4,2
R = X(J) + H(J)*AMBDA
IF (R.GT.0. .AND. R.LT.UC) GO TO 210
IF (H(J).GT.0.0) AMBDA = (UC-X(J))/H(J)
IF (H(J).LT.0.0) AMBDA = -X(J)/H(J)
210 CONTINUE
220 DO 230 I=1,N
X(I) = X(I) + AMBDA*H(I)
230 CONTINUE
IF (DEB) WRITE (IOUTD,9994) AMBDA
IF (DEB) WRITE (IOUTD,9993) (X(I),I=1,N)
C
C COMPUTE FUNCTION VALUE AND GRADIENT FOR NEW ARGUMENT
C
CALL FUNCT(N, X, F, G)
FY = F
C
C COMPUTE DIRECTIONAL DERIVATIVE DY FOR NEW ARGUMENT. TERMINATE
C SEARCH IF DY IS POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND
C
DY = 0.0
DO 240 I=1,N
DY = DY + H(I)*G(I)
240 CONTINUE
IF (DEB) WRITE (IOUTD,9989) DY, F
DYMN = AMIN1(-DX,DY)
OF = F
IF (DY) 250, 400, 260
C
C TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT
C A MINIMUM HAS BEEN PASSED
C
250 IF (FY.GE.FX) GO TO 260
IF (ABS(DX-DY).LE.EPSC) GO TO 520
C
C REPEAT SEARCH AND DOUBLE STEP SIZE FOR FURTHER SEARCHES
C
AMBDA = AMBDA + ALFA
ALFA = AMBDA
C
C END OF SEARCH LOOP
C TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE
C
IF (HNRM*AMBDA.LE.1.0) GO TO 180
C
C LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS
C
IER = 2
RETURN
C
C INTERPOLATE CUBICALLY IN THE INTERVAL DEFINED BY THE SEARCH
C ABOVE AND COMPUTE THE ARGUMENT X FOR WHICH THE INTERPOLATION
C POLYNOMIAL IS MINIMIZED
C
260 T = 0.0
270 CONTINUE
IF (DEB) WRITE (IOUTD,9987) FX, FY
IF (DEB) WRITE (IOUTD,9986) DX, DY
IF (AMBDA.EQ.0.0) GO TO 390
Z = 3.0*(FX-FY)/AMBDA + DX + DY
DALFA = Z*Z - DX*DY
IF (DALFA.LT.0.0) GO TO 500
W = SQRT(DALFA)
ALFA = DY - DX + W + W
IF (ALFA.EQ.0.0) GO TO 280
ALFA = (DY-Z+W)/ALFA
GO TO 290
280 ALFA = (Z+DY-W)/(Z+DX+Z+DY)
290 ALFA = ALFA*AMBDA
DO 300 I=1,N
X(I) = X(I) + (T-ALFA)*H(I)
300 CONTINUE
C
C TERMINATE IF THE VALUE OF THE ACTUAL FUNCTION AT X IS LESS
C THAN THE FUNCTION VALUES AT THE INTERVAL ENDS. OTHERWISE REDUCE
C THE INTERVAL BY CHOOSING ONE END POINT EQUAL TO X AND REPEAT
C THE INTERPOLATION. WHICH END POINT IS CHOSEN DEPENDS ON THE
C VALUE OF THE FUNCTION AND ITS GRADIENT AT X.
C
CALL FUNCT(N, X, F, G)
DALFA = 0.0
DO 310 I=1,N
DALFA = DALFA + H(I)*G(I)
310 CONTINUE
IF (DEB) WRITE (IOUTD,9988) DALFA, F
FYE = 1.0 - FY/F - EPS
FXE = 1.0 - FX/F - EPS
IF (DEB) WRITE (IOUTD,9985) FXE, FYE
IF (FXE.GT.0.0) GO TO 320
IF (FYE.LE.0.0) GO TO 390
320 DALFA = 0.0
DO 330 I=1,N
DALFA = DALFA + H(I)*G(I)
330 CONTINUE
IF (DALFA.GE.0.0) GO TO 360
IF (FXE) 350, 340, 360
340 IF (DX.EQ.DALFA) GO TO 390
350 FX = F
AMBDA = ALFA
T = AMBDA
DX = DALFA
GO TO 270
360 IF (FYE) 380, 370, 380
370 IF (DY.EQ.DALFA) GO TO 390
380 AMBDA = AMBDA - ALFA
DY = DALFA
FY = F
GO TO 260
C
C IF DIR. DER. AT X NOT << DIR. DER. AT INTERVAL ENDS, REPEAT INTERP
C
390 DER = ABS(DALFA/DYMN)
IF (DER.LE.1.0E-2 .OR. ABS(1.0-OF/F).LE.EPS) GO TO 400
OF = F
IF (DEB) WRITE (IOUTD,9983) DER, OLDF
IF (DALFA) 350, 350, 380
400 CONTINUE
C
C TERMINATE IF FUNCTION HAS NOT BEEN REDUCED DURING LAST ITERATION
C
IF (OLDF+EPS.LT.F) GO TO 500
C
C COMPUTE DIFFERENCE VECTORS OF ARGUMENT AND GRADIENT FROM
C TWO CONSECUTIVE ITERATIONS
C
DO 410 J=1,N
K = N + J
H(K) = G(J) - H(K)
K = N + K
H(K) = X(J) - H(K)
410 CONTINUE
IF (DEB) WRITE (IOUTD,9988) DALFA, F
C
C TEST LENGTH OF ARGUMENT DIFFERENCE VECTOR AND DIRECTION VECTOR
C IF AT LEAST N ITERATIONS HAVE BEEN EXECUTED. TERMINATE IF
C BOTH ARE LESS THAN EPS
C
IER = 0
IF (KOUNT.LT.N) GO TO 430
T = 0.0
DO 420 J=1,N
K = N2 + J
T = T + ABS(H(K))
420 CONTINUE
IF (HNRM.GT.EPS) GO TO 430
IF (T.LE.EPS) GO TO 540
C
C TERMINATE IF NUMBER OF ITERATIONS WOULD EXCEED LIMIT
C
430 IF (KOUNT.GE.LIMIT) GO TO 490
C
C PREPARE UPDATING OF MATRIX H
C
Z = 0.0
ALFA = Z
DO 460 J=1,N
W = 0.0
K = J + N3
DO 450 L=1,N
KL = N + L
W = W + H(KL)*H(K)
IF (L.GE.J) GO TO 440
K = K + N - L
GO TO 450
440 K = K + 1
450 CONTINUE
K = N + J
KN = K + N
ALFA = ALFA + W*H(K)
Z = Z + H(K)*H(KN)
H(J) = W
460 CONTINUE
IF (DEB) WRITE (IOUTD,9984) Z, ALFA
C
C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF RESULTS
C ARE NOT SATISFACTORY
C
IF (Z*ALFA.EQ.0.0) GO TO 10
C
C UPDATE MATRIX H
C
K = N31
DO 480 L=1,N
KL = N2 + L
DO 470 J=L,N
NJ = N2 + J
H(K) = H(K) + H(KL)*H(NJ)/Z - H(L)*H(J)/ALFA
K = K + 1
470 CONTINUE
480 CONTINUE
GO TO 40
C
C END OF ITERATION LOOP
C NO CONVERGENCE AFTER LIMIT ITERATIONS
C
490 IER = 1
RETURN
C
C RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS
C
500 DO 510 J=1,N
K = N2 + J
X(J) = H(K)
510 CONTINUE
CALL FUNCT(N, X, F, G)
C
C REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IN DERIVATIVE
C FAILS TO BE SUFFICIENTLY SMALL
C
IF (GNRM.LE.EPS) GO TO 530
C
C TEST FOR REPEATED FAILURE OF ITERATION
C
520 IF (IER.LT.0) GO TO 540
IER = -1
GO TO 10
530 IER = 0
540 RETURN
C
C FORMATS
C
9999 FORMAT (8H OLDF = , E11.4, 10X, 13HITERATION NO., I3)
9998 FORMAT (//(2E30.5))
9997 FORMAT (//31H FUNCTION AND GRADIENT VECTORS-)
9996 FORMAT (3(5H H(, I3, 4H) = , E13.6))
9995 FORMAT (//9H HNRM =, E12.5, 10X, 7HGNRM =, E12.5)
9994 FORMAT (//11H STEP SIZE=, E12.5)
9993 FORMAT (16H NEW ARG. VECTOR/(4E15.7))
9992 FORMAT (//9H DY(XO) =, E12.5, 10X, 7HF(XO) =, E12.5)
9991 FORMAT (//10H DY(XO+H)=, E12.5, 10X, 8HF(XO+H)=, E12.5)
9990 FORMAT (//6H ALFA=, E12.5, 10X, 3HHN=, E12.5)
9989 FORMAT (//12H DY(XO+A H)=, E12.5, 10X, 10HF(XO+A H)=, E12.5)
9988 FORMAT (//7H DYMIN=, E14.7, 10X, 4HFMN=, E14.7)
9987 FORMAT (//4H FX=, E14.7, 10X, 3HFY=, E14.7)
9986 FORMAT (//4H DX=, E14.7, 10X, 3HDY=, E14.7)
9985 FORMAT (//5H FXE=, E14.7, 10X, 4HFYE=, E14.7)
9984 FORMAT (//3H Z=, E14.7, 10X, 5HALFA=, E14.7)
9983 FORMAT (/5H DER=, E14.7, 10X, 6H OLDF=, E14.7)
9982 FORMAT (//38H ************** START STEEPEST DESCENT)
9981 FORMAT (//30X, 20H DIRECTION VECTOR H //)
END
C
C-----------------------------------------------------------------------
C FUNCTION: FS
C FUNCTION FOR SPECIFYING THE REQUIRED MAGNITUDE
C-----------------------------------------------------------------------
C
FUNCTION FS(FI)
COMMON /CM40/ FIK(10), NPK(10), NK
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
C
DO 60 K=2,NK
K1 = K - 1
C
IF (FI.LT.FIK(K1) .OR. FI.GE.FIK(K)) GO TO 60
C
GO TO (10, 20, 30, 40, 50), K1
10 FS = 1.0
GO TO 60
20 FS = 0.0
GO TO 60
30 FS = 0.0
GO TO 60
40 FS = 0.0
GO TO 60
50 FS = 0.0
60 CONTINUE
RETURN
END
C
C-----------------------------------------------------------------------
C FUNCTION: DGS
C FUNCTION FOR SPECIFYING THE REQUIRED GROUP DELAY
C-----------------------------------------------------------------------
C
FUNCTION DGS(FI)
COMMON /CM40/ FIK(10), NPK(10), NK
COMMON /UNIT/ IND, INP, IOUTD, IOUTP, PI, SPM
C
DO 60 K=2,NK
K1 = K - 1
C
IF (FI.LT.FIK(K1) .OR. FI.GE.FIK(K)) GO TO 60
C
GO TO (10, 20, 30, 40, 50), K1
10 DGS = 4.0
GO TO 60
20 DGS = 4.0
GO TO 60
30 DGS = 0.0
GO TO 60
40 DGS = 0.0
GO TO 60
50 DGS = 0.0
60 CONTINUE
RETURN
END
Return to Main Software Page