SUBROUTINE RMNGB(B, D, FX, G, IV, LIV, LV, N, V, X) C C *** CARRY OUT MNGB (SIMPLY BOUNDED MINIMIZATION) ITERATIONS, C *** USING DOUBLE-DOGLEG/BFGS STEPS. C C *** PARAMETER DECLARATIONS *** C INTEGER LIV, LV, N INTEGER IV(LIV) REAL B(2,N), D(N), FX, G(N), V(LV), X(N) C C-------------------------- PARAMETER USAGE -------------------------- C C B.... VECTOR OF LOWER AND UPPER BOUNDS ON X. C D.... SCALE VECTOR. C FX... FUNCTION VALUE. C G.... GRADIENT VECTOR. C IV... INTEGER VALUE ARRAY. C LIV.. LENGTH OF IV (AT LEAST 59) + N. C LV... LENGTH OF V (AT LEAST 71 + N*(N+19)/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 MNGB (WHICH SEE), EXCEPT THAT V CAN BE SHORTER (SINCE C THE PART OF V THAT MNGB USES FOR STORING G IS NOT NEEDED). C MOREOVER, COMPARED WITH MNGB, 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 MNGB (AND SMSNOB), IS NOT REFERENCED BY C RMNGB OR THE SUBROUTINES IT CALLS. C FX AND G NEED NOT HAVE BEEN INITIALIZED WHEN RMNGB 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 RMNGB 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 RMNGB TO IG- C NORE FX AND TRY A SMALLER STEP. THE PARAMETER NF THAT C MNGB 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 RMNGB AGAIN, HAVING CHANGED NONE OF C THE OTHER PARAMETERS EXCEPT POSSIBLY THE SCALE VECTOR D C WHEN IV(DTYPE) = 0. THE PARAMETER NF THAT MNGB PASSES C TO CALCG IS IV(NFGCAL) = IV(7). IF G(X) CANNOT BE C EVALUATED, THEN THE CALLER MAY SET IV(NFGCAL) TO 0, IN C WHICH CASE RMNGB 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 MNG FOR REFERENCES.) C C+++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C INTEGER DG1, DSTEP1, DUMMY, G01, I, I1, IPI, IPN, J, K, L, LSTGST, 1 N1, NP1, NWTST1, RSTRST, STEP1, TEMP0, TEMP1, TD1, TG1, 2 W1, X01, Z REAL GI, T, XI C C *** CONSTANTS *** C REAL NEGONE, ONE, ONEP2, ZERO C C *** NO INTRINSIC FUNCTIONS *** C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C LOGICAL STOPX REAL D7TPR, RLDST, V2NRM EXTERNAL A7SST, D7DGB, IVSET, D7TPR, I7SHFT, ITSUM, L7TVM, 1 L7UPD, L7VML, PARCK, Q7RSH, RLDST, STOPX, V2NRM, 2 V2AXY, V7CPY, V7IPR, V7SCP, V7VMP, W7ZBF C C A7SST.... ASSESSES CANDIDATE STEP. C D7DGB... COMPUTES SIMPLY BOUNDED DOUBLE-DOGLEG (CANDIDATE) STEP. C IVSET.... SUPPLIES DEFAULT IV AND V INPUT COMPONENTS. C D7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. C I7SHFT... CYCLICALLLY SHIFTS AN ARRAY OF INTEGERS. C ITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. C L7TVM... MULTIPLIES TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. C LUPDT.... UPDATES CHOLESKY FACTOR OF HESSIAN APPROXIMATION. C L7VML.... MULTIPLIES LOWER TRIANGLE TIMES VECTOR. C PARCK.... CHECKS VALIDITY OF INPUT IV AND V VALUES. C Q7RSH... CYCLICALLY SHIFTS CHOLESKY FACTOR. C RLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. C V2NRM... RETURNS THE 2-NORM OF A VECTOR. C V2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. C V7CPY.... COPIES ONE VECTOR TO ANOTHER. C V7IPR... CYCLICALLY SHIFTS A FLOATING-POINT ARRAY. C V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C V7VMP... MULTIPLIES VECTOR BY VECTOR RAISED TO POWER (COMPONENTWISE). C W7ZBF... COMPUTES W AND Z FOR L7UPD CORRESPONDING TO BFGS UPDATE. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, DG, DGNORM, DINIT, DSTNRM, F, F0, FDIF, 1 GTSTEP, INCFAC, INITH, IRC, IVNEED, KAGQT, LMAT, 2 LMAX0, LMAXS, MODE, MODEL, MXFCAL, MXITER, NC, NEXTIV, 3 NEXTV, NFCALL, NFGCAL, NGCALL, NITER, NWTSTP, PERM, 4 PREDUC, RADFAC, RADINC, RADIUS, RAD0, RELDX, RESTOR, STEP, 4 STGLIM, STLSTG, TOOBIG, TUNER4, TUNER5, VNEED, XIRC, X0 C C *** IV SUBSCRIPT VALUES *** C C *** (NOTE THAT NC IS STORED IN IV(G0)) *** C C/6 C DATA CNVCOD/55/, DG/37/, INITH/25/, IRC/29/, IVNEED/3/, KAGQT/33/, C 1 MODE/35/, MODEL/5/, MXFCAL/17/, MXITER/18/, NC/48/, C 2 NEXTIV/46/, NEXTV/47/, NFCALL/6/, NFGCAL/7/, NGCALL/30/, C 3 NITER/31/, NWTSTP/34/, PERM/58/, RADINC/8/, RESTOR/9/, C 4 STEP/40/, STGLIM/11/, STLSTG/41/, TOOBIG/2/, XIRC/13/, X0/43/ C/7 PARAMETER (CNVCOD=55, DG=37, INITH=25, IRC=29, IVNEED=3, KAGQT=33, 1 MODE=35, MODEL=5, MXFCAL=17, MXITER=18, NC=48, 2 NEXTIV=46, NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, 3 NITER=31, NWTSTP=34, PERM=58, RADINC=8, RESTOR=9, 4 STEP=40, STGLIM=11, STLSTG=41, TOOBIG=2, XIRC=13, 5 X0=43) C/ C C *** V SUBSCRIPT VALUES *** C C/6 C DATA DGNORM/1/, DINIT/38/, DSTNRM/2/, F/10/, F0/13/, FDIF/11/, C 1 GTSTEP/4/, INCFAC/23/, LMAT/42/, LMAX0/35/, LMAXS/36/, C 2 PREDUC/7/, RADFAC/16/, RADIUS/8/, RAD0/9/, RELDX/17/, C 3 TUNER4/29/, TUNER5/30/, VNEED/4/ C/7 PARAMETER (DGNORM=1, DINIT=38, DSTNRM=2, F=10, F0=13, FDIF=11, 1 GTSTEP=4, INCFAC=23, LMAT=42, LMAX0=35, LMAXS=36, 2 PREDUC=7, RADFAC=16, RADIUS=8, RAD0=9, RELDX=17, 3 TUNER4=29, TUNER5=30, VNEED=4) C/ C C/6 C DATA NEGONE/-1.E+0/, ONE/1.E+0/, ONEP2/1.2E+0/, ZERO/0.E+0/ C/7 PARAMETER (NEGONE=-1.E+0, ONE=1.E+0, ONEP2=1.2E+0, ZERO=0.E+0) C/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C I = IV(1) IF (I .EQ. 1) GO TO 70 IF (I .EQ. 2) GO TO 80 C C *** CHECK VALIDITY OF IV AND V INPUT VALUES *** C IF (IV(1) .EQ. 0) CALL IVSET(2, IV, LIV, LV, V) IF (IV(1) .LT. 12) GO TO 10 IF (IV(1) .GT. 13) GO TO 10 IV(VNEED) = IV(VNEED) + N*(N+19)/2 IV(IVNEED) = IV(IVNEED) + N 10 CALL PARCK(2, D, IV, LIV, LV, N, V) I = IV(1) - 2 IF (I .GT. 12) GO TO 999 GO TO (250, 250, 250, 250, 250, 250, 190, 150, 190, 20, 20, 30), I C C *** STORAGE ALLOCATION *** C 20 L = IV(LMAT) IV(X0) = L + N*(N+1)/2 IV(STEP) = IV(X0) + 2*N IV(STLSTG) = IV(STEP) + 2*N IV(NWTSTP) = IV(STLSTG) + N IV(DG) = IV(NWTSTP) + 2*N IV(NEXTV) = IV(DG) + 2*N IV(NEXTIV) = IV(PERM) + N IF (IV(1) .NE. 13) GO TO 30 IV(1) = 14 GO TO 999 C C *** INITIALIZATION *** C 30 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 IV(NC) = N V(RAD0) = ZERO C C *** CHECK CONSISTENCY OF B AND INITIALIZE IP ARRAY *** C IPI = IV(PERM) DO 40 I = 1, N IV(IPI) = I IPI = IPI + 1 IF (B(1,I) .GT. B(2,I)) GO TO 410 40 CONTINUE C IF (V(DINIT) .GE. ZERO) CALL V7SCP(N, D, V(DINIT)) IF (IV(INITH) .NE. 1) GO TO 60 C C *** SET THE INITIAL HESSIAN APPROXIMATION TO DIAG(D)**-2 *** C L = IV(LMAT) CALL V7SCP(N*(N+1)/2, V(L), ZERO) K = L - 1 DO 50 I = 1, N K = K + I T = D(I) IF (T .LE. ZERO) T = ONE V(K) = T 50 CONTINUE C C *** GET INITIAL FUNCTION VALUE *** C 60 IV(1) = 1 GO TO 440 C 70 V(F) = FX IF (IV(MODE) .GE. 0) GO TO 250 V(F0) = FX IV(1) = 2 IF (IV(TOOBIG) .EQ. 0) GO TO 999 IV(1) = 63 GO TO 430 C C *** MAKE SURE GRADIENT COULD BE COMPUTED *** C 80 IF (IV(TOOBIG) .EQ. 0) GO TO 90 IV(1) = 65 GO TO 430 C C *** CHOOSE INITIAL PERMUTATION *** C 90 IPI = IV(PERM) IPN = IPI + N N1 = N NP1 = N + 1 L = IV(LMAT) W1 = IV(NWTSTP) + N K = N - IV(NC) DO 120 I = 1, N IPN = IPN - 1 J = IV(IPN) IF (B(1,J) .GE. B(2,J)) GO TO 100 XI = X(J) GI = G(J) IF (XI .LE. B(1,J) .AND. GI .GT. ZERO) GO TO 100 IF (XI .GE. B(2,J) .AND. GI .LT. ZERO) GO TO 100 C *** DISALLOW CONVERGENCE IF X(J) HAS JUST BEEN FREED *** IF (I .LE. K) IV(CNVCOD) = 0 GO TO 120 100 I1 = NP1 - I IF (I1 .GE. N1) GO TO 110 CALL I7SHFT(N1, I1, IV(IPI)) CALL Q7RSH(I1, N1, .FALSE., G, V(L), V(W1)) 110 N1 = N1 - 1 120 CONTINUE C IV(NC) = N1 V(DGNORM) = ZERO IF (N1 .LE. 0) GO TO 130 DG1 = IV(DG) CALL V7VMP(N, V(DG1), G, D, -1) CALL V7IPR(N, IV(IPI), V(DG1)) V(DGNORM) = V2NRM(N1, V(DG1)) 130 IF (IV(CNVCOD) .NE. 0) GO TO 420 IF (IV(MODE) .EQ. 0) GO TO 370 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 140 CALL ITSUM(D, G, IV, LIV, LV, N, V, X) 150 K = IV(NITER) IF (K .LT. IV(MXITER)) GO TO 160 IV(1) = 10 GO TO 430 C C *** UPDATE RADIUS *** C 160 IV(NITER) = K + 1 IF (K .EQ. 0) GO TO 170 T = V(RADFAC) * V(DSTNRM) IF (V(RADFAC) .LT. ONE .OR. T .GT. V(RADIUS)) V(RADIUS) = T C C *** INITIALIZE FOR START OF NEXT ITERATION *** C 170 X01 = IV(X0) V(F0) = V(F) IV(IRC) = 4 IV(KAGQT) = -1 C C *** COPY X TO X0 *** C CALL V7CPY(N, V(X01), X) C C *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** C 180 IF (.NOT. STOPX(DUMMY)) GO TO 200 IV(1) = 11 GO TO 210 C C *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. C 190 IF (V(F) .GE. V(F0)) GO TO 200 V(RADFAC) = ONE K = IV(NITER) GO TO 160 C 200 IF (IV(NFCALL) .LT. IV(MXFCAL)) GO TO 220 IV(1) = 9 210 IF (V(F) .GE. V(F0)) GO TO 430 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 360 C C. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . C 220 STEP1 = IV(STEP) DG1 = IV(DG) NWTST1 = IV(NWTSTP) W1 = NWTST1 + N DSTEP1 = STEP1 + N IPI = IV(PERM) L = IV(LMAT) TG1 = DG1 + N X01 = IV(X0) TD1 = X01 + N CALL D7DGB(B, D, V(DG1), V(DSTEP1), G, IV(IPI), IV(KAGQT), 1 V(L), LV, N, IV(NC), V(NWTST1), V(STEP1), V(TD1), 2 V(TG1), V, V(W1), V(X01)) IF (IV(IRC) .NE. 6) GO TO 230 IF (IV(RESTOR) .NE. 2) GO TO 250 RSTRST = 2 GO TO 260 C C *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** C 230 IV(TOOBIG) = 0 IF (V(DSTNRM) .LE. ZERO) GO TO 250 IF (IV(IRC) .NE. 5) GO TO 240 IF (V(RADFAC) .LE. ONE) GO TO 240 IF (V(PREDUC) .GT. ONEP2 * V(FDIF)) GO TO 240 IF (IV(RESTOR) .NE. 2) GO TO 250 RSTRST = 0 GO TO 260 C C *** COMPUTE F(X0 + STEP) *** C 240 CALL V2AXY(N, X, ONE, V(STEP1), V(X01)) IV(NFCALL) = IV(NFCALL) + 1 IV(1) = 1 GO TO 440 C C. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . C 250 RSTRST = 3 260 X01 = IV(X0) V(RELDX) = RLDST(N, D, X, V(X01)) CALL A7SST(IV, LIV, LV, V) STEP1 = IV(STEP) LSTGST = IV(STLSTG) I = IV(RESTOR) + 1 GO TO (300, 270, 280, 290), I 270 CALL V7CPY(N, X, V(X01)) GO TO 300 280 CALL V7CPY(N, V(LSTGST), X) GO TO 300 290 CALL V7CPY(N, X, V(LSTGST)) CALL V2AXY(N, V(STEP1), NEGONE, V(X01), X) V(RELDX) = RLDST(N, D, X, V(X01)) IV(RESTOR) = RSTRST C 300 K = IV(IRC) GO TO (310,340,340,340,310,320,330,330,330,330,330,330,400,370), K C C *** RECOMPUTE STEP WITH CHANGED RADIUS *** C 310 V(RADIUS) = V(RADFAC) * V(DSTNRM) GO TO 180 C C *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST. C 320 V(RADIUS) = V(LMAXS) GO TO 220 C C *** CONVERGENCE OR FALSE CONVERGENCE *** C 330 IV(CNVCOD) = K - 4 IF (V(F) .GE. V(F0)) GO TO 420 IF (IV(XIRC) .EQ. 14) GO TO 420 IV(XIRC) = 14 C C. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . C 340 X01 = IV(X0) STEP1 = IV(STEP) CALL V2AXY(N, V(STEP1), NEGONE, V(X01), X) IF (IV(IRC) .NE. 3) GO TO 360 C C *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** C C *** USE X0 AS TEMPORARY... C IPI = IV(PERM) CALL V7CPY(N, V(X01), V(STEP1)) CALL V7IPR(N, IV(IPI), V(X01)) L = IV(LMAT) CALL L7TVM(N, V(X01), V(L), V(X01)) CALL L7VML(N, V(X01), V(L), V(X01)) C C *** UNPERMUTE X0 INTO TEMP1 *** C TEMP1 = IV(STLSTG) TEMP0 = TEMP1 - 1 DO 350 I = 1, N J = IV(IPI) IPI = IPI + 1 K = TEMP0 + J V(K) = V(X01) X01 = X01 + 1 350 CONTINUE C C *** SAVE OLD GRADIENT, COMPUTE NEW ONE *** C 360 G01 = IV(NWTSTP) + N CALL V7CPY(N, V(G01), G) IV(NGCALL) = IV(NGCALL) + 1 IV(TOOBIG) = 0 IV(1) = 2 GO TO 999 C C *** INITIALIZATIONS -- G0 = G - G0, ETC. *** C 370 G01 = IV(NWTSTP) + N CALL V2AXY(N, V(G01), NEGONE, V(G01), G) STEP1 = IV(STEP) TEMP1 = IV(STLSTG) IF (IV(IRC) .NE. 3) GO TO 390 C C *** SET V(RADFAC) BY GRADIENT TESTS *** C C *** SET TEMP1 = DIAG(D)**-1 * (HESSIAN*STEP + (G(X0)-G(X))) *** C CALL V2AXY(N, V(TEMP1), NEGONE, V(G01), V(TEMP1)) CALL V7VMP(N, V(TEMP1), V(TEMP1), D, -1) C C *** DO GRADIENT TESTS *** C IF ( V2NRM(N, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4)) 1 GO TO 380 IF ( D7TPR(N, G, V(STEP1)) 1 .GE. V(GTSTEP) * V(TUNER5)) GO TO 390 380 V(RADFAC) = V(INCFAC) C C *** UPDATE H, LOOP *** C 390 W1 = IV(NWTSTP) Z = IV(X0) L = IV(LMAT) IPI = IV(PERM) CALL V7IPR(N, IV(IPI), V(STEP1)) CALL V7IPR(N, IV(IPI), V(G01)) CALL W7ZBF(V(L), N, V(STEP1), V(W1), V(G01), V(Z)) C C ** USE THE N-VECTORS STARTING AT V(STEP1) AND V(G01) FOR SCRATCH.. CALL L7UPD(V(TEMP1), V(STEP1), V(L), V(G01), V(L), N, V(W1), 1 V(Z)) IV(1) = 2 GO TO 140 C C. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . C C *** BAD PARAMETERS TO ASSESS *** C 400 IV(1) = 64 GO TO 430 C C *** INCONSISTENT B *** C 410 IV(1) = 82 GO TO 430 C C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** C 420 IV(1) = IV(CNVCOD) IV(CNVCOD) = 0 430 CALL ITSUM(D, G, IV, LIV, LV, N, V, X) GO TO 999 C C *** PROJECT X INTO FEASIBLE REGION (PRIOR TO COMPUTING F OR G) *** C 440 DO 450 I = 1, N IF (X(I) .LT. B(1,I)) X(I) = B(1,I) IF (X(I) .GT. B(2,I)) X(I) = B(2,I) 450 CONTINUE C 999 RETURN C C *** LAST CARD OF RMNGB FOLLOWS *** END .