C SUBROUTINE OPTMIZ MAXIMIZES THE LOG LIKELIHOOD FUNCTION C SUBROUTINE OPTMZ2 IS CALLED. C FUNCTION - A QUASI-NEWTON ALGORITHM FOR FINDING THE C MINIMUM OF A FUNCTION OF N VARIABLES. C USAGE - CALL OPTMIZ(FUNCT,N,NSIG,MAXFN,IOPT,X,H,G,F, C W,IER) C PARAMETERS FUNCT - A USER SUPPLIED SUBROUTINE WHICH CALCULATES C THE FUNCTION F FOR GIVEN PARAMETER VALUES C X(1),X(2),...,X(N). C THE CALLING SEQUENCE HAS THE FOLLOWING FORM C CALL FUNCT(N,X,F) C WHERE X IS A VECTOR OF LENGTH N. C FUNCT MUST APPEAR IN AN EXTERNAL STATEMENT C IN THE CALLING PROGRAM. FUNCT MUST NOT C ALTER THE VALUES OF X(I),I=1,...,N OR N. C N - THE NUMBER OF PARAMETERS (I.E., THE LENGTH C OF X) (INPUT) C NSIG - CONVERGENCE CRITERION. (INPUT). THE NUMBER C OF DIGITS OF ACCURACY REQUIRED IN THE C PARAMETER ESTIMATES. C THIS CONVERGENCE CONDITION IS SATISIFIED IF C ON TWO SUCCESSIVE ITERATIONS, THE PARAMETER C ESTIMATES (I.E.,X(I), I=1,...,N) AGREE, C COMPONENT BY COMPONENT, TO NSIG DIGITS. C MAXFN - MAXIMUM NUMBER OF FUNCTION EVALUATIONS (I.E., C CALLS TO SUBROUTINE FUNCT) ALLOWED. (INPUT) C IOPT - INPUT OPTIONS SELECTOR. C IOPT = 0 CAUSES OPTMIZ TO INITIALIZE THE C HESSIAN MATRIX H TO THE IDENTITY MATRIX. C IOPT = 1 INDICATES THAT H HAS BEEN INITIALIZED C BY THE USER TO A POSITIVE DEFINITE MATRIX. C X - VECTOR OF LENGTH N CONTAINING PARAMETER C VALUES. C ON INPUT, X MUST CONTAIN THE INITIAL C PARAMETER ESTIMATES. C ON OUTPUT, X CONTAINS THE FINAL PARAMETER C ESTIMATES AS DETERMINED BY OPTMIZ. C H - VECTOR OF LENGTH N*(N+1)/2 CONTAINING AN C ESTIMATE OF THE HESSIAN MATRIX C D**2F/(DX(I)DX(J)), I,J=1,...,N. C H IS STORED IN SYMMETRIC STORAGE MODE. C ON INPUT, IF IOPT = 0, OPTMIZ INITIALIZES H C TO THE IDENTITY MATRIX. AN INITIAL SETTING C OF H BY THE USER IS INDICATED BY IOPT=1. C H MUST BE POSITIVE DEFINITE. IF IT IS NOT, C A TERMINAL ERROR OCCURS. C ON OUTPUT, H CONTAINS AN ESTIMATE OF THE C HESSIAN AT THE FINAL PARAMETER ESTIMATES C (I.E., AT X(1),X(2),...,X(N)) C G - A VECTOR OF LENGTH N CONTAINING AN ESTIMATE C OF THE GRADIENT DF/DX(I),I=1,...,N AT THE C FINAL PARAMETER ESTIMATES. (OUTPUT) C F - A SCALAR CONTAINING THE VALUE OF THE FUNCTION C AT THE FINAL PARAMETER ESTIMATES. (OUTPUT) C W - A VECTOR OF LENGTH 3*N USED AS WORKING SPACE. C ON OUTPUT, WORK(I), CONTAINS FOR C I = 1, THE NORM OF THE GRADIENT (I.E., C SQRT(G(1)**2+G(2)**2+...+G(N)**2)) C I = 2, THE NUMBER OF FUNCTION EVALUATIONS C PERFORMED. C I = 3, AN ESTIMATE OF THE NUMBER OF C SIGNIFICANT DIGITS IN THE FINAL C PARAMETER ESTIMATES. C IER - ERROR PARAMETER. C IER = 0 IMPLIES THAT CONVERGENCE WAS C ACHIEVED AND NO ERRORS OCCURRED. C TERMINAL ERROR C IER = 129 IMPLIES THAT THE INITIAL HESSIAN C MATRIX IS NOT POSITIVE DEFINITE. THIS C CAN OCCUR ONLY FOR IOPT = 1. C IER = 130 IMPLIES THAT THE ITERATION WAS C TERMINATED DUE TO ROUNDING ERRORS C BECOMING DOMINANT. THE PARAMETER C ESTIMATES HAVE NOT BEEN DETERMINED TO C NSIG DIGITS. C IER = 131 IMPLIES THAT THE ITERATION WAS C TERMINATED BECAUSE MAXFN WAS EXCEEDED. C PRECISION - SINGLE C REQD. ROUTINES - OPTMZ2 C C ****************************************************************** C SUBROUTINE OPTMIZ (FUNCT,N,NSIG,MAXFN,IOPT,X,H,G,F,W,IER) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION X(N), G(N), H(*), W(*), REPS, AX, ZERO, ONE, HALF 1 , SEVEN, FIVE, TWELVE, TEN, HH, EPS, HJJ, V, DF, RELX, GS0, DIFF, 2 AEPS, ALPHA, FF, F, TOT, F1, F2, Z, GYS, DGS, SIG, ZZ, GNRM, P1 INTEGER N,NSIG,MAXFN,IOPT,IER,I,J,IJ,IR,NM1,NP1,JJ,L,KJ, 1 ITN,IFN,LINK,IG,JP1,JB,JNT,IDIFF,IGG,IS,K,II,IM1,NJ EXTERNAL FUNCT COMMON /ERROR/ IERR(4) DATA REPS /1.45519152284D-11/, AX /0.1D0/ DATA ZERO/0.0D0/,ONE/1.0D0/,HALF/0.5D0/,SEVEN/7.0D0/,FIVE/5.0D0/, 1 TWELVE /12.0D0/, TEN /10.0D0/, P1 /0.1D0/ C INITIALIZATION IER=0 HH=SQRT(REPS) EPS=TEN**(-NSIG) IG=N IGG=N+N IS=IGG IDIFF=1 IR=N W(1)=-ONE W(2)=ZERO W(3)=ZERO IF (IOPT.EQ.1) GO TO 30 C SET H TO THE IDENTITY MATRIX IJ=0 DO 20 I=1,N DO 10 J=1,I IJ=IJ+1 H(IJ)=ZERO 10 CONTINUE H(IJ)=ONE 20 CONTINUE GO TO 100 C FACTOR H TO L*D*L-TRANSPOSE 30 IR=N IF (N.GT.1) GO TO 40 IF (H(1).GT.ZERO) GO TO 100 H(1)=ZERO IR=0 GO TO 90 40 NM1=N-1 JJ=0 DO 80 J=1,N JJ=JJ+J HJJ=H(JJ) IF (HJJ.GT.ZERO) GO TO 50 H(JJ)=ZERO IR=IR-1 GO TO 80 50 IF (J.EQ.N) GO TO 80 IJ=JJ L=0 DO 70 I=JP1,N L=L+1 IJ=IJ+I-1 V=H(IJ)/HJJ KJ=IJ DO 60 K=I,N H(KJ+L)=H(KJ+L)-H(KJ)*V KJ=KJ+K 60 CONTINUE H(IJ)=V 70 CONTINUE 80 CONTINUE 90 IF (IR.EQ.N) GO TO 100 IER=129 GO TO 500 100 ITN=0 C EVALUATE FUNCTION AT STARTING POINT CALL FUNCT (N,X,F) IFN=1 DF=-ONE C EVALUATE GRADIENT W(IG+I),I=1,...,N 110 LINK=1 GO TO 460 120 CONTINUE C BEGIN ITERATION LOOP CALL PARPRT(ITN,X,F) IF (IFN.GE.MAXFN) GO TO 390 ITN=ITN+1 DO 130 I=1,N W(I)=-W(IG+I) 130 CONTINUE C DETERMINE SEARCH DIRECTION W C BY SOLVING H*W = -G WHERE C H = L*D*L-TRANSPOSE IF (IR.LT.N) GO TO 190 C N .EQ. 1 G(1)=W(1) IF (N.GT.1) GO TO 140 W(1)=W(1)/H(1) GO TO 190 C N .GT. 1 140 II=1 C SOLVE L*W = -G DO 160 I=2,N IJ=II II=II+I V=W(I) IM1=I-1 DO 150 J=1,IM1 IJ=IJ+1 V=V-H(IJ)*W(J) 150 CONTINUE G(I)=V W(I)=V 160 CONTINUE C SOLVE (D*LT)*Z = W WHERE C LT = L-TRANSPOSE W(N)=W(N)/H(II) JJ=II NM1=N-1 DO 180 NJ=1,NM1 C J = N-1,N-2,...,1 J=N-NJ JP1=J+1 JJ=JJ-JP1 V=W(J)/H(JJ) IJ=JJ DO 170 I=JP1,N IJ=IJ+I-1 V=V-H(IJ)*W(I) 170 CONTINUE W(J)=V 180 CONTINUE C DETERMINE STEP LENGTH ALPHA 190 RELX=ZERO GS0=ZERO DO 200 I=1,N W(IS+I)=W(I) DIFF=ABS(W(I))/MAX(ABS(X(I)),AX) RELX=MAX(RELX,DIFF) GS0=GS0+W(IG+I)*W(I) 200 CONTINUE IF (RELX.EQ.ZERO) GO TO 400 AEPS=EPS/RELX IER=130 IF (GS0.GE.ZERO) GO TO 400 IF (DF.EQ.ZERO) GO TO 400 IER=0 ALPHA=(-DF-DF)/GS0 IF (ALPHA.LE.ZERO) ALPHA=ONE ALPHA=MIN(ALPHA,ONE) IF (IDIFF.EQ.2) ALPHA=MAX(P1,ALPHA) FF=F TOT=ZERO JNT=0 C SEARCH ALONG X+ALPHA*W 210 IF (IFN.GE.MAXFN) GO TO 390 DO 220 I=1,N W(I)=X(I)+ALPHA*W(IS+I) 220 CONTINUE CALL FUNCT (N,W,F1) IFN=IFN+1 IF (F1.GE.F) GO TO 270 F2=F TOT=TOT+ALPHA 230 IER=0 F=F1 DO 240 I=1,N X(I)=W(I) 240 CONTINUE IF (JNT-1) 250,310,320 250 IF (IFN.GE.MAXFN) GO TO 390 DO 260 I=1,N W(I)=X(I)+ALPHA*W(IS+I) 260 CONTINUE CALL FUNCT (N,W,F1) IFN=IFN+1 IF (F1.GE.F) GO TO 320 IF (F1+F2.GE.F+F.AND.SEVEN*F1+FIVE*F2.GT.TWELVE*F) JNT=2 TOT=TOT+ALPHA ALPHA=ALPHA+ALPHA GO TO 230 270 CONTINUE IF (F.EQ.FF.AND.IDIFF.EQ.2.AND.RELX.GT.EPS) IER=130 IF (ALPHA.LT.AEPS) GO TO 400 IF (IFN.GE.MAXFN) GO TO 390 ALPHA=HALF*ALPHA DO 280 I=1,N W(I)=X(I)+ALPHA*W(IS+I) 280 CONTINUE CALL FUNCT (N,W,F2) IFN=IFN+1 IF (F2.GE.F) GO TO 300 TOT=TOT+ALPHA IER=0 F=F2 DO 290 I=1,N X(I)=W(I) 290 CONTINUE GO TO 310 300 Z=P1 IF (F1+F.GT.F2+F2) Z=ONE+HALF*(F-F1)/(F+F1-F2-F2) Z=MAX(P1,Z) ALPHA=Z*ALPHA JNT=1 GO TO 210 310 IF (TOT.LT.AEPS) GO TO 400 320 ALPHA=TOT C SAVE OLD GRADIENT DO 330 I=1,N W(I)=W(IG+I) 330 CONTINUE C EVALUATE GRADIENT W(IG+I), I=1,...,N LINK=2 GO TO 460 340 IF (IFN.GE.MAXFN) GO TO 390 GYS=ZERO DO 350 I=1,N GYS=GYS+W(IG+I)*W(IS+I) W(IGG+I)=W(I) 350 CONTINUE DF=FF-F DGS=GYS-GS0 IF (DGS.LE.ZERO) GO TO 120 IF (DGS+ALPHA*GS0.GT.ZERO) GO TO 370 C UPDATE HESSIAN H USING C COMPLEMENTARY DFP FORMULA SIG=ONE/GS0 IR=-IR CALL OPTMZ2 (H,N,W,SIG,G,IR,0,ZERO) DO 360 I=1,N G(I)=W(IG+I)-W(IGG+I) 360 CONTINUE SIG=ONE/(ALPHA*DGS) IR=-IR CALL OPTMZ2 (H,N,G,SIG,W,IR,0,ZERO) GO TO 120 C UPDATE HESSIAN USING C DFP FORMULA 370 ZZ=ALPHA/(DGS-ALPHA*GS0) SIG=-ZZ CALL OPTMZ2 (H,N,W,SIG,G,IR,0,REPS) Z=DGS*ZZ-ONE DO 380 I=1,N G(I)=W(IG+I)+Z*W(IGG+I) 380 CONTINUE SIG=ONE/(ZZ*DGS*DGS) CALL OPTMZ2 (H,N,G,SIG,W,IR,0,ZERO) GO TO 120 390 IER=131 C MAXFN FUNCTION EVALUATIONS GO TO 410 400 IF (IDIFF.EQ.2) GO TO 410 C CHANGE TO CENTRAL DIFFERENCES IDIFF=2 GO TO 110 410 IF (RELX.GT.EPS.AND.IER.EQ.0) GO TO 110 C MOVE GRADIENT TO G AND RETURN GNRM=ZERO DO 420 I=1,N G(I)=W(IG+I) GNRM=GNRM+G(I)*G(I) 420 CONTINUE GNRM=SQRT(GNRM) W(1)=GNRM W(2)=IFN W(3)=-LOG10(MAX(REPS,RELX)) C COMPUTE H = L*D*L-TRANSPOSE IF (N.EQ.1) GO TO 500 NP1=N+1 NM1=N-1 JJ=(N*(NP1))/2 DO 450 JB=1,NM1 JP1=NP1-JB JJ=JJ-JP1 HJJ=H(JJ) IJ=JJ L=0 DO 440 I=JP1,N L=L+1 IJ=IJ+I-1 V=H(IJ)*HJJ KJ=IJ DO 430 K=I,N H(KJ+L)=H(KJ+L)+H(KJ)*V KJ=KJ+K 430 CONTINUE H(IJ)=V 440 CONTINUE HJJ=H(JJ) 450 CONTINUE GO TO 500 C EVALUATE GRADIENT 460 IF (IDIFF.EQ.2) GO TO 480 C FORWARD DIFFERENCES C GRADIENT = W(IG+I), I=1,...,N DO 470 I=1,N Z=HH*MAX(ABS(X(I)),AX) ZZ=X(I) X(I)=ZZ+Z CALL FUNCT (N,X,F1) W(IG+I)=(F1-F)/Z X(I)=ZZ 470 CONTINUE IFN=IFN+N GO TO (120,340), LINK GO TO 120 C CENTRAL DIFFERENCES C GRADIENT = W(IG+I), I=1,...,N 480 DO 490 I=1,N Z=HH*MAX(ABS(X(I)),AX) ZZ=X(I) X(I)=ZZ+Z CALL FUNCT (N,X,F1) X(I)=ZZ-Z CALL FUNCT (N,X,F2) W(IG+I)=(F1-F2)/(Z+Z) X(I)=ZZ 490 CONTINUE IFN=IFN+N+N GO TO (120,340), LINK GO TO 120 500 CONTINUE RETURN END