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