C......MAIN PROGRAM                                                     DES00010
      IMPLICIT REAL*8 (A-H,O-Z)                                         DES00020
      EXTERNAL FUNCT                                                    DES00030
      DIMENSION H(184),X(16),G(16)                                      DES00040
      COMMON/RAW/W(100),Y(100),M,ICALL                                  DES00050
C......READ IN DATA                                                     DES00070
      WRITE(6,1000)                                                     DES00080
1000  FORMAT( ' '/                                                      DES00090
     &' THIS IS A DESIGN PROGRAM TO FIND LOCAL MINIMUM OF A FUNCTION'/  DES00100
     &' AND ITS COEFFICIENTS. FIRST ENTER POINTS AND AMPLITUDES IN'/    DES00110
     &' FLOATING POINT 10 FORM(EX. 0.1        1.0 8 SPACES BETWEEN #S'/ DES00120
     &' ENTER AS MANY OF THESE AS YOU WANT TO CREATE ANY FUNCTION'/     DES00130
     &' ALL POINTS< 1 EXCEPT LAST PT.=1 AMP. CAN BE >0 AND <10'/        DES00140
     &' THEN ENTER NUMBER OF VARIABLES <= 16 AND DIVISIBLE BY 4 AND'/   DES00150
     &' THE LIMIT IN INTEGER5 NOTATION(5 SPACES BETWEEN START OF #S)'/  DES00160
     &' THEN ENTER IN FLOATINGPT 10 THE EST. OF MINIMUM FUNC. VALUE'/   DES00170
     &' (.2) AND THE EXPECTED ABS. ERROR (.000000001)')                 DES00180
      WRITE(6,975)                                                      DES00190
 975  FORMAT( ' THE OUTPUT CAN BE FOUND ON DISK A IN A FILE CALLED'/    DES00200
     $' DESIGN OUTPUT A.'/                                              DES00210
     $' FILTER COEFFICIENTS CAN BE FOUND ON DISK D IN A FILE CALLED'/   DES00220
     $' DESIGN CARDS D.')                                               DES00230
      M=0                                                               DES00240
   30 M=M+1                                                             DES00250
      WRITE(22,51)                                                      DES00260
   51 FORMAT(' INPUT POINTS AND AMPLITUDES')                            DES00270
      READ(23,21,END=999)W(M),Y(M)                                      DES00280
   21 FORMAT (2F10.5)                                                   DES00290
      WRITE(8,22)M,W(M),Y(M)                                            DES00300
   22 FORMAT(' I=',I3,' W=',D15.8,' Y=',D15.8)                          DES00310
      IF(W(M).LT.1.D0)GOTO30                                            DES00320
      DO 15 J=1,16                                                      DES00330
   15 X(J)=0.D0                                                         DES00340
      X(4)=.25                                                          DES00350
      WRITE(22,900)                                                     DES00360
  900 FORMAT( ' ENTER # OF VAR., MAX. # OF ITER, EST, EPS.')            DES00370
   99 READ(23,60,END=999)N,LIMIT,EST,EPS                                DES00380
   60 FORMAT(2I5,2F10.5)                                                DES00390
      WRITE(8,61)N,LIMIT,EST,EPS                                        DES00400
   61 FORMAT(' N=',I3,' LIMIT=',I5,' EST=',D15.8,' EPS=',D15.8)         DES00410
      K=N/4                                                             DES00420
      WRITE(5,1233)K                                                    DES00430
 1233 FORMAT(' DESIGN FILTER WITH ',I3,' SECTIONS ')                    DES00440
      WRITE(5,1235)N,LIMIT,EST,EPS                                      DES00450
 1235 FORMAT('N ',I3,' LIMIT ',I5,' EST ',D15.8,' EPS ',D15.8)          DES00460
      ICALL=0                                                           DES00470
   98 CALL DFMFP(FUNCT,N,X,F,G,EST,EPS,  25 ,IER,H)                     DES00480
      CALL ROOTS(N,X)                                                   DES00490
      CALL INSIDE(N,X,KFLAG)                                            DES00500
      WRITE(8,26)IER,KFLAG,ICALL                                        DES00510
   26 FORMAT(' IER=',I5,' KFLAG=',I5,' ICALL=',I5)                      DES00520
      IF((KFLAG.NE.0.OR.IER.NE.0).AND.ICALL.LE.LIMIT)GOTO98             DES00530
      CALL ROOTS(N,X)                                                   DES00540
      ICALL=-10                                                         DES00550
      CALL FUNCT(N,X,F,G)                                               DES00560
      GOTO99                                                            DES00570
999   STOP                                                              DES00580
      END                                                               DES00590
      SUBROUTINE FUNCT(N,X,F,G)                                         DES00600
      IMPLICIT REAL*8 (A-H,O-Z)                                         DES00610
      DIMENSION H(184),X(16),G(16),YHT(100),E(100),IGRAPH(60)           DES00620
      COMPLEX*16 Z(100),NUM(100,4),DEN(100,4),Q,QBAR,ZCUR,ZCUR2         DES00630
      COMMON/RAW/W(100),Y(100),M,ICALL                                  DES00640
c     CHAR ISTAR ,IBLANK
      DATA ISTAR,IBLANK/'*',' '/                                        DES00650
      PI=-3.14159265358979                                              DES00660
      K=N/4                                                             DES00670
      IF(ICALL.NE.0)GOTO101                                             DES00680
      DO 102 I=1,M                                                      DES00690
  102 Z(I)=CDEXP(DCMPLX(0.D0,W(I)*PI))                                  DES00700
  101 A1=0.D0                                                           DES00710
      A2=0.D0                                                           DES00720
      DO 40 I=1,M                                                       DES00730
      ZCUR=Z(I)                                                         DES00740
      ZCUR2=ZCUR*ZCUR                                                   DES00750
      Q=DCMPLX(1.D0,0.D0)                                               DES00760
      DO 33 J=1,K                                                       DES00770
      J4=(J-1)*4                                                        DES00780
      NUM(I,J)=1.D0+X(J4+1)*ZCUR+X(J4+2)*ZCUR2                          DES00790
      DEN(I,J)=1.D0+X(J4+3)*ZCUR+X(J4+4)*ZCUR2                          DES00800
   33 Q=Q*NUM(I,J)/DEN(I,J)                                             DES00810
      QBAR=DCONJG(Q)                                                    DES00820
      YHT(I)=Q*QBAR                                                     DES00830
      A2=A2+YHT(I)                                                      DES00840
      YHT(I)=DSQRT(YHT(I))                                              DES00850
   40 A1=A1+YHT(I)*Y(I)                                                 DES00860
      A=A1/A2                                                           DES00870
      F=0.D0                                                            DES00880
      DO 57 J=1,16                                                      DES00890
   57 G(J)=0.D0                                                         DES00900
      DO 42 I=1,M                                                       DES00910
      ZCUR=Z(I)                                                         DES00920
      ZCUR2=ZCUR*ZCUR                                                   DES00930
      YHT(I)=A*YHT(I)                                                   DES00940
      E(I)=YHT(I)-Y(I)                                                  DES00950
      F=F+E(I)**2                                                       DES00960
      DO 42 J=1,K                                                       DES00970
      J4=(J-1)*4                                                        DES00980
      Q=2.D0*E(I)*YHT(I)/NUM(I,J)                                       DES00990
      G(J4+1)=G(J4+1)+Q*ZCUR                                            DES01000
      G(J4+2)=G(J4+2)+Q*ZCUR2                                           DES01010
      Q=-2.D0*E(I)*YHT(I)/DEN(I,J)                                      DES01020
      G(J4+3)=G(J4+3)+Q*ZCUR                                            DES01030
   42 G(J4+4)=G(J4+4)+Q*ZCUR2                                           DES01040
      ICALL=ICALL+1                                                     DES01050
      IF((ICALL/10)*10.EQ.ICALL-1)WRITE(6,25)ICALL,F                    DES01060
   25 FORMAT(' CALL NO.',I4,' F=',D15.8)                                DES01070
      IF(ICALL.GT.0)RETURN                                              DES01080
C......PRINT OUT                                                        DES01090
      WRITE(8,50)F                                                      DES01100
   50 FORMAT(' FINAL FUNCTION VALUE =',D15.8)                           DES01110
      WRITE(5,51)(X(J),J=1,N)                                           DES01120
   51 FORMAT(5X,4D15.8)                                                 DES01130
      WRITE(5,52)A                                                      DES01140
   52 FORMAT(7X,D15.8)                                                  DES01150
      WRITE(8,54)(G(J),J=1,N)                                           DES01160
   54 FORMAT(' FINAL GRADIENT ='/(' ',4D15.8))                          DES01170
      DO 55 I=1,M                                                       DES01180
   55 WRITE(8,56)I,W(I),Y(I),YHT(I),E(I)                                DES01190
   56 FORMAT(' I=',I3,' W=',D15.8,' Y=',D15.8,' YHT=',D15.8,' E=',D15.8)DES01200
      WRITE(8,59)N                                                      DES01210
   59 FORMAT(' FINAL TABLE FOR N =',I3)                                 DES01220
      DO 60 I=1,201                                                     DES01230
      FREQ=.005D0*DFLOAT(I-1)                                           DES01240
      ZCUR=CDEXP(DCMPLX(0.D0,FREQ*PI))                                  DES01250
      ZCUR2=ZCUR*ZCUR                                                   DES01260
      Q=DCMPLX(1.D0,0.D0)                                               DES01270
      PHASE=0.D0                                                        DES01280
      DO 61 J=1,K                                                       DES01290
      J4=(J-1)*4                                                        DES01300
      QBAR=1.D0+X(J4+1)*ZCUR+X(J4+2)*ZCUR2                              DES01310
      A1=QBAR                                                           DES01320
      A2=(0.D0,-1.D0)*QBAR                                              DES01330
      PHASE=PHASE+DATAN2(A2,A1)                                         DES01340
      Q=Q*QBAR                                                          DES01350
      QBAR=1.D0+X(J4+3)*ZCUR+X(J4+4)*ZCUR2                              DES01360
      A1=QBAR                                                           DES01370
      A2=(0.D0,-1.D0)*QBAR                                              DES01380
      PHASE=PHASE-DATAN2(A2,A1)                                         DES01390
   61 Q=Q/QBAR                                                          DES01400
      A1=Q*DCONJG(Q)                                                    DES01410
      A1=A*DSQRT(A1)                                                    DES01420
      PHASE=-1.D0*PHASE/PI                                              DES01430
      LOC=(A1*60.0D0)+.5D0                                              DES01440
      DO 100 JZAP=1,60                                                  DES01450
100   IGRAPH(JZAP)=IBLANK                                               DES01460
      DO 200 JFILL=1,LOC                                                DES01470
200   IGRAPH(JFILL)=ISTAR                                               DES01480
60    WRITE(8,62)FREQ,PHASE,A1,ISTAR,(IGRAPH(LLL),LLL=1,60)             DES01490
62    FORMAT( ' W=',D13.6,' PHASE/PI=',D13.6,' YHT=',D13.6,2X,A1,60A1)  DES01500
      RETURN                                                            DES01510
      END                                                               DES01520
      SUBROUTINE INSIDE(N,X,KFLAG)                                      DES01530
      IMPLICIT REAL*8 (A-H,O-Z)                                         DES01540
      DIMENSION X(16)                                                   DES01550
      J=-1                                                              DES01560
      KFLAG=0                                                           DES01570
   10 J=J+2                                                             DES01580
      IF(J.GT.N)RETURN                                                  DES01590
      B=-.5D0*X(J)                                                      DES01600
      C=X(J+1)                                                          DES01610
      DISC=B*B-C                                                        DES01620
      IF(DISC.LE.0.D0)GOTO20                                            DES01630
C......REAL ROOTS                                                       DES01640
      DISC=DSQRT(DISC)                                                  DES01650
      R1=B+DISC                                                         DES01660
      R2=B-DISC                                                         DES01670
      DR1=DABS(R1)                                                      DES01680
      DR2=DABS(R2)                                                      DES01690
      IF(DR1.LE.1.D0.AND.DR2.LE.1.D0)GOTO10                             DES01700
      KFLAG=1                                                           DES01710
      IF(DR1.GT.1.D0)R1=1.D0/R1                                         DES01720
      IF(DR2.GT.1.D0)R2=1.D0/R2                                         DES01730
      X(J)=-1.D0*(R1+R2)                                                DES01740
      X(J+1)=R1*R2                                                      DES01750
      GOTO10                                                            DES01760
C......COMPLEX ROOTS                                                    DES01770
   20 IF(C.LE.1.D0)GOTO10                                               DES01780
      KFLAG=1                                                           DES01790
      C=1.D0/C                                                          DES01800
      X(J+1)=C                                                          DES01810
      X(J)=X(J)*C                                                       DES01820
      GOTO10                                                            DES01830
      END                                                               DES01840
      SUBROUTINE ROOTS(N,X)                                             DES01850
      IMPLICIT REAL*8 (A-H,O-Z)                                         DES01860
      DIMENSION X(16)                                                   DES01870
      WRITE(8,40)                                                       DES01880
   40 FORMAT(' ROOTS'/6X,'REAL',11X,'IMAG',11X,'REAL',11X,'IMAG')       DES01890
      J=-1                                                              DES01900
   10 J=J+2                                                             DES01910
      IF(J.GT.N)RETURN                                                  DES01920
      B=-.5D0*X(J)                                                      DES01930
      C=X(J+1)                                                          DES01940
      DISC=B*B-C                                                        DES01950
      IF(DISC.LE.0.D0)GOTO20                                            DES01960
C......REAL ROOTS                                                       DES01970
      DISC=DSQRT(DISC)                                                  DES01980
      R1=B+DISC                                                         DES01990
      R2=B-DISC                                                         DES02000
      WRITE(8,30)R1,R2                                                  DES02010
   30 FORMAT(' ',D15.8,15X,D15.8)                                       DES02020
      GOTO10                                                            DES02030
C......COMPLEX ROOTS                                                    DES02040
   20 DISC=DSQRT(-1.D0*DISC)                                            DES02050
      DISCM=-1.D0*DISC                                                  DES02060
      WRITE(8,50)B,DISC,B,DISCM                                         DES02070
   50 FORMAT(' ',4D15.8)                                                DES02080
      GOTO10                                                            DES02090
      END                                                               DES02100
C                                                                       DES02110
C     ..................................................................DES02120
C                                                                       DES02130
C        SUBROUTINE DFMFP                                               DES02140
C                                                                       DES02150
C        PURPOSE                                                        DES02160
C           TO FIND A LOCAL MINIMUM OF A FUNCTION OF SEVERAL VARIABLES  DES02170
C           BY THE METHOD OF FLETCHER AND POWELL                        DES02180
C                                                                       DES02190
C        USAGE                                                          DES02200
C           CALL DFMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)               DES02210
C                                                                       DES02220
C        DESCRIPTION OF PARAMETERS                                      DES02230
C           FUNCT  - USER-WRITTEN SUBROUTINE CONCERNING THE FUNCTION TO DES02240
C                    BE MINIMIZED. IT MUST BE OF THE FORM               DES02250
C                    SUBROUTINE FUNCT(N,ARG,VAL,GRAD)                   DES02260
C                    AND MUST SERVE THE FOLLOWING PURPOSE               DES02270
C                    FOR EACH N-DIMENSIONAL ARGUMENT VECTOR  ARG,       DES02280
C                    FUNCTION VALUE AND GRADIENT VECTOR MUST BE COMPUTEDDES02290
C                    AND, ON RETURN, STORED IN VAL AND GRAD RESPECTIVELYDES02300
C                    ARG,VAL AND GRAD MUST BE OF DOUBLE PRECISION.      DES02310
C           N      - NUMBER OF VARIABLES                                DES02320
C           X      - VECTOR OF DIMENSION N CONTAINING THE INITIAL       DES02330
C                    ARGUMENT WHERE THE ITERATION STARTS. ON RETURN,    DES02340
C                    X HOLDS THE ARGUMENT CORRESPONDING TO THE          DES02350
C                    COMPUTED MINIMUM FUNCTION VALUE                    DES02360
C                    DOUBLE PRECISION VECTOR.                           DES02370
C           F      - SINGLE VARIABLE CONTAINING THE MINIMUM FUNCTION    DES02380
C                    VALUE ON RETURN, I.E. F=F(X).                      DES02390
C                    DOUBLE PRECISION VARIABLE.                         DES02400
C           G      - VECTOR OF DIMENSION N CONTAINING THE GRADIENT      DES02410
C                    VECTOR CORRESPONDING TO THE MINIMUM ON RETURN,     DES02420
C                    I.E. G=G(X).                                       DES02430
C                    DOUBLE PRECISION VECTOR.                           DES02440
C           EST    - IS AN ESTIMATE OF THE MINIMUM FUNCTION VALUE.      DES02450
C                    SINGLE PRECISION VARIABLE.                         DES02460
C           EPS    - TESTVALUE REPRESENTING THE EXPECTED ABSOLUTE ERROR.DES02470
C                    A REASONABLE CHOICE IS 10**(-16), I.E.             DES02480
C                    SOMEWHAT GREATER THAN 10**(-D), WHERE D IS THE     DES02490
C                    NUMBER OF SIGNIFICANT DIGITS IN FLOATING POINT     DES02500
C                    REPRESENTATION.                                    DES02510
C                    SINGLE PRECISION VARIABLE.                         DES02520
C           LIMIT  - MAXIMUM NUMBER OF ITERATIONS.                      DES02530
C           IER    - ERROR PARAMETER                                    DES02540
C                    IER = 0 MEANS CONVERGENCE WAS OBTAINED             DES02550
C                    IER = 1 MEANS NO CONVERGENCE IN LIMIT ITERATIONS   DES02560
C                    IER =-1 MEANS ERRORS IN GRADIENT CALCULATION       DES02570
C                    IER = 2 MEANS LINEAR SEARCH TECHNIQUE INDICATES    DES02580
C                    IT IS LIKELY THAT THERE EXISTS NO MINIMUM.         DES02590
C           H      - WORKING STORAGE OF DIMENSION N*(N+7)/2.            DES02600
C                    DOUBLE PRECISION ARRAY.                            DES02610
C                                                                       DES02620
C        REMARKS                                                        DES02630
C            I) THE SUBROUTINE NAME REPLACING THE DUMMY ARGUMENT  FUNCT DES02640
C               MUST BE DECLARED AS EXTERNAL IN THE CALLING PROGRAM.    DES02650
C           II) IER IS SET TO 2 IF , STEPPING IN ONE OF THE COMPUTED    DES02660
C               DIRECTIONS, THE FUNCTION WILL NEVER INCREASE WITHIN     DES02670
C               A TOLERABLE RANGE OF ARGUMENT.                          DES02680
C               IER = 2 MAY OCCUR ALSO IF THE INTERVAL WHERE F          DES02690
C               INCREASES IS SMALL AND THE INITIAL ARGUMENT WAS         DES02700
C               RELATIVELY FAR AWAY FROM THE MINIMUM SUCH THAT THE      DES02710
C               MINIMUM WAS OVERLEAPED. THIS IS DUE TO THE SEARCH       DES02720
C               TECHNIQUE WHICH DOUBLES THE STEPSIZE UNTIL A POINT      DES02730
C               IS FOUND WHERE THE FUNCTION INCREASES.                  DES02740
C                                                                       DES02750
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED                  DES02760
C           FUNCT                                                       DES02770
C                                                                       DES02780
C        METHOD                                                         DES02790
C           THE METHOD IS DESCRIBED IN THE FOLLOWING ARTICLE            DES02800
C           R. FLETCHER AND M.J.D. POWELL, A RAPID DESCENT METHOD FOR   DES02810
C           MINIMIZATION,                                               DES02820
C           COMPUTER JOURNAL VOL.6, ISS. 2, 1963, PP.163-168.           DES02830
C                                                                       DES02840
C     ..................................................................DES02850
C                                                                       DES02860
      SUBROUTINE DFMFP(FUNCT,N,X,F,G,EST,EPS,LIMIT,IER,H)               DES02870
C                                                                       DES02880
C        DIMENSIONED DUMMY VARIABLES                                    DES02890
      DIMENSION H(1),X(1),G(1)                                          DES02900
      DOUBLE PRECISION X,F,FX,FY,OLDF,HNRM,GNRM,H,G,DX,DY,ALFA,DALFA,   DES02910
     1AMBDA,T,Z,W                                                       DES02920
C                                                                       DES02930
C        COMPUTE FUNCTION VALUE AND GRADIENT VECTOR FOR INITIAL ARGUMENTDES02940
      CALL FUNCT(N,X,F,G)                                               DES02950
C                                                                       DES02960
C        RESET ITERATION COUNTER AND GENERATE IDENTITY MATRIX           DES02970
      IER=0                                                             DES02980
      KOUNT=0                                                           DES02990
      N2=N+N                                                            DES03000
      N3=N2+N                                                           DES03010
      N31=N3+1                                                          DES03020
    1 K=N31                                                             DES03030
      DO 4 J=1,N                                                        DES03040
      H(K)=1.D0                                                         DES03050
      NJ=N-J                                                            DES03060
      IF(NJ)5,5,2                                                       DES03070
    2 DO 3 L=1,NJ                                                       DES03080
      KL=K+L                                                            DES03090
    3 H(KL)=0.D0                                                        DES03100
    4 K=KL+1                                                            DES03110
C                                                                       DES03120
C        START ITERATION LOOP                                           DES03130
    5 KOUNT=KOUNT +1                                                    DES03140
C                                                                       DES03150
C        SAVE FUNCTION VALUE, ARGUMENT VECTOR AND GRADIENT VECTOR       DES03160
      OLDF=F                                                            DES03170
      DO 9 J=1,N                                                        DES03180
      K=N+J                                                             DES03190
      H(K)=G(J)                                                         DES03200
      K=K+N                                                             DES03210
      H(K)=X(J)                                                         DES03220
C                                                                       DES03230
C        DETERMINE DIRECTION VECTOR H                                   DES03240
      K=J+N3                                                            DES03250
      T=0.D0                                                            DES03260
      DO 8 L=1,N                                                        DES03270
      T=T-G(L)*H(K)                                                     DES03280
      IF(L-J)6,7,7                                                      DES03290
    6 K=K+N-L                                                           DES03300
      GO TO 8                                                           DES03310
    7 K=K+1                                                             DES03320
    8 CONTINUE                                                          DES03330
    9 H(J)=T                                                            DES03340
C                                                                       DES03350
C        CHECK WHETHER FUNCTION WILL DECREASE STEPPING ALONG H.         DES03360
      DY=0.D0                                                           DES03370
      HNRM=0.D0                                                         DES03380
      GNRM=0.D0                                                         DES03390
C                                                                       DES03400
C        CALCULATE DIRECTIONAL DERIVATIVE AND TESTVALUES FOR DIRECTION  DES03410
C        VECTOR H AND GRADIENT VECTOR G.                                DES03420
      DO 10 J=1,N                                                       DES03430
      HNRM=HNRM+DABS(H(J))                                              DES03440
      GNRM=GNRM+DABS(G(J))                                              DES03450
   10 DY=DY+H(J)*G(J)                                                   DES03460
C                                                                       DES03470
C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTIONAL  DES03480
C        DERIVATIVE APPEARS TO BE POSITIVE OR ZERO.                     DES03490
      IF(DY)11,51,51                                                    DES03500
C                                                                       DES03510
C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DIRECTION    DES03520
C        VECTOR H IS SMALL COMPARED TO GRADIENT VECTOR G.               DES03530
   11 IF(HNRM/GNRM-EPS)51,51,12                                         DES03540
C                                                                       DES03550
C        SEARCH MINIMUM ALONG DIRECTION H                               DES03560
C                                                                       DES03570
C        SEARCH ALONG H FOR POSITIVE DIRECTIONAL DERIVATIVE             DES03580
   12 FY=F                                                              DES03590
      ALFA=2.D0*(EST-F)/DY                                              DES03600
      AMBDA=1.D0                                                        DES03610
C                                                                       DES03620
C        USE ESTIMATE FOR STEPSIZE ONLY IF IT IS POSITIVE AND LESS THAN DES03630
C        1. OTHERWISE TAKE 1. AS STEPSIZE                               DES03640
      IF(ALFA)15,15,13                                                  DES03650
   13 IF(ALFA-AMBDA)14,15,15                                            DES03660
   14 AMBDA=ALFA                                                        DES03670
   15 ALFA=0.D0                                                         DES03680
C                                                                       DES03690
C        SAVE FUNCTION AND DERIVATIVE VALUES FOR OLD ARGUMENT           DES03700
   16 FX=FY                                                             DES03710
      DX=DY                                                             DES03720
C                                                                       DES03730
C        STEP ARGUMENT ALONG H                                          DES03740
      DO 17 I=1,N                                                       DES03750
   17 X(I)=X(I)+AMBDA*H(I)                                              DES03760
C                                                                       DES03770
C        COMPUTE FUNCTION VALUE AND GRADIENT FOR NEW ARGUMENT           DES03780
      CALL FUNCT(N,X,F,G)                                               DES03790
      FY=F                                                              DES03800
C                                                                       DES03810
C        COMPUTE DIRECTIONAL DERIVATIVE DY FOR NEW ARGUMENT.  TERMINATE DES03820
C        SEARCH, IF DY IS POSITIVE. IF DY IS ZERO THE MINIMUM IS FOUND  DES03830
      DY=0.D0                                                           DES03840
      DO 18 I=1,N                                                       DES03850
   18 DY=DY+G(I)*H(I)                                                   DES03860
      IF(DY)19,36,22                                                    DES03870
C                                                                       DES03880
C        TERMINATE SEARCH ALSO IF THE FUNCTION VALUE INDICATES THAT     DES03890
C        A MINIMUM HAS BEEN PASSED                                      DES03900
   19 IF(FY-FX)20,22,22                                                 DES03910
C                                                                       DES03920
C        REPEAT SEARCH AND DOUBLE STEPSIZE FOR FURTHER SEARCHES         DES03930
   20 AMBDA=AMBDA+ALFA                                                  DES03940
      ALFA=AMBDA                                                        DES03950
C        END OF SEARCH LOOP                                             DES03960
C                                                                       DES03970
C        TERMINATE IF THE CHANGE IN ARGUMENT GETS VERY LARGE            DES03980
      IF(HNRM*AMBDA-1.D10)16,16,21                                      DES03990
C                                                                       DES04000
C        LINEAR SEARCH TECHNIQUE INDICATES THAT NO MINIMUM EXISTS       DES04010
   21 IER=2                                                             DES04020
      RETURN                                                            DES04030
C                                                                       DES04040
C        INTERPOLATE CUBICALLY IN THE INTERVAL DEFINED BY THE SEARCH    DES04050
C        ABOVE AND COMPUTE THE ARGUMENT X FOR WHICH THE INTERPOLATION   DES04060
C        POLYNOMIAL IS MINIMIZED                                        DES04070
   22 T=0.D0                                                            DES04080
   23 IF(AMBDA)24,36,24                                                 DES04090
   24 Z=3.D0*(FX-FY)/AMBDA+DX+DY                                        DES04100
      ALFA=DMAX1(DABS(Z),DABS(DX),DABS(DY))                             DES04110
      DALFA=Z/ALFA                                                      DES04120
      DALFA=DALFA*DALFA-DX/ALFA*DY/ALFA                                 DES04130
      IF(DALFA)51,25,25                                                 DES04140
   25 W=ALFA*DSQRT(DALFA)                                               DES04150
      ALFA=DY-DX+W+W                                                    DES04160
      IF (ALFA) 250,251,250                                             DES04170
  250 ALFA=(DY-Z+W)/ALFA                                                DES04180
      GO TO 252                                                         DES04190
 251  ALFA=(Z+DY-W)/(Z+DX+Z+DY)                                         DES04200
 252  ALFA=ALFA*AMBDA                                                   DES04210
      DO 26 I=1,N                                                       DES04220
   26 X(I)=X(I)+(T-ALFA)*H(I)                                           DES04230
C                                                                       DES04240
C        TERMINATE, IF THE VALUE OF THE ACTUAL FUNCTION AT X IS LESS    DES04250
C        THAN THE FUNCTION VALUES AT THE INTERVAL ENDS. OTHERWISE REDUCEDES04260
C        THE INTERVAL BY CHOOSING ONE END-POINT EQUAL TO X AND REPEAT   DES04270
C        THE INTERPOLATION.  WHICH END-POINT IS CHOOSEN DEPENDS ON THE  DES04280
C        VALUE OF THE FUNCTION AND ITS GRADIENT AT X                    DES04290
C                                                                       DES04300
      CALL FUNCT(N,X,F,G)                                               DES04310
      IF(F-FX)27,27,28                                                  DES04320
   27 IF(F-FY)36,36,28                                                  DES04330
   28 DALFA=0.D0                                                        DES04340
      DO 29 I=1,N                                                       DES04350
   29 DALFA=DALFA+G(I)*H(I)                                             DES04360
      IF(DALFA)30,33,33                                                 DES04370
   30 IF(F-FX)32,31,33                                                  DES04380
   31 IF(DX-DALFA)32,36,32                                              DES04390
   32 FX=F                                                              DES04400
      DX=DALFA                                                          DES04410
      T=ALFA                                                            DES04420
      AMBDA=ALFA                                                        DES04430
      GO TO 23                                                          DES04440
   33 IF(FY-F)35,34,35                                                  DES04450
   34 IF(DY-DALFA)35,36,35                                              DES04460
   35 FY=F                                                              DES04470
      DY=DALFA                                                          DES04480
      AMBDA=AMBDA-ALFA                                                  DES04490
      GO TO 22                                                          DES04500
C                                                                       DES04510
C     TERMINATE, IF FUNCTION HAS NOT DECREASED DURING LAST ITERATION    DES04520
   36 IF(OLDF-F+EPS)51,38,38                                            DES04530
C                                                                       DES04540
C     COMPUTE DIFFERENCE VECTORS OF ARGUMENT AND GRADIENT FROM          DES04550
C       TWO CONSECUTIVE ITERATIONS                                      DES04560
   38 DO 37 J=1,N                                                       DES04570
      K=N+J                                                             DES04580
      H(K)=G(J)-H(K)                                                    DES04590
      K=N+K                                                             DES04600
   37 H(K)=X(J)-H(K)                                                    DES04610
C                                                                       DES04620
C        TEST LENGTH OF ARGUMENT DIFFERENCE VECTOR AND DIRECTION VECTOR DES04630
C        IF AT LEAST N ITERATIONS HAVE BEEN EXECUTED. TERMINATE, IF     DES04640
C        BOTH ARE LESS THAN  EPS                                        DES04650
      IER=0                                                             DES04660
      IF(KOUNT-N)42,39,39                                               DES04670
   39 T=0.D0                                                            DES04680
      Z=0.D0                                                            DES04690
      DO 40 J=1,N                                                       DES04700
      K=N+J                                                             DES04710
      W=H(K)                                                            DES04720
      K=K+N                                                             DES04730
      T=T+DABS(H(K))                                                    DES04740
   40 Z=Z+W*H(K)                                                        DES04750
      IF(HNRM-EPS)41,41,42                                              DES04760
   41 IF(T-EPS)56,56,42                                                 DES04770
C                                                                       DES04780
C        TERMINATE, IF NUMBER OF ITERATIONS WOULD EXCEED  LIMIT         DES04790
   42 IF(KOUNT-LIMIT)43,50,50                                           DES04800
C                                                                       DES04810
C        PREPARE UPDATING OF MATRIX H                                   DES04820
   43 ALFA=0.D0                                                         DES04830
      DO 47 J=1,N                                                       DES04840
      K=J+N3                                                            DES04850
      W=0.D0                                                            DES04860
      DO 46 L=1,N                                                       DES04870
      KL=N+L                                                            DES04880
      W=W+H(KL)*H(K)                                                    DES04890
      IF(L-J)44,45,45                                                   DES04900
   44 K=K+N-L                                                           DES04910
      GO TO 46                                                          DES04920
   45 K=K+1                                                             DES04930
   46 CONTINUE                                                          DES04940
      K=N+J                                                             DES04950
      ALFA=ALFA+W*H(K)                                                  DES04960
   47 H(J)=W                                                            DES04970
C                                                                       DES04980
C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF RESULTS      DES04990
C        ARE NOT SATISFACTORY                                           DES05000
      IF(Z*ALFA)48,1,48                                                 DES05010
C                                                                       DES05020
C        UPDATE MATRIX H                                                DES05030
   48 K=N31                                                             DES05040
      DO 49 L=1,N                                                       DES05050
      KL=N2+L                                                           DES05060
      DO 49 J=L,N                                                       DES05070
      NJ=N2+J                                                           DES05080
      H(K)=H(K)+H(KL)*H(NJ)/Z-H(L)*H(J)/ALFA                            DES05090
   49 K=K+1                                                             DES05100
      GO TO 5                                                           DES05110
C        END OF ITERATION LOOP                                          DES05120
C                                                                       DES05130
C        NO CONVERGENCE AFTER  LIMIT  ITERATIONS                        DES05140
   50 IER=1                                                             DES05150
      RETURN                                                            DES05160
C                                                                       DES05170
C        RESTORE OLD VALUES OF FUNCTION AND ARGUMENTS                   DES05180
   51 DO 52 J=1,N                                                       DES05190
      K=N2+J                                                            DES05200
   52 X(J)=H(K)                                                         DES05210
      CALL FUNCT(N,X,F,G)                                               DES05220
C                                                                       DES05230
C        REPEAT SEARCH IN DIRECTION OF STEEPEST DESCENT IF DERIVATIVE   DES05240
C        FAILS TO BE SUFFICIENTLY SMALL                                 DES05250
      IF(GNRM-EPS)55,55,53                                              DES05260
C                                                                       DES05270
C        TEST FOR REPEATED FAILURE OF ITERATION                         DES05280
   53 IF(IER)56,54,54                                                   DES05290
   54 IER=-1                                                            DES05300
      GOTO 1                                                            DES05310
   55 IER=0                                                             DES05320
   56 RETURN                                                            DES05330
      END                                                               DES05340
