SUBROUTINE G7ITB(B, D, G, IV, LIV, LV, P, PS, V, X, Y) C C *** CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR *** C *** REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE) *** C *** HAVING SIMPLE BOUNDS ON THE PARAMETERS BEING ESTIMATED. *** C C *** PARAMETER DECLARATIONS *** C INTEGER LIV, LV, P, PS INTEGER IV(LIV) REAL B(2,P), D(P), G(P), V(LV), X(P), Y(P) C C-------------------------- PARAMETER USAGE -------------------------- C C B.... VECTOR OF LOWER AND UPPER BOUNDS ON X. C D.... SCALE VECTOR. C IV... INTEGER VALUE ARRAY. C LIV.. LENGTH OF IV. MUST BE AT LEAST 80. C LH... LENGTH OF H = P*(P+1)/2. C LV... LENGTH OF V. MUST BE AT LEAST P*(3*P + 19)/2 + 7. C G.... GRADIENT AT X (WHEN IV(1) = 2). C HC... GAUSS-NEWTON HESSIAN AT X (WHEN IV(1) = 2). C P.... NUMBER OF PARAMETERS (COMPONENTS IN X). C PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S. C V.... FLOATING-POINT VALUE ARRAY. C X.... PARAMETER VECTOR. C Y.... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE). C C *** DISCUSSION *** C C G7ITB IS SIMILAR TO G7LIT, EXCEPT FOR THE EXTRA PARAMETER B C -- G7ITB ENFORCES THE BOUNDS B(1,I) .LE. X(I) .LE. B(2,I), C I = 1(1)P. C G7ITB PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF C REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES C IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED C FIRST-ORDER TERM AND A SECOND-ORDER TERM. THE CALLER SUPPLIES C THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED C COMPACTLY BY ROWS), AND G7ITB BUILDS AN APPROXIMATION, S, TO THE C SECOND-ORDER TERM. THE CALLER ALSO PROVIDES THE FUNCTION VALUE, C GRADIENT, AND PART OF THE YIELD VECTOR USED IN UPDATING S. C G7ITB DECIDES DYNAMICALLY WHETHER OR NOT TO USE S WHEN CHOOSING C THE NEXT STEP TO TRY... THE HESSIAN APPROXIMATION USED IS EITHER C HC ALONE (GAUSS-NEWTON MODEL) OR HC + S (AUGMENTED MODEL). C IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT C CONSTANT. THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO C 1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM, OR THE CALLER SETS C IV(INITS) TO 3 OR 4 AND THE FINITE-DIFFERENCE INITIAL S THEN C COMPUTED HAS NONZERO VALUES IN THESE ROWS. C C IF IV(INITS) IS 3 OR 4, THEN THE INITIAL S IS COMPUTED BY C FINITE DIFFERENCES. 3 MEANS USE FUNCTION DIFFERENCES, 4 MEANS C USE GRADIENT DIFFERENCES. FINITE DIFFERENCING IS DONE THE SAME C WAY AS IN COMPUTING A COVARIANCE MATRIX (WITH IV(COVREQ) = -1, -2, C 1, OR 2). C C FOR UPDATING S, G7ITB ASSUMES THAT THE GRADIENT HAS THE FORM C OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE C GRADIENT WITH RESPECT TO X. THE TRUE SECOND-ORDER TERM THEN IS C THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)). IF X = X0 + STEP, C THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF C RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))). THE CALLER MUST SUPPLY C PART OF THIS IN Y, NAMELY THE SUM OVER I OF C RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING G7ITB WITH IV(1) = 2 AND C IV(MODE) = 0 (WHERE MODE = 38). G THEN CONTANS THE OTHER PART, C SO THAT THE DESIRED YIELD VECTOR IS G - Y. IF PS .LT. P, THEN C THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF C GRAD(R(I,X)), STEP, AND Y. C C PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING C ONES TO N2GB (AND NL2SOL), EXCEPT THAT V CAN BE SHORTER C (SINCE THE PART OF V THAT N2GB USES FOR STORING D, J, AND R IS C NOT NEEDED). MOREOVER, COMPARED WITH N2GB (AND NL2SOL), IV(1) C MAY HAVE THE TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE C EXPLAINED BELOW, AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). C THE VALUES IV(D), IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM C N2GB (AND N2FB), ARE NOT REFERENCED BY G7ITB OR THE C SUBROUTINES IT CALLS. C C WHEN G7ITB IS FIRST CALLED, I.E., WHEN G7ITB IS CALLED WITH C IV(1) = 0 OR 12, V(F), G, AND HC NEED NOT BE INITIALIZED. TO C OBTAIN THESE STARTING VALUES, G7ITB RETURNS FIRST WITH IV(1) = 1, C THEN WITH IV(1) = 2, WITH IV(MODE) = -1 IN BOTH CASES. ON C SUBSEQUENT RETURNS WITH IV(1) = 2, IV(MODE) = 0 IMPLIES THAT C Y MUST ALSO BE SUPPLIED. (NOTE THAT Y IS USED FOR SCRATCH -- ITS C INPUT CONTENTS ARE LOST. BY CONTRAST, HC IS NEVER CHANGED.) C ONCE CONVERGENCE HAS BEEN OBTAINED, IV(RDREQ) AND IV(COVREQ) MAY C IMPLY THAT A FINITE-DIFFERENCE HESSIAN SHOULD BE COMPUTED FOR USE C IN COMPUTING A COVARIANCE MATRIX. IN THIS CASE G7ITB WILL MAKE C A NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE. C WHEN IV(MODE) IS POSITIVE, Y SHOULD NOT BE CHANGED. C C IV(1) = 1 MEANS THE CALLER SHOULD SET V(F) (I.E., V(10)) TO F(X), THE C FUNCTION VALUE AT X, AND CALL G7ITB AGAIN, HAVING CHANGED C NONE OF THE OTHER PARAMETERS. AN EXCEPTION OCCURS IF F(X) C CANNOT BE EVALUATED (E.G. IF OVERFLOW WOULD OCCUR), WHICH C MAY HAPPEN BECAUSE OF AN OVERSIZED STEP. IN THIS CASE C THE CALLER SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL C CAUSE G7ITB TO IGNORE V(F) AND TRY A SMALLER STEP. NOTE C THAT THE CURRENT FUNCTION EVALUATION COUNT IS AVAILABLE C IN IV(NFCALL) = IV(6). THIS MAY BE USED TO IDENTIFY C WHICH COPY OF SAVED INFORMATION SHOULD BE USED IN COM- C PUTING G, HC, AND Y THE NEXT TIME G7ITB RETURNS WITH C IV(1) = 2. SEE MLPIT FOR AN EXAMPLE OF THIS. C IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT OF F AT C X. THE CALLER SHOULD ALSO SET HC TO THE GAUSS-NEWTON C HESSIAN AT X. IF IV(MODE) = 0, THEN THE CALLER SHOULD C ALSO COMPUTE THE PART OF THE YIELD VECTOR DESCRIBED ABOVE. C THE CALLER SHOULD THEN CALL G7ITB AGAIN (WITH IV(1) = 2). C THE CALLER MAY ALSO CHANGE D AT THIS TIME, BUT SHOULD NOT C CHANGE X. NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE C VALUE THAT IV(NFCALL) HAD DURING THE RETURN WITH C IV(1) = 1 IN WHICH X HAD THE SAME VALUE AS IT NOW HAS. C IV(NFGCAL) IS EITHER IV(NFCALL) OR IV(NFCALL) - 1. MLPIT C IS AN EXAMPLE WHERE THIS INFORMATION IS USED. IF G OR HC C CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET C IV(NFGCAL) TO 0, IN WHICH CASE G7ITB WILL RETURN WITH C IV(1) = 15. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C (SEE NL2SOL FOR REFERENCES.) C C+++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C LOGICAL HAVQTR, HAVRM INTEGER DUMMY, DIG1, G01, H1, HC1, I, I1, IPI, IPIV0, IPIV1, 1 IPIV2, IPN, J, K, L, LMAT1, LSTGST, P1, P1LEN, PP1, PP1O2, 2 QTR1, RMAT1, RSTRST, STEP1, STPMOD, S1, TD1, TEMP1, TEMP2, 3 TG1, W1, WLM1, X01 REAL E, GI, STTSST, T, T1, XI C C *** CONSTANTS *** C REAL HALF, NEGONE, ONE, ONEP2, ZERO C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C LOGICAL STOPX REAL D7TPR, RLDST, V2NRM EXTERNAL A7SST, D7TPR, F7DHB, G7QSB,I7COPY, I7PNVR, I7SHFT, 1 ITSUM, L7MSB, L7SQR, L7TVM, L7VML, PARCK, Q7RSH, 2 RLDST, S7DMP, S7IPR, S7LUP, S7LVM, STOPX, V2NRM, 3 V2AXY, V7CPY, V7IPR, V7SCP, V7VMP C C A7SST.... ASSESSES CANDIDATE STEP. C D7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. C F7DHB... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR INIT. S MATRIX). C G7QSB... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). C I7COPY.... COPIES ONE INTEGER VECTOR TO ANOTHER. C I7PNVR... INVERTS PERMUTATION ARRAY. C I7SHFT... SHIFTS AN INTEGER VECTOR. C ITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. C L7MSB... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). C L7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L. C L7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C L7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C PARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS. C Q7RSH... SHIFTS A QR FACTORIZATION. C RLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. C S7DMP... MULTIPLIES A SYM. MATRIX FORE AND AFT BY A DIAG. MATRIX. C S7IPR... APPLIES PERMUTATION TO (LOWER TRIANG. OF) SYM. MATRIX. C S7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI- C ANGLE OF A SYMMETRIC MATRIX. C S7LVM... MULTIPLIES COMPACTLY STORED SYM. MATRIX TIMES VECTOR. 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... APPLIES A PERMUTATION TO A VECTOR. C V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C V7VMP... MULTIPLIES (DIVIDES) VECTORS COMPONENTWISE. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, COSMIN, COVMAT, COVREQ, DGNORM, DIG, 1 DSTNRM, F, FDH, FDIF, FUZZ, F0, GTSTEP, H, HC, IERR, 2 INCFAC, INITS, IPIVOT, IRC, IVNEED, KAGQT, KALM, LMAT, 3 LMAX0, LMAXS, MODE, MODEL, MXFCAL, MXITER, NEXTIV, NEXTV, 4 NFCALL, NFGCAL, NFCOV, NGCOV, NGCALL, NITER, NVSAVE, P0, 5 PC, PERM, PHMXFC, PREDUC, QTR, RADFAC, RADINC, RADIUS, 6 RAD0, RDREQ, REGD, RELDX, RESTOR, RMAT, S, SIZE, STEP, 7 STGLIM, STPPAR, SUSED, SWITCH, TOOBIG, TUNER4, TUNER5, 8 VNEED, VSAVE, W, WSCALE, XIRC, X0 C C *** IV SUBSCRIPT VALUES *** C C *** (NOTE THAT P0 AND PC ARE STORED IN IV(G0) AND IV(STLSTG) RESP.) C C/6 C DATA CNVCOD/55/, COVMAT/26/, COVREQ/15/, DIG/37/, FDH/74/, H/56/, C 1 HC/71/, IERR/75/, INITS/25/, IPIVOT/76/, IRC/29/, IVNEED/3/, C 2 KAGQT/33/, KALM/34/, LMAT/42/, MODE/35/, MODEL/5/, C 3 MXFCAL/17/, MXITER/18/, NEXTIV/46/, NEXTV/47/, NFCALL/6/, C 4 NFGCAL/7/, NFCOV/52/, NGCOV/53/, NGCALL/30/, NITER/31/, C 5 P0/48/, PC/41/, PERM/58/, QTR/77/, RADINC/8/, RDREQ/57/, C 6 REGD/67/, RESTOR/9/, RMAT/78/, S/62/, STEP/40/, STGLIM/11/, C 7 SUSED/64/, SWITCH/12/, TOOBIG/2/, VNEED/4/, VSAVE/60/, W/65/, C 8 XIRC/13/, X0/43/ C/7 PARAMETER (CNVCOD=55, COVMAT=26, COVREQ=15, DIG=37, FDH=74, H=56, 1 HC=71, IERR=75, INITS=25, IPIVOT=76, IRC=29, IVNEED=3, 2 KAGQT=33, KALM=34, LMAT=42, MODE=35, MODEL=5, 3 MXFCAL=17, MXITER=18, NEXTIV=46, NEXTV=47, NFCALL=6, 4 NFGCAL=7, NFCOV=52, NGCOV=53, NGCALL=30, NITER=31, 5 P0=48, PC=41, PERM=58, QTR=77, RADINC=8, RDREQ=57, 6 REGD=67, RESTOR=9, RMAT=78, S=62, STEP=40, STGLIM=11, 7 SUSED=64, SWITCH=12, TOOBIG=2, VNEED=4, VSAVE=60, W=65, 8 XIRC=13, X0=43) C/ C C *** V SUBSCRIPT VALUES *** C C/6 C DATA COSMIN/47/, DGNORM/1/, DSTNRM/2/, F/10/, FDIF/11/, FUZZ/45/, C 1 F0/13/, GTSTEP/4/, INCFAC/23/, LMAX0/35/, LMAXS/36/, C 2 NVSAVE/9/, PHMXFC/21/, PREDUC/7/, RADFAC/16/, RADIUS/8/, C 3 RAD0/9/, RELDX/17/, SIZE/55/, STPPAR/5/, TUNER4/29/, C 4 TUNER5/30/, WSCALE/56/ C/7 PARAMETER (COSMIN=47, DGNORM=1, DSTNRM=2, F=10, FDIF=11, FUZZ=45, 1 F0=13, GTSTEP=4, INCFAC=23, LMAX0=35, LMAXS=36, 2 NVSAVE=9, PHMXFC=21, PREDUC=7, RADFAC=16, RADIUS=8, 3 RAD0=9, RELDX=17, SIZE=55, STPPAR=5, TUNER4=29, 4 TUNER5=30, WSCALE=56) C/ C C C/6 C DATA HALF/0.5E+0/, NEGONE/-1.E+0/, ONE/1.E+0/, ONEP2/1.2E+0/, C 1 ZERO/0.E+0/ C/7 PARAMETER (HALF=0.5E+0, NEGONE=-1.E+0, ONE=1.E+0, ONEP2=1.2E+0, 1 ZERO=0.E+0) C/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C I = IV(1) IF (I .EQ. 1) GO TO 50 IF (I .EQ. 2) GO TO 60 C IF (I .LT. 12) GO TO 10 IF (I .GT. 13) GO TO 10 IV(VNEED) = IV(VNEED) + P*(3*P + 25)/2 + 7 IV(IVNEED) = IV(IVNEED) + 4*P 10 CALL PARCK(1, D, IV, LIV, LV, P, V) I = IV(1) - 2 IF (I .GT. 12) GO TO 999 GO TO (360, 360, 360, 360, 360, 360, 240, 190, 240, 20, 20, 30), I C C *** STORAGE ALLOCATION *** C 20 PP1O2 = P * (P + 1) / 2 IV(S) = IV(LMAT) + PP1O2 IV(X0) = IV(S) + PP1O2 IV(STEP) = IV(X0) + 2*P IV(DIG) = IV(STEP) + 3*P IV(W) = IV(DIG) + 2*P IV(H) = IV(W) + 4*P + 7 IV(NEXTV) = IV(H) + PP1O2 IV(IPIVOT) = IV(PERM) + 3*P IV(NEXTIV) = IV(IPIVOT) + P 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(STGLIM) = 2 IV(TOOBIG) = 0 IV(CNVCOD) = 0 IV(COVMAT) = 0 IV(NFCOV) = 0 IV(NGCOV) = 0 IV(RADINC) = 0 IV(PC) = P V(RAD0) = ZERO V(STPPAR) = ZERO V(RADIUS) = V(LMAX0) / (ONE + V(PHMXFC)) C C *** CHECK CONSISTENCY OF B AND INITIALIZE IP ARRAY *** C IPI = IV(IPIVOT) DO 40 I = 1, P IV(IPI) = I IPI = IPI + 1 IF (B(1,I) .GT. B(2,I)) GO TO 680 40 CONTINUE C C *** SET INITIAL MODEL AND S MATRIX *** C IV(MODEL) = 1 IV(1) = 1 IF (IV(S) .LT. 0) GO TO 710 IF (IV(INITS) .GT. 1) IV(MODEL) = 2 S1 = IV(S) IF (IV(INITS) .EQ. 0 .OR. IV(INITS) .GT. 2) 1 CALL V7SCP(P*(P+1)/2, V(S1), ZERO) GO TO 710 C C *** NEW FUNCTION VALUE *** C 50 IF (IV(MODE) .EQ. 0) GO TO 360 IF (IV(MODE) .GT. 0) GO TO 590 C IF (IV(TOOBIG) .EQ. 0) GO TO 690 IV(1) = 63 GO TO 999 C C *** MAKE SURE GRADIENT COULD BE COMPUTED *** C 60 IF (IV(TOOBIG) .EQ. 0) GO TO 70 IV(1) = 65 GO TO 999 C C *** NEW GRADIENT *** C 70 IV(KALM) = -1 IV(KAGQT) = -1 IV(FDH) = 0 IF (IV(MODE) .GT. 0) GO TO 590 IF (IV(HC) .LE. 0 .AND. IV(RMAT) .LE. 0) GO TO 670 C C *** CHOOSE INITIAL PERMUTATION *** C IPI = IV(IPIVOT) IPN = IPI + P - 1 IPIV2 = IV(PERM) - 1 K = IV(PC) P1 = P PP1 = P + 1 RMAT1 = IV(RMAT) HAVRM = RMAT1 .GT. 0 QTR1 = IV(QTR) HAVQTR = QTR1 .GT. 0 C *** MAKE SURE V(QTR1) IS LEGAL (EVEN WHEN NOT REFERENCED) *** W1 = IV(W) IF (.NOT. HAVQTR) QTR1 = W1 + P C DO 100 I = 1, P I1 = IV(IPN) IPN = IPN - 1 IF (B(1,I1) .GE. B(2,I1)) GO TO 80 XI = X(I1) GI = G(I1) IF (XI .LE. B(1,I1) .AND. GI .GT. ZERO) GO TO 80 IF (XI .GE. B(2,I1) .AND. GI .LT. ZERO) GO TO 80 C *** DISALLOW CONVERGENCE IF X(I1) HAS JUST BEEN FREED *** J = IPIV2 + I1 IF (IV(J) .GT. K) IV(CNVCOD) = 0 GO TO 100 80 IF (I1 .GE. P1) GO TO 90 I1 = PP1 - I CALL I7SHFT(P1, I1, IV(IPI)) IF (HAVRM) 1 CALL Q7RSH(I1, P1, HAVQTR, V(QTR1), V(RMAT1), V(W1)) 90 P1 = P1 - 1 100 CONTINUE IV(PC) = P1 C C *** COMPUTE V(DGNORM) (AN OUTPUT VALUE IF WE STOP NOW) *** C V(DGNORM) = ZERO IF (P1 .LE. 0) GO TO 110 DIG1 = IV(DIG) CALL V7VMP(P, V(DIG1), G, D, -1) CALL V7IPR(P, IV(IPI), V(DIG1)) V(DGNORM) = V2NRM(P1, V(DIG1)) 110 IF (IV(CNVCOD) .NE. 0) GO TO 580 IF (IV(MODE) .EQ. 0) GO TO 510 IV(MODE) = 0 V(F0) = V(F) IF (IV(INITS) .LE. 2) GO TO 170 C C *** ARRANGE FOR FINITE-DIFFERENCE INITIAL S *** C IV(XIRC) = IV(COVREQ) IV(COVREQ) = -1 IF (IV(INITS) .GT. 3) IV(COVREQ) = 1 IV(CNVCOD) = 70 GO TO 600 C C *** COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S *** C 120 H1 = IV(FDH) IF (H1 .LE. 0) GO TO 660 IV(CNVCOD) = 0 IV(MODE) = 0 IV(NFCOV) = 0 IV(NGCOV) = 0 IV(COVREQ) = IV(XIRC) S1 = IV(S) PP1O2 = PS * (PS + 1) / 2 HC1 = IV(HC) IF (HC1 .LE. 0) GO TO 130 CALL V2AXY(PP1O2, V(S1), NEGONE, V(HC1), V(H1)) GO TO 140 130 RMAT1 = IV(RMAT) LMAT1 = IV(LMAT) CALL L7SQR(P, V(LMAT1), V(RMAT1)) IPI = IV(IPIVOT) IPIV1 = IV(PERM) + P CALL I7PNVR(P, IV(IPIV1), IV(IPI)) CALL S7IPR(P, IV(IPIV1), V(LMAT1)) CALL V2AXY(PP1O2, V(S1), NEGONE, V(LMAT1), V(H1)) C C *** ZERO PORTION OF S CORRESPONDING TO FIXED X COMPONENTS *** C 140 DO 160 I = 1, P IF (B(1,I) .LT. B(2,I)) GO TO 160 K = S1 + I*(I-1)/2 CALL V7SCP(I, V(K), ZERO) IF (I .GE. P) GO TO 170 K = K + 2*I - 1 I1 = I + 1 DO 150 J = I1, P V(K) = ZERO K = K + J 150 CONTINUE 160 CONTINUE C 170 IV(1) = 2 C C C----------------------------- MAIN LOOP ----------------------------- C C C *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** C 180 CALL ITSUM(D, G, IV, LIV, LV, P, V, X) 190 K = IV(NITER) IF (K .LT. IV(MXITER)) GO TO 200 IV(1) = 10 GO TO 999 200 IV(NITER) = K + 1 C C *** UPDATE RADIUS *** C IF (K .EQ. 0) GO TO 220 STEP1 = IV(STEP) DO 210 I = 1, P V(STEP1) = D(I) * V(STEP1) STEP1 = STEP1 + 1 210 CONTINUE STEP1 = IV(STEP) T = V(RADFAC) * V2NRM(P, V(STEP1)) IF (V(RADFAC) .LT. ONE .OR. T .GT. V(RADIUS)) V(RADIUS) = T C C *** INITIALIZE FOR START OF NEXT ITERATION *** C 220 X01 = IV(X0) V(F0) = V(F) IV(IRC) = 4 IV(H) = -IABS(IV(H)) IV(SUSED) = IV(MODEL) C C *** COPY X TO X0 *** C CALL V7CPY(P, V(X01), X) C C *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** C 230 IF (.NOT. STOPX(DUMMY)) GO TO 250 IV(1) = 11 GO TO 260 C C *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. C 240 IF (V(F) .GE. V(F0)) GO TO 250 V(RADFAC) = ONE K = IV(NITER) GO TO 200 C 250 IF (IV(NFCALL) .LT. IV(MXFCAL) + IV(NFCOV)) GO TO 270 IV(1) = 9 260 IF (V(F) .GE. V(F0)) GO TO 999 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 500 C C. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . C 270 STEP1 = IV(STEP) TG1 = IV(DIG) TD1 = TG1 + P X01 = IV(X0) W1 = IV(W) H1 = IV(H) P1 = IV(PC) IPI = IV(PERM) IPIV1 = IPI + P IPIV2 = IPIV1 + P IPIV0 = IV(IPIVOT) IF (IV(MODEL) .EQ. 2) GO TO 280 C C *** COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE... C RMAT1 = IV(RMAT) IF (RMAT1 .LE. 0) GO TO 280 QTR1 = IV(QTR) IF (QTR1 .LE. 0) GO TO 280 LMAT1 = IV(LMAT) WLM1 = W1 + P CALL L7MSB(B, D, G, IV(IERR), IV(IPIV0), IV(IPIV1), 1 IV(IPIV2), IV(KALM), V(LMAT1), LV, P, IV(P0), 2 IV(PC), V(QTR1), V(RMAT1), V(STEP1), V(TD1), 3 V(TG1), V, V(W1), V(WLM1), X, V(X01)) C *** H IS STORED IN THE END OF W AND HAS JUST BEEN OVERWRITTEN, C *** SO WE MARK IT INVALID... IV(H) = -IABS(H1) C *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO C *** MARK INVALID THE INFORMATION G7QTS MAY HAVE STORED IN V... IV(KAGQT) = -1 GO TO 330 C 280 IF (H1 .GT. 0) GO TO 320 C C *** SET H TO D**-1 * (HC + T1*S) * D**-1. *** C P1LEN = P1*(P1+1)/2 H1 = -H1 IV(H) = H1 IV(FDH) = 0 IF (P1 .LE. 0) GO TO 320 C *** MAKE TEMPORARY PERMUTATION ARRAY *** CALL I7COPY(P, IV(IPI), IV(IPIV0)) J = IV(HC) IF (J .GT. 0) GO TO 290 J = H1 RMAT1 = IV(RMAT) CALL L7SQR(P1, V(H1), V(RMAT1)) GO TO 300 290 CALL V7CPY(P*(P+1)/2, V(H1), V(J)) CALL S7IPR(P, IV(IPI), V(H1)) 300 IF (IV(MODEL) .EQ. 1) GO TO 310 LMAT1 = IV(LMAT) S1 = IV(S) CALL V7CPY(P*(P+1)/2, V(LMAT1), V(S1)) CALL S7IPR(P, IV(IPI), V(LMAT1)) CALL V2AXY(P1LEN, V(H1), ONE, V(LMAT1), V(H1)) 310 CALL V7CPY(P, V(TD1), D) CALL V7IPR(P, IV(IPI), V(TD1)) CALL S7DMP(P1, V(H1), V(H1), V(TD1), -1) IV(KAGQT) = -1 C C *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP *** C 320 LMAT1 = IV(LMAT) CALL G7QSB(B, D, V(H1), G, IV(IPI), IV(IPIV1), IV(IPIV2), 1 IV(KAGQT), V(LMAT1), LV, P, IV(P0), P1, V(STEP1), 2 V(TD1), V(TG1), V, V(W1), X, V(X01)) IF (IV(KALM) .GT. 0) IV(KALM) = 0 C 330 IF (IV(IRC) .NE. 6) GO TO 340 IF (IV(RESTOR) .NE. 2) GO TO 360 RSTRST = 2 GO TO 370 C C *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** C 340 IV(TOOBIG) = 0 IF (V(DSTNRM) .LE. ZERO) GO TO 360 IF (IV(IRC) .NE. 5) GO TO 350 IF (V(RADFAC) .LE. ONE) GO TO 350 IF (V(PREDUC) .GT. ONEP2 * V(FDIF)) GO TO 350 STEP1 = IV(STEP) X01 = IV(X0) CALL V2AXY(P, V(STEP1), NEGONE, V(X01), X) IF (IV(RESTOR) .NE. 2) GO TO 360 RSTRST = 0 GO TO 370 C C *** COMPUTE F(X0 + STEP) *** C 350 X01 = IV(X0) STEP1 = IV(STEP) CALL V2AXY(P, X, ONE, V(STEP1), V(X01)) IV(NFCALL) = IV(NFCALL) + 1 IV(1) = 1 GO TO 710 C C. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . C 360 RSTRST = 3 370 X01 = IV(X0) V(RELDX) = RLDST(P, D, X, V(X01)) CALL A7SST(IV, LIV, LV, V) STEP1 = IV(STEP) LSTGST = X01 + P I = IV(RESTOR) + 1 GO TO (410, 380, 390, 400), I 380 CALL V7CPY(P, X, V(X01)) GO TO 410 390 CALL V7CPY(P, V(LSTGST), V(STEP1)) GO TO 410 400 CALL V7CPY(P, V(STEP1), V(LSTGST)) CALL V2AXY(P, X, ONE, V(STEP1), V(X01)) V(RELDX) = RLDST(P, D, X, V(X01)) IV(RESTOR) = RSTRST C C *** IF NECESSARY, SWITCH MODELS *** C 410 IF (IV(SWITCH) .EQ. 0) GO TO 420 IV(H) = -IABS(IV(H)) IV(SUSED) = IV(SUSED) + 2 L = IV(VSAVE) CALL V7CPY(NVSAVE, V, V(L)) 420 L = IV(IRC) - 4 STPMOD = IV(MODEL) IF (L .GT. 0) GO TO (440,450,460,460,460,460,460,460,570,510), L C C *** DECIDE WHETHER TO CHANGE MODELS *** C E = V(PREDUC) - V(FDIF) S1 = IV(S) CALL S7LVM(PS, Y, V(S1), V(STEP1)) STTSST = HALF * D7TPR(PS, V(STEP1), Y) IF (IV(MODEL) .EQ. 1) STTSST = -STTSST IF ( ABS(E + STTSST) * V(FUZZ) .GE. ABS(E)) GO TO 430 C C *** SWITCH MODELS *** C IV(MODEL) = 3 - IV(MODEL) IF (-2 .LT. L) GO TO 470 IV(H) = -IABS(IV(H)) IV(SUSED) = IV(SUSED) + 2 L = IV(VSAVE) CALL V7CPY(NVSAVE, V(L), V) GO TO 230 C 430 IF (-3 .LT. L) GO TO 470 C C *** RECOMPUTE STEP WITH DIFFERENT RADIUS *** C 440 V(RADIUS) = V(RADFAC) * V(DSTNRM) GO TO 230 C C *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST C 450 V(RADIUS) = V(LMAXS) GO TO 270 C C *** CONVERGENCE OR FALSE CONVERGENCE *** C 460 IV(CNVCOD) = L IF (V(F) .GE. V(F0)) GO TO 580 IF (IV(XIRC) .EQ. 14) GO TO 580 IV(XIRC) = 14 C C. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . C 470 IV(COVMAT) = 0 IV(REGD) = 0 C C *** SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS *** C IF (IV(IRC) .NE. 3) GO TO 500 STEP1 = IV(STEP) TEMP1 = STEP1 + P TEMP2 = IV(X0) C C *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** C HC1 = IV(HC) IF (HC1 .LE. 0) GO TO 480 CALL S7LVM(P, V(TEMP1), V(HC1), V(STEP1)) GO TO 490 480 RMAT1 = IV(RMAT) IPIV0 = IV(IPIVOT) CALL V7CPY(P, V(TEMP1), V(STEP1)) CALL V7IPR(P, IV(IPIV0), V(TEMP1)) CALL L7TVM(P, V(TEMP1), V(RMAT1), V(TEMP1)) CALL L7VML(P, V(TEMP1), V(RMAT1), V(TEMP1)) IPIV1 = IV(PERM) + P CALL I7PNVR(P, IV(IPIV1), IV(IPIV0)) CALL V7IPR(P, IV(IPIV1), V(TEMP1)) C 490 IF (STPMOD .EQ. 1) GO TO 500 S1 = IV(S) CALL S7LVM(PS, V(TEMP2), V(S1), V(STEP1)) CALL V2AXY(PS, V(TEMP1), ONE, V(TEMP2), V(TEMP1)) C C *** SAVE OLD GRADIENT AND COMPUTE NEW ONE *** C 500 IV(NGCALL) = IV(NGCALL) + 1 G01 = IV(W) CALL V7CPY(P, V(G01), G) GO TO 690 C C *** INITIALIZATIONS -- G0 = G - G0, ETC. *** C 510 G01 = IV(W) CALL V2AXY(P, V(G01), NEGONE, V(G01), G) STEP1 = IV(STEP) TEMP1 = STEP1 + P TEMP2 = IV(X0) IF (IV(IRC) .NE. 3) GO TO 540 C C *** SET V(RADFAC) BY GRADIENT TESTS *** C C *** SET TEMP1 = D**-1 * (HESSIAN * STEP + (G(X0) - G(X))) *** C K = TEMP1 L = G01 DO 520 I = 1, P V(K) = (V(K) - V(L)) / D(I) K = K + 1 L = L + 1 520 CONTINUE C C *** DO GRADIENT TESTS *** C IF ( V2NRM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4)) GO TO 530 IF ( D7TPR(P, G, V(STEP1)) 1 .GE. V(GTSTEP) * V(TUNER5)) GO TO 540 530 V(RADFAC) = V(INCFAC) C C *** COMPUTE Y VECTOR NEEDED FOR UPDATING S *** C 540 CALL V2AXY(PS, Y, NEGONE, Y, G) C C *** DETERMINE SIZING FACTOR V(SIZE) *** C C *** SET TEMP1 = S * STEP *** S1 = IV(S) CALL S7LVM(PS, V(TEMP1), V(S1), V(STEP1)) C T1 = ABS( D7TPR(PS, V(STEP1), V(TEMP1))) T = ABS( D7TPR(PS, V(STEP1), Y)) V(SIZE) = ONE IF (T .LT. T1) V(SIZE) = T / T1 C C *** SET G0 TO WCHMTD CHOICE OF FLETCHER AND AL-BAALI *** C HC1 = IV(HC) IF (HC1 .LE. 0) GO TO 550 CALL S7LVM(PS, V(G01), V(HC1), V(STEP1)) GO TO 560 C 550 RMAT1 = IV(RMAT) IPIV0 = IV(IPIVOT) CALL V7CPY(P, V(G01), V(STEP1)) I = G01 + PS IF (PS .LT. P) CALL V7SCP(P-PS, V(I), ZERO) CALL V7IPR(P, IV(IPIV0), V(G01)) CALL L7TVM(P, V(G01), V(RMAT1), V(G01)) CALL L7VML(P, V(G01), V(RMAT1), V(G01)) IPIV1 = IV(PERM) + P CALL I7PNVR(P, IV(IPIV1), IV(IPIV0)) CALL V7IPR(P, IV(IPIV1), V(G01)) C 560 CALL V2AXY(PS, V(G01), ONE, Y, V(G01)) C C *** UPDATE S *** C CALL S7LUP(V(S1), V(COSMIN), PS, V(SIZE), V(STEP1), V(TEMP1), 1 V(TEMP2), V(G01), V(WSCALE), Y) IV(1) = 2 GO TO 180 C C. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . C C *** BAD PARAMETERS TO ASSESS *** C 570 IV(1) = 64 GO TO 999 C C C *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** C 580 IF (IV(RDREQ) .EQ. 0) GO TO 660 IF (IV(FDH) .NE. 0) GO TO 660 IF (IV(CNVCOD) .GE. 7) GO TO 660 IF (IV(REGD) .GT. 0) GO TO 660 IF (IV(COVMAT) .GT. 0) GO TO 660 IF (IABS(IV(COVREQ)) .GE. 3) GO TO 640 IF (IV(RESTOR) .EQ. 0) IV(RESTOR) = 2 GO TO 600 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE *** C 590 IV(RESTOR) = 0 600 CALL F7DHB(B, D, G, I, IV, LIV, LV, P, V, X) GO TO (610, 620, 630), I 610 IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(1) = 1 GO TO 710 C 620 IV(NGCOV) = IV(NGCOV) + 1 IV(NGCALL) = IV(NGCALL) + 1 IV(NFGCAL) = IV(NFCALL) + IV(NGCOV) GO TO 690 C 630 IF (IV(CNVCOD) .EQ. 70) GO TO 120 GO TO 660 C 640 H1 = IABS(IV(H)) IV(FDH) = H1 IV(H) = -H1 HC1 = IV(HC) IF (HC1 .LE. 0) GO TO 650 CALL V7CPY(P*(P+1)/2, V(H1), V(HC1)) GO TO 660 650 RMAT1 = IV(RMAT) CALL L7SQR(P, V(H1), V(RMAT1)) C 660 IV(MODE) = 0 IV(1) = IV(CNVCOD) IV(CNVCOD) = 0 GO TO 999 C C *** SPECIAL RETURN FOR MISSING HESSIAN INFORMATION -- BOTH C *** IV(HC) .LE. 0 AND IV(RMAT) .LE. 0 C 670 IV(1) = 1400 GO TO 999 C C *** INCONSISTENT B *** C 680 IV(1) = 82 GO TO 999 C C *** SAVE, THEN INITIALIZE IPIVOT ARRAY BEFORE COMPUTING G *** C 690 IV(1) = 2 J = IV(IPIVOT) IPI = IV(PERM) CALL I7PNVR(P, IV(IPI), IV(J)) DO 700 I = 1, P IV(J) = I J = J + 1 700 CONTINUE C C *** PROJECT X INTO FEASIBLE REGION (PRIOR TO COMPUTING F OR G) *** C 710 DO 720 I = 1, P IF (X(I) .LT. B(1,I)) X(I) = B(1,I) IF (X(I) .GT. B(2,I)) X(I) = B(2,I) 720 CONTINUE IV(TOOBIG) = 0 C 999 RETURN C C *** LAST LINE OF G7ITB FOLLOWS *** END .