Make your own free website on Tripod.com
An Optimization Program for the Design of Digital Filter Transfer Functions

Return to Main Software Page

C
C-----------------------------------------------------------------------
C MAIN PROGRAM: OPTIMIZATION PROGRAM FOR THE DESIGN OF
C               DIGITAL FILTER TRANSFER FUNCTIONS
C AUTHORS:      M. T. DOLAN AND J. F. KAISER
C               BELL LABORATORIES, MURRAY HILL, NEW JERSEY 07974
C
C INPUT:        NORD  IS THE ORDER, MUST BE EVEN, .LE. 20
C               T     IS THE SAMPLING TIME
C               XK    IS THE GAIN FACTOR
C               X     IS THE COEFFICIENT VECTOR
C               NC    IS THE NUMBER OF COEFS TO REMAIN CONSTANT
C               NL    IS THE ARRAY OF INDICES OF FIXED COEFS
C                       IF NC=0, DO NOT INPUT NL ARRAY
C               NG    IS THE GAIN OPTION
C               NF    IS THE NUMBER OF FREQUENCIES FOR OPTIMIZATION
C               FREQ  IS THE ARRAY OF FREQUENCIES
C               N1    IS THE INDEX FOR TOP RANGE OF LOWER SPECS
C               SPEC1 IS THE ARRAY OF LOWER SPECS
C               TH1   IS THE ARRAY OF WEIGHTS FOR LOWER SPECS
C               N2    IS THE LOWER INDEX FOR RANGE OF UPPER SPECS
C               N3    IS THE UPPER INDEX FOR RANGE OF UPPER SPECS
C               SPEC2 IS THE ARRAY OF UPPER SPECS
C               TH2   IS THE ARRAY OF WEIGHTS FOR UPPER SPECS
C               NGID  IS THE OPTIMIZATION GUIDE OPTION
C               NOUT  IS THE OUTPUT OPTION
C               NN    IS THE INITIAL STEP OPTION
C               VBND  IS THE PARAMETER USED BY INITIAL STEP OPTIONS
C               NPEN  IS THE PENALTY FUNCTION CHOICE
C               KP    DETERMINES STOP CRITERION FOR MINOR ITERATIONS
C               EPS   IS A PARAMETER USED BY 2 MINOR STOP OPTIONS
C               KTIME IS A PARAMETER USED BY 2 MINOR STOP OPTIONS
C               NFINR DETERMINES STOP CRITERION FOR MAJOR ITERATIONS
C               PERR  IS THE PERCENTAGE VALUE USED FOR MAJOR STOP
C                       OPTION 1 WITH INTERIOR PENALTY FUNCTIONS
C               TAZ   IS THE SUM OF CONSTRAINTS-SQUARED USED FOR MAJOR
C                       STOP OPTION 1 WITH EXTERIOR PENALTY FUNCTION
C               NR    IS THE NUMBER OF R VALUES
C               INIR  IS THE R GUIDE OPTION - 1 TO GUIDE
C                                           - 0 FOR AUTOMATIC
C               RINI  IS THE INITIAL R VALUE SET BY ROUTINE ZR IF INIR=0
C               RCHNG IS THE R CHANGE FACTOR SET BY MAIN IF INIR=0
C-----------------------------------------------------------------------
C
      DIMENSION X(42), S(42), XB(42), H(42,42), G(42), Y(42), GB(42)
      DIMENSION SPEC1(42), SPEC2(42), TH1(42), TH2(42), FREQ(42), W(42)
      COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NQ, K
      COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
      COMMON /FGVARY/ T, XK, FREQ, W, SPEC1, SPEC2, TH1, TH2, R, NX, N,
     *    NF, N1, N2, N3
      COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
      COMMON /RVALUE/ RINI, INIR
      COMMON /WRITE/ IW
C
C IR AND IW ARE MACHINE DEPENDENT READ AND WRITE DEVICE NUMBERS.
C THE WRITE DEVICE NUMBER, IW, IS PASSED TO THE SUBROUTINES
C WHICH NEED IT THROUGH LABELLED COMMON.
C
      IR = I1MACH(1)
      IW = I1MACH(2)
C
C READ ORDER, SAMPLING TIME, AND GAIN FACTOR
C
      READ (IR,9999) NORD, T, XK
9999  FORMAT (I3, 2E16.8)
C
      N = NORD
      NW = 2*N
C
C READ THE COEFFICIENT VECTOR
C
      READ (IR,9998) (X(I),I=1,NW)
9998  FORMAT (4E16.8)
C
C READ NUMBER OF COEFS TO REMAIN CONSTANT
C
      READ (IR,9997) NC
9997  FORMAT (10I3)
      IF (NC.EQ.0) GO TO 10
C
C READ ARRAY OF INDICES OF FIXED COEFS IF NC .NE. 0
C
      READ (IR,9997) (NL(I),I=1,NC)
C
C READ GAIN OPTION AND NUMBER OF FREQS FOR OPTIMIZATION
C
  10  READ (IR,9997) NG, NF
C
C READ ARRAY OF FREQUENCIES
C
      READ (IR,9998) (FREQ(I),I=1,NF)
C
C READ INDEX FOR TOP RANGE OF LOWER SPECS
C
      READ (IR,9997) N1
C
C READ ARRAY OF LOWER SPECS
C
      READ (IR,9998) (SPEC1(I),I=1,N1)
C
C READ ARRAY OF WEIGHTS FOR LOWER SPECS
C
      READ (IR,9998) (TH1(I),I=1,N1)
C
C READ LOWER AND UPPER INDICES FOR RANGE OF UPPER SPECS
C
      READ (IR,9997) N2, N3
C
      N23 = N3 - N2 + 1
C
C READ ARRAY OF UPPER SPECS
C
      READ (IR,9998) (SPEC2(I),I=1,N23)
C
C READ ARRAY OF WEIGHTS FOR UPPER SPECS
C
      READ (IR,9998) (TH2(I),I=1,N23)
C
C NGID IS THE OPTIMIZATION GUIDE OPTION.  INPUT 0
C FOR AUTOMATIC SELECTION OF PARAMETERS AND NO FURTHER
C INPUT CARDS ARE NECESSARY
C READ OPTIMIZATION GUIDE OPTION
C
      READ (IR,9997) NGID
      IF (NGID.EQ.0) GO TO 20
C
C READ ADDITIONAL CARDS ONLY FOR NGID=1
C TO SIMPLIFY MATTERS, ALL DATA IS ENTERED.
C USERS SHOULD INPUT 0 FOR PARAMETERS THAT WILL NOT BE USED.
C
C IF NGID=1, READ OUTPUT OPTION, INITIAL STEP OPTION,
C PARAMETER FOR STEP OPTION, PENALTY FUNCTION CHOICE,
C STOP CRITERION FOR MINOR ITERATIONS, TWO PARAMETERS
C FOR STOP OPTION, STOP CRITERION FOR MAJOR ITERATIONS,
C THREE PARAMETERS FOR STOP OPTION, R-GUIDE OPTION,
C INITIAL R VALUE, AND THE R CHANGE FACTOR.
C
      READ (IR,9996) NOUT, NN, VBND, NPEN, KP, EPS, KTIME
      READ (IR,9995) NFINR, PERR, TAZ, NR, INIR, RINI, RCHNG
      GO TO 30
9996  FORMAT (2I3, E12.4, 2I3, E12.4, I3)
9995  FORMAT (I3, 2E12.4, 2I3, 2E12.4)
C
C AUTOMATIC SELECTION FOR NGID=0
C FOR CONSISTENCY, PARAMETERS NOT USED WILL BE SET TO 0
C
  20  NOUT = 0
      NN = 2
      VBND = 0.
      NPEN = 1
      KP = 2
C
C EPS AND KTIME WILL NOT BE USED
C
      EPS = 0.
      KTIME = 0
      NFINR = 2
C
C PERR AND TAZ WILL NOT BE USED
C
      PERR = 0.
      TAZ = 0.
      NR = 4
      INIR = 0
C
C RINI AND RCHNG WILL BE SET BY THE PROGRAM
C
      RINI = 0.
      RCHNG = 0.
C
  30  NX = NW + 2
      NVI = 0
      IF (INIR.EQ.0) RCHNG = 10.
      NXM1 = NX - 1
      X(NXM1) = SQRT(ABS(XK))
      IF (NG.NE.0) GO TO 40
      NC = NC + 1
      NL(NC) = NXM1
  40  TPI = 8.*ATAN(1.)
      NQ = NX
      DO 50 I=1,NF
        W(I) = TPI*FREQ(I)
  50  CONTINUE
      XK = X(NXM1)**2
      CALL ZR(X)
  60  NIT = 0
      NFE = 0
      NGE = 0
      K = KP
      OLDXNX = X(NX)
  70  CALL FLPWL
      IF (K) 100, 100, 80
  80  CALL STEP(NN)
      CALL QUAD
      IF (K) 100, 100, 90
  90  CALL FINISH
 100  CONTINUE
      IF (K) 110, 110, 70
 110  WRITE (IW,9994) R
      WRITE (IW,9993)
9994  FORMAT (26H0FINAL COEFFICIENTS FOR R=, E16.8)
9993  FORMAT (54H        A2              A1              B2            ,
     *    4H  B1)
      WRITE (IW,9998) (X(I),I=1,NW)
      XK = X(NXM1)**2
      WRITE (IW,9992) XK, X(NX), F
9992  FORMAT (6H GAIN=, E15.8, 7H ALPHA=, E15.8, 10H FUNCTION=, E15.8)
      CALL RINV(X, Y, NX)
      WRITE (IW,9991)
9991  FORMAT (22H0MODIFIED COEFFICIENTS)
      WRITE (IW,9998) (Y(I),I=1,NW)
      DXK = Y(NXM1)**2
      WRITE (IW,9990) DXK
9990  FORMAT (15H MODIFIED GAIN=, E15.8)
      WRITE (IW,9989)
9989  FORMAT (10H0GRADIENTS)
      WRITE (IW,9998) (G(I),I=1,NW)
      WRITE (IW,9986) NFE, NGE
      WRITE (IW,9988)
      WRITE (IW,9987)
9988  FORMAT (54H0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX,
     *    18HXXXXXXXXXXXXXXXXXX)
9987  FORMAT (1H0)
9986  FORMAT (1H , I5, 21H FUNCTION EVALUATIONS, I10, 14H GRADIENT EVAL,
     *    7HUATIONS)
      IF (NFINR.EQ.2) GO TO 120
      IF (ABS(X(NX)-OLDXNX)-ABS(PERR*.01*X(NX))) 160, 160, 130
 120  NR = NR - 1
      IF (NR.EQ.0) STOP
 130  GO TO (140, 140, 150), NPEN
 140  R = R/RCHNG
      GO TO 60
 150  R = R*RCHNG
      GO TO 60
 160  STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  ZR
C CALCULATES X(NX)=ALPHA AND INITIAL R
C-----------------------------------------------------------------------
C
      SUBROUTINE ZR(X)
C
      DIMENSION X(42), G(42), FF(42), BOUND1(42), BOUND2(42),
     *    DF(42,42), HM(42), HG(42,42)
      DIMENSION SPEC1(42), SPEC2(42), TH1(42), TH2(42), FREQ(42), W(42)
      COMMON /FGVARY/ T, XK, FREQ, W, SPEC1, SPEC2, TH1, TH2, R, NX, N,
     *    NF, N1, N2, N3
      COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
      COMMON /RVALUE/ RINI, INIR
      COMMON /WRITE/ IW
C
      NXM1 = NX - 1
      XK = X(NXM1)**2
      CALL CASC(N, X, W, NF, T, XK, HM, FF, 0, HG, DF)
      DO 10 J=1,N1
        IF (TH1(J).EQ.0.) GO TO 10
        BOUND1(J) = (FF(J)-SPEC1(J))/TH1(J)
  10  CONTINUE
      JJ = 0
      DO 20 J=N2,N3
        JJ = JJ + 1
        IF (TH2(JJ).EQ.0.) GO TO 20
        BOUND2(JJ) = (-FF(J)+SPEC2(JJ))/TH2(JJ)
  20  CONTINUE
      DO 30 J=1,N1
        IF (BOUND1(J).LT.0.) GO TO 100
  30  CONTINUE
      JJ = 0
      DO 40 J=N2,N3
        JJ = JJ + 1
        IF (BOUND2(JJ).LT.0.) GO TO 100
  40  CONTINUE
C
C ALL POSITIVE.  FIND MINIMUM HIT
C
      DO 50 J=1,N1
        IF (BOUND1(J).EQ.0.) GO TO 50
        X(NX) = BOUND1(J)
        GO TO 70
  50  CONTINUE
      JJ = 0
      DO 60 J=N2,N3
        JJ = JJ + 1
        IF (BOUND2(JJ).EQ.0.) GO TO 60
        X(NX) = BOUND2(JJ)
        GO TO 70
  60  CONTINUE
  70  DO 80 J=1,N1
        IF (BOUND1(J).EQ.0.) GO TO 80
        IF (X(NX).LE.BOUND1(J)) GO TO 80
        X(NX) = BOUND1(J)
  80  CONTINUE
      JJ = 0
      DO 90 J=N2,N3
        JJ = JJ + 1
        IF (BOUND2(JJ).EQ.0.) GO TO 90
        IF (X(NX).LE.BOUND2(JJ)) GO TO 90
        X(NX) = BOUND2(JJ)
  90  CONTINUE
      X(NX) = -X(NX)*.9999
      ZETA = X(NX)
      GO TO 150
C
C AT LEAST ONE BOUND IS NEGATIVE.  FIND THE MAXIMUM MISS.
C
 100  X(NX) = 0.
      DO 120 J=1,N1
        IF (BOUND1(J)) 110, 110, 120
 110    TXNX = ABS(BOUND1(J))
        IF (TXNX.LT.X(NX)) GO TO 120
        X(NX) = TXNX
 120  CONTINUE
      JJ = 0
      DO 140 J=N2,N3
        JJ = JJ + 1
        IF (BOUND2(JJ)) 130, 130, 140
 130    TXNX = ABS(BOUND2(JJ))
        IF (TXNX.LT.X(NX)) GO TO 140
        X(NX) = TXNX
 140  CONTINUE
      X(NX) = X(NX)*1.00001
      ZETA = X(NX)
 150  AZ = 0.
      DO 220 J=1,N1
        BOUND1(J) = FF(J) + ZETA*TH1(J) - SPEC1(J)
        IF (NPEN.EQ.3) GO TO 190
        IF (BOUND1(J)) 320, 320, 160
 160    GO TO (170, 180), NPEN
 170    BOUND1(J) = 1./BOUND1(J)
        AZ = AZ + BOUND1(J)
        GO TO 220
 180    AZ = AZ - ALOG(BOUND1(J))
        GO TO 220
 190    IF (BOUND1(J)) 200, 220, 210
 200    AZ = AZ + BOUND1(J)**2
        GO TO 220
 210    BOUND1(J) = 0.
 220  CONTINUE
      JJ = 0
      DO 290 J=N2,N3
        JJ = JJ + 1
        BOUND2(JJ) = -FF(J) + ZETA*TH2(JJ) + SPEC2(JJ)
        IF (NPEN.EQ.3) GO TO 260
        IF (BOUND2(JJ)) 320, 320, 230
 230    GO TO (240, 250), NPEN
 240    BOUND2(JJ) = 1./BOUND2(JJ)
        AZ = AZ + BOUND2(JJ)
        GO TO 290
 250    AZ = AZ - ALOG(BOUND2(JJ))
        GO TO 290
 260    IF (BOUND2(JJ)) 270, 290, 280
 270    AZ = AZ + BOUND2(JJ)**2
        GO TO 290
 280    BOUND2(JJ) = 0.
 290  CONTINUE
      IF (INIR.EQ.1) GO TO 310
      IF (AZ.EQ.0.) GO TO 300
      R = ABS(ZETA/AZ)
      IF (R.EQ.0.) R = 1.
      RETURN
 300  R = 1.
      RETURN
 310  R = RINI
      RETURN
 320  WRITE (IW,9999)
9999  FORMAT (54H0TROUBLE - CHECK THAT SPECS ARE MET AT FREQUENCIES WIT,
     *    11HH 0 WEIGHTS)
      STOP
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FLPWL
C COMPUTES DIRECTION VECTOR
C-----------------------------------------------------------------------
C
      SUBROUTINE FLPWL
C
      DIMENSION X(42), S(42), XB(42), H(42,42), G(42)
      DIMENSION SIGMA(42), GB(42), Y(42), YTH(42), HY(42)
      COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NX, K
      COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
      COMMON /WRITE/ IW
C
C FIND GRADIENT USING NEWEST X-ARRAY
C
      CALL FG
      NGE = NGE + 1
      IF (IPANIC) 20, 20, 10
  10  WRITE (IW,9999)
9999  FORMAT (24H0  FATAL IPANIC IN FLPWL)
      STOP
C
C IF THIS IS THE FIRST ITERATION, SET H = UNIT MATRIX
C
  20  IF (NIT) 30, 30, 90
  30  DO 70 I=1,NX
        DO 60 J=1,NX
          IF (I-J) 40, 50, 40
  40      H(I,J) = 0.
          GO TO 60
  50      H(I,J) = 1.
  60    CONTINUE
  70  CONTINUE
C
C IN THIS CASE, SLOPE = -GRADIENT
C G NOW BECOMES THE OLD G       STORE IT IN GB
C
      DO 80 I=1,NX
        GB(I) = G(I)
        S(I) = -G(I)
  80  CONTINUE
      GO TO 210
C
C IF THIS IS NOT THE FIRST ITERATION, CONTROL REACHES HERE.
C PREPARE TO UPDATE H
C
  90  STY = 0.
      DO 100 I=1,NX
C
C SIGMA VECTOR REPRESENTS CHANGE IN THE X-ARRAY
C
        SIGMA(I) = X(I) - XB(I)
C
C Y-VECTOR REPRESENTS CHANGE IN THE GRADIENT
C
        Y(I) = G(I) - GB(I)
C
C ONCE THE NEW G-VECTOR IS USED TO FIND Y, IT BECOMES THE OLD G
C
        GB(I) = G(I)
        STY = SIGMA(I)*Y(I) + STY
 100  CONTINUE
      IF (STY) 110, 180, 110
C
C A-MATRIX = (SIGMA-VECTOR TIMES SIGMA-TRANSPOSE) DIVIDED BY
C (SIGMA-TRANSPOSE TIMES Y-VECTOR)
C B-MATRIX = -(H-MATRIX TIMES Y-VECTOR) TIMES
C (Y-TRANSPOSE TIMES H-MATRIX) ALL DIVIDED BY
C (Y-TRANSPOSE TIMES H-MATRIX TIMES Y-VECTOR)
C
 110  DO 130 I=1,NX
        HY(I) = 0.
        YTH(I) = 0.
        DO 120 J=1,NX
          HY(I) = H(I,J)*Y(J) + HY(I)
          YTH(I) = H(J,I)*Y(J) + YTH(I)
 120    CONTINUE
 130  CONTINUE
C
      YTHY = 0.
      DO 140 I=1,NX
        YTHY = YTH(I)*Y(I) + YTHY
 140  CONTINUE
      IF (YTHY) 150, 180, 150
C
C UPDATE H
C H-MATRIX = H-MATRIX + A-MATRIX + B-MATRIX
C
 150  DO 170 I=1,NX
        DO 160 J=1,NX
          H(I,J) = -HY(I)*HY(J)/YTHY + H(I,J) + (SIGMA(I)*SIGMA(J))/STY
 160    CONTINUE
 170  CONTINUE
C
C CALCULATE THE SLOPE VECTOR  S-VECTOR=-(H-MATRIX TIMES G-VECTOR)
C
 180  DO 200 I=1,NX
        S(I) = 0.
        DO 190 J=1,NX
          S(I) = -H(J,I)*G(J) + S(I)
 190    CONTINUE
 200  CONTINUE
C
 210  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  STEP
C SETS INITIAL STEP
C-----------------------------------------------------------------------
C
      SUBROUTINE STEP(N)
C
      DIMENSION X(42), S(42), XB(42), H(42,42), G(42), GB(42)
      COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NX, K
      COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
      COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
      COMMON /WRITE/ IW
C
      IF (N-2) 10, 120, 160
C
C OPTION 1
C USER PROVIDES VMAX OR VMAX=.05*(MAX ABS X)
C STEP=MIN(2.,(VMAX/(MAX ABS S))
C
  10  VMAX = VBND
      IF (VMAX) 60, 20, 60
  20  XMAX = ABS(X(1))
      DO 40 I=2,NX
        IF (ABS(X(I))-XMAX) 40, 40, 30
  30    XMAX = ABS(X(I))
  40  CONTINUE
      VMAX = .05*XMAX
      IF (VMAX) 60, 50, 60
  50  VMAX = .01
  60  SMAX = ABS(S(1))
      DO 80 I=2,NX
        IF (ABS(S(I))-SMAX) 80, 80, 70
  70    SMAX = ABS(S(I))
  80  CONTINUE
      IF (SMAX) 90, 90, 110
  90  DEL = 2.
 100  IF (NOUT.EQ.0) RETURN
      WRITE (IW,9999) DEL
9999  FORMAT (14H0INITIAL STEP=, E16.8)
      RETURN
 110  FAC2 = VMAX/SMAX
      DEL = AMIN1(2.,FAC2)
      GO TO 100
C
C OPTION 2
C USER PROVIDES ESTIMATE OF FUNCTION MINIMUM CALLED VMIN.
C STEP=MIN(2.,-2.*(F-VMIN))/(G-TRANSPOSE*S-VECTOR)
C
 120  VMIN = VBND
      DEN = 0.
      DO 130 I=1,NX
        DEN = G(I)*S(I) + DEN
 130  CONTINUE
      IF (DEN) 140, 90, 140
 140  FAC2 = -2.*(F-VMIN)/DEN
      IF (FAC2) 90, 90, 150
 150  DEL = AMIN1(2.,FAC2)
      GO TO 100
C
C OPTION 3  -  IF THIS IS THE FIRST TIME THROUGH, USE OPTION 1.
C IF THIS IS NOT THE FIRST TIME, DIVIDE PRESENT STEP BY 2.
C
C OPTION 4 -  DO NOTHING        USER PROVIDED STEP
C
 160  IF (N-4) 170, 190, 190
 170  IF (NIT) 180, 10, 180
 180  DEL = DEL/2.
      GO TO 100
 190  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  QUAD
C PERFORMS QUADRATIC INTERPOLATION
C-----------------------------------------------------------------------
C
      SUBROUTINE QUAD
C
      DIMENSION X(42), S(42), XB(42), H(42,42), G(42), GB(42)
      DIMENSION X0(42), X1(42), X2(42)
      COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NX, K
      COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
      COMMON /WRITE/ IW
C
C EVALUATE THE FUNCTION WITH THE GIVEN X-ARRAY AND CALL IT F0.
C AT THE SAME TIME, SAVE THIS ARRAY IN XB AS THE OLD X-ARRAY
C AND SAVE ITS FTN VALUE IN FB AS THE OLD FTN VALUE.  LATER
C WE WILL STORE THE NEW ARRAY IN X AND ITS FTN VALUE IN F.
C
      CALL FUNC(X, F)
      NFE = NFE + 1
      DO 10 I=1,NX
        XB(I) = X(I)
        X0(I) = X(I)
  10  CONTINUE
      FB = F
      F0 = F
C
C GENERATE X1-ARRAY AND F1
C
  20  DO 30 I=1,NX
        X1(I) = X0(I) + DEL*S(I)
  30  CONTINUE
      CALL FUNC(X1, F1)
      NFE = NFE + 1
      IF (IPANIC) 50, 50, 40
  40  DEL = DEL/5.
      GO TO 20
C
C CHECK IF F1 IS SMALLER THAN F0
C
  50  IF (F1-F0) 60, 60, 120
  60  DO 70 I=1,NX
        X2(I) = X1(I) + DEL*S(I)
  70  CONTINUE
      CALL FUNC(X2, F2)
      NFE = NFE + 1
      IF (IPANIC) 90, 90, 80
  80  DEL = DEL/5.
      GO TO 60
  90  IF (F2-F1) 100, 170, 170
C
C IF F2 IS SMALLER THAN F1, DISCARD F0, SHIFT ALL RESULTS LEFT,
C DOUBLE THE INCREMENT, AND GENERATE A NEW X2-ARRAY AND F2
C
 100  DO 110 I=1,NX
        X0(I) = X1(I)
        X1(I) = X2(I)
 110  CONTINUE
      F0 = F1
      F1 = F2
      DEL = 2.*DEL
      GO TO 60
C
C IF CONTROL REACHES HERE, F1 WAS LARGER AND F0.
C CHECK IF IT WAS MORE THAN 2.*F0
C
 120  TF0 = 2.*F0
      IF (F1-TF0) 130, 160, 160
 130  DO 140 I=1,NX
        X2(I) = X1(I)
 140  CONTINUE
      F2 = F1
      DEL = DEL/2.
      DO 150 I=1,NX
        X1(I) = X0(I) + DEL*S(I)
 150  CONTINUE
      CALL FUNC(X1, F1)
      NFE = NFE + 1
      IF (F1-F0) 170, 170, 130
C
C F1 WAS MORE THAN 2.*F0, DIVIDE INCREMENT BY 7 AND
C GENERATE NEW X1-ARRAY AND F2
C
 160  DEL = DEL/7.
      GO TO 20
C
C INTERPOLATE
 170  DO 190 I=1,NX
        IF (S(I)) 180, 190, 180
 180    B = (X1(I)-X0(I))/S(I)
        GO TO 200
 190  CONTINUE
      GO TO 270
 200  C = B + DEL
      B2 = B**2
      C2 = C**2
      FAC1 = .5*(F0*(C2-B2)-C2*F1+B2*F2)
      IF (FAC1) 210, 260, 210
 210  FAC2 = F0*(C-B) - C*F1 + B*F2
      IF (FAC2) 220, 260, 220
 220  AL = .5*(F0*(C2-B2)-C2*F1+B2*F2)/(F0*(C-B)-C*F1+B*F2)
      DO 230 I=1,NX
        X0(I) = X0(I) + AL*S(I)
 230  CONTINUE
      CALL FUNC(X0, F0)
      NFE = NFE + 1
C
C COMPARE RESULT WITH F1 AND RETURN THE BETTER ONE
C
      IF (F1-F0) 270, 240, 240
 240  DO 250 I=1,NX
        X(I) = X0(I)
 250  CONTINUE
      F = F0
      GO TO 290
 260  K = 0
      WRITE (IW,9999)
9999  FORMAT (31H0UNDERFLOW DURING INTERPOLATION)
 270  DO 280 I=1,NX
        X(I) = X1(I)
 280  CONTINUE
      F = F1
 290  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FUNC
C EVALUATES THE PENALTY FUNCTION
C-----------------------------------------------------------------------
C
      SUBROUTINE FUNC(X, F)
C
      DIMENSION SPEC1(42), SPEC2(42), TH1(42), TH2(42), FREQ(42), W(42)
      DIMENSION HG(42,42), DF(42,42), S(42), XB(42), H(42,42), GB(42)
      DIMENSION X(42), G(42), FF(42), BOUND1(42), BOUND2(42), HM(42)
      COMMON /FGVARY/ T, XK, FREQ, W, SPEC1, SPEC2, TH1, TH2, R, NX, N,
     *    NF, N1, N2, N3
      COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
      COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
C
      NXM1 = NX - 1
      XK = X(NXM1)**2
      CALL CASC(N, X, W, NF, T, XK, HM, FF, 0, HG, DF)
      AZ = 0.
      IPANIC = 0
      ZETA = X(NX)
      DO 70 J=1,N1
        BOUND1(J) = FF(J) + ZETA*TH1(J) - SPEC1(J)
        IF (NPEN.EQ.3) GO TO 40
        IF (BOUND1(J)) 150, 150, 10
  10    GO TO (20, 30), NPEN
  20    BOUND1(J) = 1./BOUND1(J)
        AZ = AZ + BOUND1(J)
        GO TO 70
  30    AZ = AZ - ALOG(BOUND1(J))
        GO TO 70
  40    IF (BOUND1(J)) 50, 70, 60
  50    AZ = AZ + BOUND1(J)**2
        GO TO 70
  60    BOUND1(J) = 0.
  70  CONTINUE
      JJ = 0
      DO 140 J=N2,N3
        JJ = JJ + 1
        BOUND2(JJ) = -FF(J) + ZETA*TH2(JJ) + SPEC2(JJ)
        IF (NPEN.EQ.3) GO TO 110
        IF (BOUND2(JJ)) 150, 150, 80
  80    GO TO (90, 100), NPEN
  90    BOUND2(JJ) = 1./BOUND2(JJ)
        AZ = AZ + BOUND2(JJ)
        GO TO 140
 100    AZ = AZ - ALOG(BOUND2(JJ))
        GO TO 140
 110    IF (BOUND2(JJ)) 120, 140, 130
 120    AZ = AZ + BOUND2(JJ)**2
        GO TO 140
 130    BOUND2(JJ) = 0.
 140  CONTINUE
      F = ZETA + R*AZ
      GO TO 160
 150  IPANIC = 1
      NVI = NVI + 1
 160  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FG
C EVALUATES BOTH THE PENALTY FUNCTION AND ITS GRADIENT
C-----------------------------------------------------------------------
C
      SUBROUTINE FG
C
C EVALUATES THE PENALTY FUNCTION AND ITS GRADIENT FOR A FUNCTION
C WITH VALUES FF(J) AND DERIVATIVE VALUES DF(K,J) AT POINTS J FOR
C ELEMENT NUMBER K. THE ELEMENT VALUES ARE STORED AS X(K), THE
C PENALTY FUNCTION AS F, AND THE GRADIENT AS G(K).
C
      DIMENSION X(42), G(42), FF(42), BOUND1(42), BOUND2(42),
     *    DF(42,42), HM(42), HG(42,42)
      DIMENSION SPEC1(42), SPEC2(42), TH1(42), TH2(42), FREQ(42), W(42)
      DIMENSION S(42), XB(42), H(42,42), GG(42), GB(42)
      COMMON /FGVARY/ T, XK, FREQ, W, SPEC1, SPEC2, TH1, TH2, R, NX, N,
     *    NF, N1, N2, N3
      COMMON /GRAD/ TAZ, NC, NL(42), NOUT, NVI, NPEN, NFINR
      COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
      COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NQ, KDUM
      COMMON /WRITE/ IW
C
      CONS = -20./ALOG(10.)
      NXM1 = NX - 1
      NXM2 = NX - 2
      XK = X(NXM1)**2
      CALL CASC(N, X, W, NF, T, XK, HM, FF, 1, HG, DF)
      AZ = 0.
      ZETA = X(NX)
      IPANIC = 0
      DO 70 J=1,N1
        BOUND1(J) = FF(J) + ZETA*TH1(J) - SPEC1(J)
        IF (NPEN.EQ.3) GO TO 40
        IF (BOUND1(J)) 340, 340, 10
  10    GO TO (20, 30), NPEN
  20    BOUND1(J) = 1./BOUND1(J)
        AZ = AZ + BOUND1(J)
        GO TO 70
  30    AZ = AZ - ALOG(BOUND1(J))
        GO TO 70
  40    IF (BOUND1(J)) 50, 70, 60
  50    AZ = AZ + BOUND1(J)**2
        GO TO 70
  60    BOUND1(J) = 0.
  70  CONTINUE
      JJ = 0
      DO 140 J=N2,N3
        JJ = JJ + 1
        BOUND2(JJ) = -FF(J) + ZETA*TH2(JJ) + SPEC2(JJ)
        IF (NPEN.EQ.3) GO TO 110
        IF (BOUND2(JJ)) 340, 340, 80
  80    GO TO (90, 100), NPEN
  90    BOUND2(JJ) = 1./BOUND2(JJ)
        AZ = AZ + BOUND2(JJ)
        GO TO 140
 100    AZ = AZ - ALOG(BOUND2(JJ))
        GO TO 140
 110    IF (BOUND2(JJ)) 120, 140, 130
 120    AZ = AZ + BOUND2(JJ)**2
        GO TO 140
 130    BOUND2(JJ) = 0.
 140  CONTINUE
      F = ZETA + R*AZ
      IF ((NPEN.EQ.3) .AND. (NFINR.EQ.1)) GO TO 150
      GO TO 160
 150  IF (AZ.GT.TAZ) GO TO 160
      KDUM = 0
C
C THIS NEXT SECTION DETERMINES THE GRADIENT OF THE PENALTY
C FUNCTION WITH RESPECT TO THE VARIABLES X(K), BUT NOT INCLUDING
C IS X(NX).  NOR DOES IT INCLUDE XK WHICH IS X(NX-1)
C
 160  NXM2 = NX - 2
      GO TO (170, 200, 230), NPEN
 170  DO 180 J=1,N1
        BOUND1(J) = BOUND1(J)*BOUND1(J)
 180  CONTINUE
      N4 = N3 - N2 + 1
      DO 190 J=1,N4
        BOUND2(J) = BOUND2(J)*BOUND2(J)
 190  CONTINUE
      GO TO 260
 200  DO 210 J=1,N1
        BOUND1(J) = 1./BOUND1(J)
 210  CONTINUE
      N4 = N3 - N2 + 1
      DO 220 J=1,N4
        BOUND2(J) = 1./BOUND2(J)
 220  CONTINUE
      GO TO 260
C
C SIGN CHANGES BELOW SO CODING AT 2 CAN BE USED
C
 230  DO 240 J=1,N1
        BOUND1(J) = -2.*BOUND1(J)
 240  CONTINUE
      N4 = N3 - N2 + 1
      DO 250 J=1,N4
        BOUND2(J) = -2.*BOUND2(J)
 250  CONTINUE
 260  DO 290 K=1,NXM2
        AZ = 0.0
        DO 270 J=1,N1
          AZ = AZ - DF(K,J)*BOUND1(J)
 270    CONTINUE
        JJ = 0
        DO 280 J=N2,N3
          JJ = JJ + 1
          AZ = AZ + DF(K,J)*BOUND2(JJ)
 280    CONTINUE
        G(K) = R*AZ
 290  CONTINUE
C
C THE PARTIAL OF THE FTN WRT XK - SAME FOR RECIP LOG & EXT
C BECAUSE OF WHAT IS CURRENTLY STORED IN BOUNDS FOR EACH
C
      AZ = 0.
      DO 300 J=1,N1
        AZ = AZ - BOUND1(J)
 300  CONTINUE
      JJ = 0
      DO 310 J=N2,N3
        JJ = JJ + 1
        AZ = AZ + BOUND2(JJ)
 310  CONTINUE
      G(NXM1) = R*AZ*CONS*2./X(NXM1)
C
C THE GRADIENT OF THE PENALTY FUNCTION WITH RESPECT TO ZETA
C SAME CODING FOR ALL PENALTY FTNS
C
      AZ = 0.0
      DO 320 J=1,N1
        AZ = AZ + TH1(J)*BOUND1(J)
 320  CONTINUE
      JJ = 0
      DO 330 J=N2,N3
        JJ = JJ + 1
        AZ = AZ + TH2(JJ)*BOUND2(JJ)
 330  CONTINUE
      G(NX) = 1. - R*AZ
      GO TO 350
 340  IPANIC = 1
      WRITE (IW,9999)
9999  FORMAT (29H0  CONSTRAINTS VIOLATED IN FG)
 350  IF (NIT.EQ.0) GO TO 370
      IF (NOUT.EQ.0) GO TO 380
      WRITE (IW,9998) NVI
9998  FORMAT (23H0  CONSTRAINTS VIOLATED, I5, 22H TIMES DURING QUADRATI,
     *    5HC FIT)
      NVI = 0
      WRITE (IW,9997)
9997  FORMAT (13H0COEFFICIENTS)
 360  WRITE (IW,9996)
9996  FORMAT (54H        A2              A1              B2            ,
     *    4H  B1)
      WRITE (IW,9995) (X(I),I=1,NXM2)
9995  FORMAT (1X, 4E16.8)
      XK = X(NXM1)**2
      WRITE (IW,9994) XK, X(NX), F
9994  FORMAT (6H GAIN=, E15.8, 7H ALPHA=, E15.8, 10H FUNCTION=, E15.8)
      GO TO 380
 370  WRITE (IW,9993) R
9993  FORMAT (28H0INITIAL COEFFICIENTS FOR R=, E16.8)
      GO TO 360
 380  IF (NC.EQ.0) GO TO 400
      DO 390 I=1,NC
        IFIX = NL(I)
        G(IFIX) = 0.
 390  CONTINUE
 400  RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  CASC
C CALCULATES THE FREQUENCY RESPONSE AND ITS GRADIENT
C-----------------------------------------------------------------------
C
      SUBROUTINE CASC(N, AL, W, NF, T, XK, HM, HL, JJ, HG, HGL)
C
C NOTES FOR USE OUTSIDE THE PACKAGE:
C CASC CAN BE USED TO CALCULATE THE FREQUENCY RESPONSE AND ITS
C GRADIENT FOR ANY NUMBER OF FREQUENCY VALUES, NF.
C THE 42 OF THE CURRENT DIMENSION STATEMENT LIMITS THE TRANSFER
C FUNCTION TO 10 SECTIONS EACH HAVING 4 COEFFICIENTS.
C THE EXTRA 2 LOCATIONS ARE USED BY THE PACKAGE AND ARE NOT
C NECESSARY IF CASC IS USED INDEPENDENTLY.  TO USE CASC FOR
C A TRANSFER FUNCTION WITH M SECTIONS WHERE M IS GREATER
C THAN 10, CHANGE THE 42 TO 4*M AND CHANGE THE SIZE OF AH TO M.
C THERE ARE 4 COEFFICIENTS PER SECTION BECAUSE THE CONSTANTS MUST
C BE 1.
C           N   - INTEGER INPUT - ORDER OF THE FILTER .LE. 20
C           AL  - REAL INPUT - ARRAY OF COEFFICIENTS  OF THE
C                 TRANSFER FUNCTION.  SEE "DESCRIPTION OF
C                 REQUIRED DATA".
C           W   - REAL INPUT - ARRAY OF FREQUENCY VALUES MULTIPLIED
C                 BY 2 PI
C           NF  - INTEGER INPUT - NUMBER OF FREQUENCY VALUES
C           T   - REAL INPUT - SAMPLING TIME
C           XK  - REAL INPUT - GAIN FACTOR
C           HM  - REAL OUTPUT - ARRAY OF MAGNITUDE VALUES
C           HL  - REAL OUTPUT - ARRAY OF INSERTION LOSS VALUES
C           JJ  - INTEGER INPUT OPTION FOR CALCULATION OF THE
C                 GRADIENT  JJ=0 FOR NO GRADIENT
C                           JJ=1 FOR GRADIENT
C           HG  - REAL OUTPUT MATRIX - GRADIENT OF THE MAGNITUDE
C                 WHERE HG(I,J) IS THE DERIVATIVE AT FREQUENCY J
C                 FOR COEFFICIENT I
C           HGL - REAL OUTPUT MATRIX - GRADIENT OF THE INSERTION LOSS
C                 WHERE HGL(I,J) IS THE DERIVATIVE AT FREQUENCY J
C                 FOR COEFFICIENT I
C
      DIMENSION AL(42), W(NF), HM(NF), HL(NF), HG(42,NF), HGL(42,NF),
     *    AH(10)
      CONS = -20./ALOG(10.)
      NU = N/2
C
      DO 100 I=1,NF
C
C CALCULATE QUANTITIES CONSTANT WITH FREQUENCY
C
        WT = W(I)*T
        CWT = COS(WT)
        TWT = 2.*WT
        CTWT = COS(TWT)
        SWT = SIN(WT)
        STWT = SIN(TWT)
C
C INITIALIZE MAGNITUDE OF H
C
        HM(I) = 1.
C
        DO 30 K=1,NU
C
C CALCULATE QUANTITIES CONSTANT WITH ALPHA(K)
C
          K4 = 4*K
          K4M3 = K4 - 3
          K4M2 = K4 - 2
          K4M1 = K4 - 1
          A = AL(K4M3)
          B = AL(K4M2)
          C = AL(K4M1)
          D = AL(K4)
          XM1 = A*CTWT + B*CWT + 1.
          XM2 = C*CTWT + D*CWT + 1.
          XMU1 = A*STWT + B*SWT
          XMU2 = C*STWT + D*SWT
          XM12 = XM1*XM2
          XMU12 = XMU1*XMU2
          XM1MU2 = XM1*XMU2
          XM2MU1 = XM2*XMU1
          DW = XM2**2 + XMU2**2
C
C SAVE MAGNITUDE OF H(K)
          RE = (XM12+XMU12)/DW
          XIM = (XM1MU2-XM2MU1)/DW
          AH(K) = SQRT(RE**2+XIM**2)
          IF (JJ) 20, 20, 10
C
C IF NECESSARY, MAKE PARTIAL CALCULATION OF GRADIENT
C
  10      DW2 = DW**2
          CST = XM2*CTWT + XMU2*STWT
          CS = XM2*CWT + XMU2*SWT
          XXR = XM12 + XMU12
          XXI = XM1MU2 - XM2MU1
          PRA = CST/DW
          PRB = CS/DW
          PRC = (XM1*CTWT+XMU1*STWT)/DW - 2.*XXR*CST/DW2
          PRD = (XM1*CWT+XMU1*SWT)/DW - 2.*XXR*CS/DW2
          PIA = (XMU2*CTWT-XM2*STWT)/DW
          PIB = (XMU2*CWT-XM2*SWT)/DW
          PIC = (XM1*STWT-XMU1*CTWT)/DW - 2.*XXI*CST/DW2
          PID = (XM1*SWT-XMU1*CWT)/DW - 2.*XXI*CS/DW2
          HG(K4M3,I) = (RE*PRA+XIM*PIA)/AH(K)
          HG(K4M2,I) = (RE*PRB+XIM*PIB)/AH(K)
          HG(K4M1,I) = (RE*PRC+XIM*PIC)/AH(K)
          HG(K4,I) = (RE*PRD+XIM*PID)/AH(K)
C
C MULTIPLY MAGNITUDE OF H BY LATEST MAGNITUDE OF H(K)
C
  20      HM(I) = HM(I)*AH(K)
  30    CONTINUE
C
        HM(I) = HM(I)*XK
C
C FIND INSERTION LOSS
C
        HL(I) = -20.*ALOG10(HM(I))
C
        IF (JJ) 100, 100, 40
C
C IF NECESSARY, COMPLETE CALCULATION OF GRADIENT
C
  40    DO 70 K=1,NU
          K4 = 4*K
          DO 60 L=1,NU
            IF (L-K) 50, 60, 50
  50        HG(K4-3,I) = HG(K4-3,I)*AH(L)
            HG(K4-2,I) = HG(K4-2,I)*AH(L)
            HG(K4-1,I) = HG(K4-1,I)*AH(L)
            HG(K4,I) = HG(K4,I)*AH(L)
  60      CONTINUE
  70    CONTINUE
        DO 80 K=1,NU
          K4 = 4*K
          HG(K4-3,I) = HG(K4-3,I)*XK
          HG(K4-2,I) = HG(K4-2,I)*XK
          HG(K4-1,I) = HG(K4-1,I)*XK
          HG(K4,I) = HG(K4,I)*XK
  80    CONTINUE
C
C FIND GRADIENT OF INSERTION LOSS
C
        CH = CONS/HM(I)
        DO 90 K=1,NU
          K4 = 4*K
          HGL(K4-3,I) = CH*HG(K4-3,I)
          HGL(K4-2,I) = CH*HG(K4-2,I)
          HGL(K4-1,I) = CH*HG(K4-1,I)
          HGL(K4,I) = CH*HG(K4,I)
  90    CONTINUE
C
 100  CONTINUE
C
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  RINV
C INVERTS POLES OUTSIDE UNIT CIRCLE
C-----------------------------------------------------------------------
C
      SUBROUTINE RINV(X, Y, NX)
C
      DIMENSION X(42), Y(42)
C
C SQRT OF CORRECTION ON GAIN IS USED BECAUSE GAIN
C IS SQUARED BY CALLING PROGRAM. OPTIMIZER WORKS
C WITH SQRT OF GAIN TO KEEP IT POSITIVE
      NW = NX - 2
      NXM1 = NX - 1
      DO 10 I=1,NX
        Y(I) = X(I)
  10  CONTINUE
      DO 100 KK=3,NW,4
        IF (Y(KK).NE.0.) GO TO 20
        GO TO 40
  20    DISC = Y(KK+1)**2 - 4.*Y(KK)
        IF (DISC.GE.0.) GO TO 50
        IF (Y(KK).GT.1.) GO TO 30
        GO TO 100
  30    SAVE = Y(KK)
        Y(KK) = 1./Y(KK)
        Y(KK+1) = Y(KK+1)/SAVE
        Y(NXM1) = Y(NXM1)/SQRT(ABS(SAVE))
        GO TO 100
  40    IF (ABS(Y(KK+1)).LE.1.) GO TO 100
C
C FOR REAL ROOTS MULT GAIN BY ROOT BEFORE IT'S FLIPPED
C
        SAVE = Y(KK+1)
        Y(KK+1) = 1./Y(KK+1)
        Y(NXM1) = Y(NXM1)/SQRT(ABS(SAVE))
        GO TO 100
  50    R1 = (-Y(KK+1)+SQRT(DISC))/(2.*Y(KK))
        R2 = (-Y(KK+1)-SQRT(DISC))/(2.*Y(KK))
        IF (ABS(R1).LT.1.) GO TO 60
        GO TO 70
  60    Y(NXM1) = Y(NXM1)*SQRT(ABS(R1))
        R1 = 1./R1
  70    IF (ABS(R2).LT.1.) GO TO 80
        GO TO 90
  80    Y(NXM1) = Y(NXM1)*SQRT(ABS(R2))
        R2 = 1./R2
C
C RECONSTRUCT THE QUADRATIC
C
  90    Y(KK) = 1./(R1*R2)
        Y(KK+1) = -(R1+R2)/(R1*R2)
 100  CONTINUE
      RETURN
      END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FINISH
C CHECKS THE STOP CRITERION
C-----------------------------------------------------------------------
C
      SUBROUTINE FINISH
C
      COMMON X, F, DEL, VBND, EPS, KTIME, NIT, NFE, NGE, NX, N
      COMMON /SUBC/ S, XB, GB, FB, H, G, IPANIC
      DIMENSION X(42), S(42), XB(42), H(42,42), G(42), XX(42), GB(42)
C
C OPTIONS 1,2,4,5 INITIALIZE THE COUNTER KT IN OPTION 3 IN
C CASE THERE IS SWITCHING TO AND FROM OPTION 3
C
C IF N=0 ON RETURN FROM THIS SUBROUTINE, THEN THE STOP
C CRITERION HAS BEEN SATISFIED.     NIT, THE NUMBER OF
C ITERATIONS IS RESET TO ZERO.
C
C N NEGATIVE ON RETURN IS A WARNING THAT OVER 1000 ITERATIONS
C HAVE OCCURRED.
C
C N POSITIVE ON RETURN IS NORMAL AND INDICATES THAT THE
C ITERATIONS ARE TO CONTINUE.
C
      GO TO (10, 50, 60, 110, 200), N
C
C OPTION 1
C STOP WHEN GRADIENT-TRANSPOSE TIMES GRADIENT (GTG) IS LESS
C THAN SOME NUMBER GIVEN BY USER OR SET TO GTG/1.E+06
C
  10  KT = 0
      GTG = 0.
      DO 20 I=1,NX
        GTG = G(I)*G(I) + GTG
  20  CONTINUE
      IF (EPS) 30, 30, 40
  30  EPS = GTG/1.E+06
  40  IF (GTG-EPS) 210, 220, 220
C
C OPTION 2
C STOP WHEN PRESENT FUNCTION VALUE MINUS PREVIOUS FUNCTION
C VALUE IS LESS THAN .0001 TIMES PRESENT VALUE
C
  50  KT = 0
      IF (ABS(F-FB)-ABS(.0001*F)) 210, 210, 220
C
C OPTION 3
C STOP WHEN OPTION 2 IS SATISFIED A GIVEN NUMBER OF TIMES.
C IF NUMBER IS NOT SPECIFIED IT IS SET TO 4.
C
  60  IF (KTIME) 70, 70, 80
  70  KTIME = 4
  80  IF (ABS(F-FB)-ABS(.0001*F)) 100, 100, 90
  90  KT = 0
      GO TO 220
 100  KT = KT + 1
      IF (KT-KTIME) 220, 210, 210
C
C OPTION 4
C STOP WHEN MAX ABS (X-NEW MINUS X-OLD) IS LESS THAN A GIVEN
C NUMBER.  IF NOT SPECIFIED, THE NUMBER IS SET TO .0005 MAX X
C
 110  KT = 0
      DELTA = EPS
      IF (DELTA) 120, 120, 160
 120  XMAX = ABS(X(1))
      DO 140 I=2,NX
        IF (ABS(X(I))-XMAX) 140, 140, 130
 130    XMAX = ABS(X(I))
 140  CONTINUE
      DELTA = .0005*XMAX
      IF (DELTA) 150, 150, 160
 150  DELTA = .0005
 160  DO 170 I=1,NX
        XX(I) = X(I) - XB(I)
 170  CONTINUE
      DMAX = ABS(XX(1))
      DO 190 I=2,NX
        IF (ABS(XX(I))-DMAX) 190, 190, 180
 180    DMAX = ABS(XX(I))
 190  CONTINUE
      IF (DMAX-DELTA) 210, 210, 220
C
C OPTION 5
C STOP AFTER A GIVEN NUMBER OF ITERATIONS
C
 200  KT = 0
      IF (NIT-KTIME) 220, 210, 210
C
 210  N = 0
      NIT = 0
      RETURN
 220  NIT = NIT + 1
      IF (NIT-1000) 240, 230, 230
 230  N = -1
 240  RETURN
      END
Return to Main Software Page