SUBROUTINE DRMNG(D, FX, G, IV, LIV, LV, N, V, X) C C *** CARRY OUT DMNG (UNCONSTRAINED MINIMIZATION) ITERATIONS, USING C *** DOUBLE-DOGLEG/BFGS STEPS. C C *** PARAMETER DECLARATIONS *** C INTEGER LIV, LV, N INTEGER IV(LIV) DOUBLE PRECISION D(N), FX, G(N), V(LV), X(N) C C-------------------------- PARAMETER USAGE -------------------------- C C D.... SCALE VECTOR. C FX... FUNCTION VALUE. C G.... GRADIENT VECTOR. C IV... INTEGER VALUE ARRAY. C LIV.. LENGTH OF IV (AT LEAST 60). C LV... LENGTH OF V (AT LEAST 71 + N*(N+13)/2). C N.... NUMBER OF VARIABLES (COMPONENTS IN X AND G). C V.... FLOATING-POINT VALUE ARRAY. C X.... VECTOR OF PARAMETERS TO BE OPTIMIZED. C C *** DISCUSSION *** C C PARAMETERS IV, N, V, AND X ARE THE SAME AS THE CORRESPONDING C ONES TO DMNG (WHICH SEE), EXCEPT THAT V CAN BE SHORTER (SINCE C THE PART OF V THAT DMNG USES FOR STORING G IS NOT NEEDED). C MOREOVER, COMPARED WITH DMNG, IV(1) MAY HAVE THE TWO ADDITIONAL C OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, AS IS THE USE C OF IV(TOOBIG) AND IV(NFGCAL). THE VALUE IV(G), WHICH IS AN C OUTPUT VALUE FROM DMNG (AND DMNF), IS NOT REFERENCED BY C DRMNG OR THE SUBROUTINES IT CALLS. C FX AND G NEED NOT HAVE BEEN INITIALIZED WHEN DRMNG IS CALLED C WITH IV(1) = 12, 13, OR 14. C C IV(1) = 1 MEANS THE CALLER SHOULD SET FX TO F(X), THE FUNCTION VALUE C AT X, AND CALL DRMNG AGAIN, HAVING CHANGED NONE OF THE C OTHER PARAMETERS. AN EXCEPTION OCCURS IF F(X) CANNOT BE C (E.G. IF OVERFLOW WOULD OCCUR), WHICH MAY HAPPEN BECAUSE C OF AN OVERSIZED STEP. IN THIS CASE THE CALLER SHOULD SET C IV(TOOBIG) = IV(2) TO 1, WHICH WILL CAUSE DRMNG TO IG- C NORE FX AND TRY A SMALLER STEP. THE PARAMETER NF THAT C DMNG PASSES TO CALCF (FOR POSSIBLE USE BY CALCG) IS A C COPY OF IV(NFCALL) = IV(6). C IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT VECTOR C OF F AT X, AND CALL DRMNG AGAIN, HAVING CHANGED NONE OF C THE OTHER PARAMETERS EXCEPT POSSIBLY THE SCALE VECTOR D C WHEN IV(DTYPE) = 0. THE PARAMETER NF THAT DMNG PASSES C TO CALCG IS IV(NFGCAL) = IV(7). IF G(X) CANNOT BE C EVALUATED, THEN THE CALLER MAY SET IV(TOOBIG) TO 0, IN C WHICH CASE DRMNG WILL RETURN WITH IV(1) = 65. C. C *** GENERAL *** C C CODED BY DAVID M. GAY (DECEMBER 1979). REVISED SEPT. 1982. C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED C IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS C MCS-7600324 AND MCS-7906671. C C (SEE DMNG FOR REFERENCES.) C C+++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C INTEGER DG1, DUMMY, G01, I, K, L, LSTGST, NWTST1, RSTRST, STEP1, 1 TEMP1, W, X01, Z DOUBLE PRECISION T C C *** CONSTANTS *** C DOUBLE PRECISION HALF, NEGONE, ONE, ONEP2, ZERO C C *** NO INTRINSIC FUNCTIONS *** C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C LOGICAL STOPX DOUBLE PRECISION DD7TPR, DRLDST, DV2NRM EXTERNAL DA7SST,DD7DOG,DIVSET, DD7TPR,DITSUM, DL7ITV, DL7IVM, 1 DL7TVM, DL7UPD,DL7VML,DPARCK, DRLDST, STOPX,DV2AXY, 2 DV7CPY, DV7SCP, DV7VMP, DV2NRM, DW7ZBF C C DA7SST.... ASSESSES CANDIDATE STEP. C DD7DOG.... COMPUTES DOUBLE-DOGLEG (CANDIDATE) STEP. C DIVSET.... SUPPLIES DEFAULT IV AND V INPUT COMPONENTS. C DD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. C DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. C DL7ITV... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. C DL7IVM... MULTIPLIES INVERSE OF LOWER TRIANGLE TIMES VECTOR. C DL7TVM... MULTIPLIES TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. C LUPDT.... UPDATES CHOLESKY FACTOR OF HESSIAN APPROXIMATION. C DL7VML.... MULTIPLIES LOWER TRIANGLE TIMES VECTOR. C DPARCK.... CHECKS VALIDITY OF INPUT IV AND V VALUES. C DRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. C DV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C DV7VMP... MULTIPLIES VECTOR BY VECTOR RAISED TO POWER (COMPONENTWISE). C DV2NRM... RETURNS THE 2-NORM OF A VECTOR. C DW7ZBF... COMPUTES W AND Z FOR DL7UPD CORRESPONDING TO BFGS UPDATE. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, DG, DGNORM, DINIT, DSTNRM, DST0, F, F0, FDIF, 1 GTHG, GTSTEP, G0, INCFAC, INITH, IRC, KAGQT, LMAT, LMAX0, 2 LMAXS, MODE, MODEL, MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL, 3 NGCALL, NITER, NREDUC, NWTSTP, PREDUC, RADFAC, RADINC, 4 RADIUS, RAD0, RELDX, RESTOR, STEP, STGLIM, STLSTG, TOOBIG, 5 TUNER4, TUNER5, VNEED, XIRC, X0 C C *** IV SUBSCRIPT VALUES *** C C/6 C DATA CNVCOD/55/, DG/37/, G0/48/, INITH/25/, IRC/29/, KAGQT/33/, C 1 MODE/35/, MODEL/5/, MXFCAL/17/, MXITER/18/, NFCALL/6/, C 2 NFGCAL/7/, NGCALL/30/, NITER/31/, NWTSTP/34/, RADINC/8/, C 3 RESTOR/9/, STEP/40/, STGLIM/11/, STLSTG/41/, TOOBIG/2/, C 4 VNEED/4/, XIRC/13/, X0/43/ C/7 PARAMETER (CNVCOD=55, DG=37, G0=48, INITH=25, IRC=29, KAGQT=33, 1 MODE=35, MODEL=5, MXFCAL=17, MXITER=18, NFCALL=6, 2 NFGCAL=7, NGCALL=30, NITER=31, NWTSTP=34, RADINC=8, 3 RESTOR=9, STEP=40, STGLIM=11, STLSTG=41, TOOBIG=2, 4 VNEED=4, XIRC=13, X0=43) C/ C C *** V SUBSCRIPT VALUES *** C C/6 C DATA DGNORM/1/, DINIT/38/, DSTNRM/2/, DST0/3/, F/10/, F0/13/, C 1 FDIF/11/, GTHG/44/, GTSTEP/4/, INCFAC/23/, LMAT/42/, C 2 LMAX0/35/, LMAXS/36/, NEXTV/47/, NREDUC/6/, PREDUC/7/, C 3 RADFAC/16/, RADIUS/8/, RAD0/9/, RELDX/17/, TUNER4/29/, C 4 TUNER5/30/ C/7 PARAMETER (DGNORM=1, DINIT=38, DSTNRM=2, DST0=3, F=10, F0=13, 1 FDIF=11, GTHG=44, GTSTEP=4, INCFAC=23, LMAT=42, 2 LMAX0=35, LMAXS=36, NEXTV=47, NREDUC=6, PREDUC=7, 3 RADFAC=16, RADIUS=8, RAD0=9, RELDX=17, TUNER4=29, 4 TUNER5=30) C/ C C/6 C DATA HALF/0.5D+0/, NEGONE/-1.D+0/, ONE/1.D+0/, ONEP2/1.2D+0/, C 1 ZERO/0.D+0/ C/7 PARAMETER (HALF=0.5D+0, NEGONE=-1.D+0, ONE=1.D+0, ONEP2=1.2D+0, 1 ZERO=0.D+0) C/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C I = IV(1) IF (I .EQ. 1) GO TO 50 IF (I .EQ. 2) GO TO 60 C C *** CHECK VALIDITY OF IV AND V INPUT VALUES *** C IF (IV(1) .EQ. 0) CALL DIVSET(2, IV, LIV, LV, V) IF (IV(1) .EQ. 12 .OR. IV(1) .EQ. 13) 1 IV(VNEED) = IV(VNEED) + N*(N+13)/2 CALL DPARCK(2, D, IV, LIV, LV, N, V) I = IV(1) - 2 IF (I .GT. 12) GO TO 999 GO TO (190, 190, 190, 190, 190, 190, 120, 90, 120, 10, 10, 20), I C C *** STORAGE ALLOCATION *** C 10 L = IV(LMAT) IV(X0) = L + N*(N+1)/2 IV(STEP) = IV(X0) + N IV(STLSTG) = IV(STEP) + N IV(G0) = IV(STLSTG) + N IV(NWTSTP) = IV(G0) + N IV(DG) = IV(NWTSTP) + N IV(NEXTV) = IV(DG) + N IF (IV(1) .NE. 13) GO TO 20 IV(1) = 14 GO TO 999 C C *** INITIALIZATION *** C 20 IV(NITER) = 0 IV(NFCALL) = 1 IV(NGCALL) = 1 IV(NFGCAL) = 1 IV(MODE) = -1 IV(MODEL) = 1 IV(STGLIM) = 1 IV(TOOBIG) = 0 IV(CNVCOD) = 0 IV(RADINC) = 0 V(RAD0) = ZERO IF (V(DINIT) .GE. ZERO) CALL DV7SCP(N, D, V(DINIT)) IF (IV(INITH) .NE. 1) GO TO 40 C C *** SET THE INITIAL HESSIAN APPROXIMATION TO DIAG(D)**-2 *** C L = IV(LMAT) CALL DV7SCP(N*(N+1)/2, V(L), ZERO) K = L - 1 DO 30 I = 1, N K = K + I T = D(I) IF (T .LE. ZERO) T = ONE V(K) = T 30 CONTINUE C C *** COMPUTE INITIAL FUNCTION VALUE *** C 40 IV(1) = 1 GO TO 999 C 50 V(F) = FX IF (IV(MODE) .GE. 0) GO TO 190 V(F0) = FX IV(1) = 2 IF (IV(TOOBIG) .EQ. 0) GO TO 999 IV(1) = 63 GO TO 350 C C *** MAKE SURE GRADIENT COULD BE COMPUTED *** C 60 IF (IV(TOOBIG) .EQ. 0) GO TO 70 IV(1) = 65 GO TO 350 C 70 DG1 = IV(DG) CALL DV7VMP(N, V(DG1), G, D, -1) V(DGNORM) = DV2NRM(N, V(DG1)) C IF (IV(CNVCOD) .NE. 0) GO TO 340 IF (IV(MODE) .EQ. 0) GO TO 300 C C *** ALLOW FIRST STEP TO HAVE SCALED 2-NORM AT MOST V(LMAX0) *** C V(RADIUS) = V(LMAX0) C IV(MODE) = 0 C C C----------------------------- MAIN LOOP ----------------------------- C C C *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** C 80 CALL DITSUM(D, G, IV, LIV, LV, N, V, X) 90 K = IV(NITER) IF (K .LT. IV(MXITER)) GO TO 100 IV(1) = 10 GO TO 350 C C *** UPDATE RADIUS *** C 100 IV(NITER) = K + 1 IF (K .GT. 0) V(RADIUS) = V(RADFAC) * V(DSTNRM) C C *** INITIALIZE FOR START OF NEXT ITERATION *** C G01 = IV(G0) X01 = IV(X0) V(F0) = V(F) IV(IRC) = 4 IV(KAGQT) = -1 C C *** COPY X TO X0, G TO G0 *** C CALL DV7CPY(N, V(X01), X) CALL DV7CPY(N, V(G01), G) C C *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** C 110 IF (.NOT. STOPX(DUMMY)) GO TO 130 IV(1) = 11 GO TO 140 C C *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. C 120 IF (V(F) .GE. V(F0)) GO TO 130 V(RADFAC) = ONE K = IV(NITER) GO TO 100 C 130 IF (IV(NFCALL) .LT. IV(MXFCAL)) GO TO 150 IV(1) = 9 140 IF (V(F) .GE. V(F0)) GO TO 350 C C *** IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH C *** IMPROVED V(F), EVALUATE THE GRADIENT AT X. C IV(CNVCOD) = IV(1) GO TO 290 C C. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . C 150 STEP1 = IV(STEP) DG1 = IV(DG) NWTST1 = IV(NWTSTP) IF (IV(KAGQT) .GE. 0) GO TO 160 L = IV(LMAT) CALL DL7IVM(N, V(NWTST1), V(L), G) V(NREDUC) = HALF * DD7TPR(N, V(NWTST1), V(NWTST1)) CALL DL7ITV(N, V(NWTST1), V(L), V(NWTST1)) CALL DV7VMP(N, V(STEP1), V(NWTST1), D, 1) V(DST0) = DV2NRM(N, V(STEP1)) CALL DV7VMP(N, V(DG1), V(DG1), D, -1) CALL DL7TVM(N, V(STEP1), V(L), V(DG1)) V(GTHG) = DV2NRM(N, V(STEP1)) IV(KAGQT) = 0 160 CALL DD7DOG(V(DG1), LV, N, V(NWTST1), V(STEP1), V) IF (IV(IRC) .NE. 6) GO TO 170 IF (IV(RESTOR) .NE. 2) GO TO 190 RSTRST = 2 GO TO 200 C C *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** C 170 IV(TOOBIG) = 0 IF (V(DSTNRM) .LE. ZERO) GO TO 190 IF (IV(IRC) .NE. 5) GO TO 180 IF (V(RADFAC) .LE. ONE) GO TO 180 IF (V(PREDUC) .GT. ONEP2 * V(FDIF)) GO TO 180 IF (IV(RESTOR) .NE. 2) GO TO 190 RSTRST = 0 GO TO 200 C C *** COMPUTE F(X0 + STEP) *** C 180 X01 = IV(X0) STEP1 = IV(STEP) CALL DV2AXY(N, X, ONE, V(STEP1), V(X01)) IV(NFCALL) = IV(NFCALL) + 1 IV(1) = 1 GO TO 999 C C. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . C 190 RSTRST = 3 200 X01 = IV(X0) V(RELDX) = DRLDST(N, D, X, V(X01)) CALL DA7SST(IV, LIV, LV, V) STEP1 = IV(STEP) LSTGST = IV(STLSTG) I = IV(RESTOR) + 1 GO TO (240, 210, 220, 230), I 210 CALL DV7CPY(N, X, V(X01)) GO TO 240 220 CALL DV7CPY(N, V(LSTGST), V(STEP1)) GO TO 240 230 CALL DV7CPY(N, V(STEP1), V(LSTGST)) CALL DV2AXY(N, X, ONE, V(STEP1), V(X01)) V(RELDX) = DRLDST(N, D, X, V(X01)) IV(RESTOR) = RSTRST C 240 K = IV(IRC) GO TO (250,280,280,280,250,260,270,270,270,270,270,270,330,300), K C C *** RECOMPUTE STEP WITH CHANGED RADIUS *** C 250 V(RADIUS) = V(RADFAC) * V(DSTNRM) GO TO 110 C C *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST. C 260 V(RADIUS) = V(LMAXS) GO TO 150 C C *** CONVERGENCE OR FALSE CONVERGENCE *** C 270 IV(CNVCOD) = K - 4 IF (V(F) .GE. V(F0)) GO TO 340 IF (IV(XIRC) .EQ. 14) GO TO 340 IV(XIRC) = 14 C C. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . C 280 IF (IV(IRC) .NE. 3) GO TO 290 STEP1 = IV(STEP) TEMP1 = IV(STLSTG) C C *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** C L = IV(LMAT) CALL DL7TVM(N, V(TEMP1), V(L), V(STEP1)) CALL DL7VML(N, V(TEMP1), V(L), V(TEMP1)) C C *** COMPUTE GRADIENT *** C 290 IV(NGCALL) = IV(NGCALL) + 1 IV(1) = 2 GO TO 999 C C *** INITIALIZATIONS -- G0 = G - G0, ETC. *** C 300 G01 = IV(G0) CALL DV2AXY(N, V(G01), NEGONE, V(G01), G) STEP1 = IV(STEP) TEMP1 = IV(STLSTG) IF (IV(IRC) .NE. 3) GO TO 320 C C *** SET V(RADFAC) BY GRADIENT TESTS *** C C *** SET TEMP1 = DIAG(D)**-1 * (HESSIAN*STEP + (G(X0)-G(X))) *** C CALL DV2AXY(N, V(TEMP1), NEGONE, V(G01), V(TEMP1)) CALL DV7VMP(N, V(TEMP1), V(TEMP1), D, -1) C C *** DO GRADIENT TESTS *** C IF (DV2NRM(N, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4)) 1 GO TO 310 IF (DD7TPR(N, G, V(STEP1)) 1 .GE. V(GTSTEP) * V(TUNER5)) GO TO 320 310 V(RADFAC) = V(INCFAC) C C *** UPDATE H, LOOP *** C 320 W = IV(NWTSTP) Z = IV(X0) L = IV(LMAT) CALL DW7ZBF(V(L), N, V(STEP1), V(W), V(G01), V(Z)) C C ** USE THE N-VECTORS STARTING AT V(STEP1) AND V(G01) FOR SCRATCH.. CALL DL7UPD(V(TEMP1), V(STEP1), V(L), V(G01), V(L), N, V(W), V(Z)) IV(1) = 2 GO TO 80 C C. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . C C *** BAD PARAMETERS TO ASSESS *** C 330 IV(1) = 64 GO TO 350 C C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** C 340 IV(1) = IV(CNVCOD) IV(CNVCOD) = 0 350 CALL DITSUM(D, G, IV, LIV, LV, N, V, X) C 999 RETURN C C *** LAST LINE OF DRMNG FOLLOWS *** END