SUBROUTINE GLG(N, P, PS, X, RHO, RHOI, RHOR, IV, LIV, LV, V, 1 CALCRJ, UI, UR, UF) C C *** GENERALIZED LINEAR REGRESSION A LA NL2SOL *** C C *** PARAMETERS *** C INTEGER N, P, PS, LIV, LV INTEGER IV(LIV), RHOI(*), UI(*) REAL X(*), RHOR(*), V(LV), UR(*) EXTERNAL CALCRJ, RHO, UF C C *** PARAMETER USAGE *** C C N....... TOTAL NUMBER OF RESIDUALS. C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS...... NUMBER OF NON-NUISANCE PARAMETERS (THOSE INVOLVED IN S). C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C RHO..... SUBROUTINE FOR COMPUTING LOSS FUNCTIONS AND THEIR DERIVS. C SEE RGLG FOR DETAILS ABOUT RHO. C RHOI.... PASSED WITHOUT CHANGE TO RHO. C RHOR.... PASSED WITHOUT CHANGE TO RHO. C IV...... INTEGER VALUES ARRAY. C LIV..... LENGTH OF IV, AT LEAST 90 + P. C LV...... LENGTH OF V, AT LEAST C 105 + P*(2*P + 17) + 4*PS C + N*(P + 3 + (P-PS+1)*(P-PS+2)/2). C V....... FLOATING-POINT VALUES ARRAY. C CALCRJ.. SUBROUTINE FOR COMPUTING RESIDUAL VECTOR AND JACOBIAN MATRIX. C UI...... PASSED UNCHANGED TO CALCRJ. C UR...... PASSED UNCHANGED TO CALCRJ. C UF...... PASSED UNCHANGED TO CALCRJ. C C *** CALCRJ CALLING SEQUENCE... C C CALL CALCRJ(N, PS, X, NF, NEED, R, RP, UI, UR, UF) C C PARAMETERS N, PS, X, UI, UR, AND UF ARE AS ABOVE. C R AND RP ARE FLOATING-POINT ARRAYS DIMENSIONED R(N) AND RP(PS,N). C NEED IS AN INTEGER ARRAY OF LENGTH 2... C NEED(1) = 1 MEANS CALCRJ SHOULD COMPUTE THE RESIDUAL VECTOR R, C AND NEED(2) IS THE VALUE NF HAD AT THE LAST X WHERE C CALCRJ MIGHT BE CALLED WITH NEED(1) = 2. C NEED(1) = 2 MEANS CALCRJ SHOULD COMPUTE THE JACOBIAN MATRIX RP, C WHERE RP(J,I) = DERIVATIVE OF R(I) WITH RESPECT TO X(J). C (CALCRJ SHOULD NOT CHANGE NEED AND SHOULD CHANGE AT MOST ONE OF R C AND RP. IF R OR RP, AS APPROPRIATE, CANNOT BE COMPUTED, THEN CALCRJ C SHOULD SET NF TO 0. OTHERWISE IT SHOULD NOT CHANGE NF.) C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C *** EXTERNAL SUBROUTINES *** C EXTERNAL IVSET, RGLG C C IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C RGLG ... CARRIES OUT OPTIMIZATION ITERATIONS. C C C *** LOCAL VARIABLES *** C INTEGER D1, DR1, I, IV1, NEED1(2), NEED2(2), NF, R1, RD1 C C *** IV COMPONENTS *** C INTEGER D, J, NEXTV, NFCALL, NFGCAL, R, REGD, REGD0, TOOBIG, VNEED PARAMETER (D=27, J=70, NEXTV=47, NFCALL=6, NFGCAL=7, R=61, 1 REGD=67, REGD0=82, TOOBIG=2, VNEED=4) SAVE NEED1, NEED2 DATA NEED1(1)/1/, NEED1(2)/0/, NEED2(1)/2/, NEED2(2)/0/ C C--------------------------------- BODY ------------------------------ C IF (IV(1) .EQ. 0) CALL IVSET(1, IV, LIV, LV, V) IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 10 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 IF (IV1 .EQ. 12) IV(1) = 13 I = (P-PS+2)*(P-PS+1)/2 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+1+I) CALL RGLG(X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, RHO, RHOI, 1 RHOR, V, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(D) = IV(NEXTV) IV(R) = IV(D) + P IV(REGD0) = IV(R) + (P - PS + 1)*N IV(J) = IV(REGD0) + ((P-PS+2)*(P-PS+1)/2)*N IV(NEXTV) = IV(J) + N*PS IF (IV1 .EQ. 13) GO TO 999 C 10 D1 = IV(D) DR1 = IV(J) R1 = IV(R) RD1 = IV(REGD0) C 20 CALL RGLG(V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS, V(R1), 1 V(RD1), RHO, RHOI, RHOR, V, X) IF (IV(1)-2) 30, 50, 60 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 30 NF = IV(NFCALL) NEED1(2) = IV(NFGCAL) CALL CALCRJ(N, PS, X, NF, NEED1, V(R1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 40 IV(TOOBIG) = 1 GO TO 20 40 IF (IV(1) .GT. 0) GO TO 20 C C *** COMPUTE DR = GRADIENT OF R COMPONENTS *** C 50 CALL CALCRJ(N, PS, X, IV(NFGCAL), NEED2, V(R1), V(DR1), UI, UR,UF) IF (IV(NFGCAL) .EQ. 0) IV(TOOBIG) = 1 GO TO 20 C C *** INDICATE WHETHER THE REGRESSION DIAGNOSTIC ARRAY WAS COMPUTED C *** AND PRINT IT IF SO REQUESTED... C 60 IF (IV(REGD) .GT. 0) IV(REGD) = RD1 C 999 RETURN C C *** LAST LINE OF GLG FOLLOWS *** END SUBROUTINE GLF(N, P, PS, X, RHO, RHOI, RHOR, IV, LIV, LV, V, 1 CALCRJ, UI, UR, UF) C C *** GENERALIZED LINEAR REGRESSION, FINITE-DIFFERENCE JACOBIAN *** C C *** PARAMETERS *** C INTEGER N, P, PS, LIV, LV INTEGER IV(LIV), RHOI(*), UI(*) REAL X(*), V(LV), RHOR(*), UR(*) EXTERNAL CALCRJ, RHO, UF C C *** PARAMETER USAGE *** C C N....... TOTAL NUMBER OF RESIDUALS. C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS...... NUMBER OF NON-NUISANCE PARAMETERS (THOSE INVOLVED IN S). C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C RHO..... SUBROUTINE FOR COMPUTING LOSS FUNCTIONS AND THEIR DERIVS. C SEE RGLG FOR DETAILS ABOUT RHO. C RHOI.... PASSED WITHOUT CHANGE TO RHO. C RHOR.... PASSED WITHOUT CHANGE TO RHO. C IV...... INTEGER VALUES ARRAY. C LIV..... LENGTH OF IV, AT LEAST 90 + P. C LV...... LENGTH OF V, AT LEAST C 105 + P*(2*P + 17) + 4*PS C + N*(P + 5 + (P-PS+1)*(P-PS+2)/2). C V....... FLOATING-POINT VALUES ARRAY. C CALCRJ.. SUBROUTINE FOR COMPUTING RESIDUAL VECTOR. C UI...... PASSED UNCHANGED TO CALCRJ. C UR...... PASSED UNCHANGED TO CALCRJ. C UF...... PASSED UNCHANGED TO CALCRJ. C C *** CALCRJ CALLING SEQUENCE... C C CALL CALCRJ(N, PS, X, NF, NEED, R, RP, UI, UR, UF) C C PARAMETERS N, PS, X, UI, UR, AND UF ARE AS ABOVE. C R AND RP ARE FLOATING-POINT ARRAYS DIMENSIONED R(N) AND RP(PS,N). C NEED MAY BE REGARDED AS AN INTEGER THAT ALWAYS HAS THE VALUE 1 C WHEN GLF CALLS CALCRJ. THIS MEANS CALCRJ SHOULD COMPUTE THE C RESIDUAL VECTOR R. (CALCRJ SHOULD NOT CHANGE NEED OR RP. IF R C CANNOT BE COMPUTED, THEN CALCRJ SHOULD SET NF TO 0. OTHERWISE IT C SHOULD NOT CHANGE NF. FOR COMPATIBILITY WITH GLG, NEED IS A C VECTOR OF LENGTH 2.) C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ C C *** EXTERNAL SUBROUTINES *** C EXTERNAL IVSET, RGLG, V7CPY C C IVSET... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C RGLG... CARRIES OUT OPTIMIZATION ITERATIONS. C V7CPY... COPIES ONE VECTOR TO ANOTHER. C C *** LOCAL VARIABLES *** C INTEGER D1, DK, DR1, I, I1, IV1, J1K, J1K0, K, NEED(2), NF, 1 NG, RD1, R1, R21, RN, RS1 REAL H, H0, HLIM, NEGPT5, ONE, XK, ZERO C C *** IV AND V COMPONENTS *** C INTEGER COVREQ, D, DINIT, DLTFDJ, J, MODE, NEXTV, NFCALL, NFGCAL, 1 NGCALL, NGCOV, R, RDREQ, REGD, REGD0, TOOBIG, VNEED PARAMETER (COVREQ=15, D=27, DINIT=38, DLTFDJ=43, J=70, MODE=35, 1 NEXTV=47, NFCALL=6, NFGCAL=7, NGCALL=30, NGCOV=53, 2 R=61, RDREQ=57, REGD=67, REGD0=82, TOOBIG=2, VNEED=4) SAVE NEED DATA HLIM/0.1E+0/, NEGPT5/-0.5E+0/, ONE/1.E+0/, ZERO/0.E+0/ DATA NEED(1)/1/, NEED(2)/0/ C C--------------------------------- BODY ------------------------------ C IF (IV(1) .EQ. 0) CALL IVSET(1, IV, LIV, LV, V) IV(COVREQ) = -IABS(IV(COVREQ)) IF (IV(COVREQ) .EQ. 0 .AND. IV(RDREQ) .GT. 0) IV(COVREQ) = -1 IV1 = IV(1) IF (IV1 .EQ. 14) GO TO 10 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 IF (IV1 .EQ. 12) IV(1) = 13 I = (P-PS+2)*(P-PS+1)/2 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+3+I) CALL RGLG(X, V, IV, LIV, LV, N, PS, N, P, PS, V, V, RHO, RHOI, 1 RHOR, V, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(D) = IV(NEXTV) IV(R) = IV(D) + P IV(REGD0) = IV(R) + (P - PS + 3)*N IV(J) = IV(REGD0) + ((P-PS+2)*(P-PS+1)/2)*N IV(NEXTV) = IV(J) + N*PS IF (IV1 .EQ. 13) GO TO 999 C 10 D1 = IV(D) DR1 = IV(J) R1 = IV(R) RD1 = IV(REGD0) R21 = RD1 - N RS1 = R21 - N RN = RS1 + N - 1 C 20 CALL RGLG(V(D1), V(DR1), IV, LIV, LV, N, PS, N, P, PS, V(R1), 1 V(RD1), RHO, RHOI, RHOR, V, X) IF (IV(1)-2) 30, 50, 120 C C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** C 30 NF = IV(NFCALL) CALL CALCRJ(N, PS, X, NF, NEED, V(R1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 40 IV(TOOBIG) = 1 GO TO 20 40 CALL V7CPY(N, V(RS1), V(R1)) IF (IV(1) .GT. 0) GO TO 20 C C *** COMPUTE FINITE-DIFFERENCE APPROXIMATION TO DR = GRAD. OF R *** C C *** INITIALIZE D IF NECESSARY *** C 50 IF (IV(MODE) .LT. 0 .AND. V(DINIT) .EQ. ZERO) 1 CALL V7SCP(P, V(D1), ONE) C DK = D1 NG = IV(NGCALL) - 1 IF (IV(1) .EQ. (-1)) IV(NGCOV) = IV(NGCOV) - 1 J1K0 = DR1 NF = IV(NFCALL) IF (NF .EQ. IV(NFGCAL)) GO TO 70 NG = NG + 1 CALL CALCRJ(N, PS, X, NF, NEED, V(RS1), V(DR1), UI, UR, UF) IF (NF .GT. 0) GO TO 70 60 IV(TOOBIG) = 1 IV(NGCALL) = NG GO TO 20 70 DO 110 K = 1, PS XK = X(K) H = V(DLTFDJ) * MAX( ABS(XK), ONE/V(DK)) H0 = H DK = DK + 1 80 X(K) = XK + H NG = NG + 1 NF = -NG CALL CALCRJ(N, PS, X, NF, NEED, V(R21), V(DR1), UI, UR, UF) IF (NF .LT. 0) GO TO 90 H = NEGPT5 * H IF ( ABS(H/H0) .GE. HLIM) GO TO 80 GO TO 60 90 X(K) = XK IV(NGCALL) = NG I1 = R21 J1K = J1K0 J1K0 = J1K0 + 1 DO 100 I = RS1, RN V(J1K) = (V(I1) - V(I)) / H I1 = I1 + 1 J1K = J1K + PS 100 CONTINUE 110 CONTINUE GO TO 20 C 120 IF (IV(REGD) .GT. 0) IV(REGD) = RD1 C 999 RETURN C C *** LAST LINE OF GLF FOLLOWS *** END SUBROUTINE RGLG(D, DR, IV, LIV, LV, N, ND, NN, P, PS, R, 1 RD, RHO, RHOI, RHOR, V, X) C C *** ITERATION DRIVER FOR GENERALIZED (NON)LINEAR MODELS (ETC.) C INTEGER LIV, LV, N, ND, NN, P, PS INTEGER IV(LIV), RHOI(*) REAL D(P), DR(ND,N), R(*), RD(*), RHOR(*), 1 V(LV), X(*) C DIMENSION RD(N, (P-PS)*(P-PS+1)/2 + 1) EXTERNAL RHO C C-------------------------- PARAMETER USAGE -------------------------- C C D....... SCALE VECTOR. C DR...... DERIVATIVES OF R AT X. C IV...... INTEGER VALUES ARRAY. C LIV..... LENGTH OF IV... LIV MUST BE AT LEAST P + 90. C LV...... LENGTH OF V... LV MUST BE AT LEAST C 105 + P*(2*P+16) + 2*N + 4*PS. C N....... TOTAL NUMBER OF RESIDUALS. C ND...... LEADING DIMENSION OF DR -- MUST BE AT LEAST PS. C NN...... LEAD DIMENSION OF R, RD. C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C PS...... NUMBER OF NON-NUISANCE PARAMETERS. C R....... RESIDUALS (OR MEANS -- FUNCTIONS OF X) WHEN RGLG IS CALLED C WITH IV(1) = 1. C RD...... RD(I) = HALF * (G(I)**T * H(I)**-1 * G(I)) ON OUTPUT WHEN C IV(RDREQ) IS 2, 3, 4, OR 5. RGLG SETS IV(REGD) = 1 IF RD C IS SUCCESSFULLY COMPUTED, TO 0 IF NO ATTEMPT WAS MADE C TO COMPUTE IT, AND TO -1 IF H (THE FINITE-DIFFERENCE HESSIAN) C WAS INDEFINITE. BEFORE CONVERGENCE, RD IS ALSO USED AS C TEMPORARY STORAGE. C RHO..... COMPUTES INFO ABOUT OBJECTIVE FUNCTION. C RHOI.... PASSED WITHOUT CHANGE TO RHO. C RHOR.... PASSED WITHOUT CHANGE TO RHO. C V....... FLOATING-POINT VALUES ARRAY. C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C C *** CALLING SEQUENCE FOR RHO... C C CALL RHO(NEED, F, N, NF, XN, R, RD, RHOI, RHOR, W) C C PARAMETER DECLARATIONS FOR RHO... C C INTEGER NEED(2), N, NF, RHOI(*) C FLOATING-POINT F, XN(*), R(*), RD(N,*), RHOR(*), W(N) C C RHOI AND RHOR ARE FOR RHO TO USE AS IT SEES FIT. THEY ARE PASSED C TO RHO WITHOUT CHANGE. IF IV(RDREQ) IS AT LEAST 4 (MOD 6), I.E., IF C MORE THAN THE SIMPLEST REGRESSION DIAGNOSTIC INFORMATION IS TO BE C COMPUTED, THEN SOME COMPONENTS OF RHOI AND RHOR MUST CONVEY SOME EXTRA C DETAILS, AS DESCRIBED BELOW. C F, R, RD, AND W ARE EXPLAINED BELOW WITH NEED. C XN IS THE VECTOR OF NUISANCE PARAMETERS (OF LENGTH P - PS). IF C RHO NEEDS TO KNOW THE LENGTH OF XN, THEN THIS LENGTH SHOULD BE C COMMUNICATED THROUGH RHOI (OR THROUGH COMMON). RHO SHOULD NOT CHANGE C XN. C NEED(1) = 1 MEANS RHO SHOULD SET F TO THE SUM OF THE LOSS FUNCTION C VALUES AT THE RESIDUALS R(I). NF IS THE CURRENT FUNCTION INVOCATION C COUNT (A VALUE THAT IS INCREMENTED EACH TIME A NEW PARAMETER ESTIMATE C X IS CONSIDERED). NEED(2) IS THE VALUE NF HAD AT THE LAST R WHERE C RHO MIGHT BE CALLED WITH NEED(1) = 2. IF RHO SAVES INTERMEDIATE C RESULTS FOR USE IN CALLS WITH NEED(1) = 2, THEN IT CAN USE NF TO TELL C WHICH INTERMEDIATE RESULTS ARE APPROPRIATE, AND IT CAN SAVE SOME OF C THESE RESULTS IN R. C NEED(1) = 2 MEANS RHO SHOULD SET R(I) TO THE LOSS FUNCTION C DERIVATIVE WITH RESPECT TO THE RESIDUALS THAT WERE PASSED TO RHO WHEN C NF HAD THE SAME VALUE IT DOES NOW (AND NEED(1) WAS 1). RHO SHOULD C ALSO SET W(I) TO THE APPROXIMATION OF THE SECOND DERIVATIVE OF THE C LOSS FUNCTION (WITH RESPECT TO THE I-TH RESIDUAL) THAT SHOULD BE USED C IN THE GAUSS-NEWTON MODEL. WHEN THERE ARE NUISANCE PARAMETERS (I.E., C WHEN PS .LT. P) RHO SHOULD ALSO SET R(I+K*N) TO THE DERIVATIVE OF THE C LOSS FUNCTION WITH RESPECT TO THE I-TH RESIDUAL AND XN(K), AND IT C SHOULD SET RD(I,J + K*(K+1)/2 + 1) TO THE SECOND PARTIAL DERIVATIVE C OF THE I-TH RESIDUAL WITH RESPECT TO XN(J) AND XN(K), 0 .LE. J .LE. K C AND 1 .LE. K .LE. P - PS, WHERE XN(0) MEANS THE I-TH RESIDUAL ITSELF. C IN ANY EVENT, RHO SHOULD ALSO SET RD(I,1) TO THE (TRUE) SECOND C DERIVATIVE OF THE LOSS FUNCTION WITH RESPECT TO THE I-TH RESIDUAL. C NF (THE FUNCTION INVOCATION COUNT WHOSE NORMAL USE IS EXPLAINED C ABOVE) SHOULD NOT BE CHANGED UNLESS RHO CANNOT CARRY OUT THE REQUESTED C TASK, IN WHICH CASE RHO SHOULD SET NF TO 0. C C C *** REGRESSION DIAGNOSTICS *** C C IV(RDREQ) INDICATES WHETHER A COVARIANCE MATRIX AND REGRESSION C DIAGNOSTIC VECTOR ARE TO BE COMPUTED. IV(RDREQ) (MOD 6) HAS THE FORM C IV(RDREQ) = CVR +2*RDR, WHERE CVR = 0 OR 1 AND RDR = 0, 1, OR 2, C SO THAT C C CVR = MOD(IV(RDREQ), 2) C RDR = MOD(IV(RDREQ)/2, 3). C C CVR = 0 FOR NO COVARIANCE MATRIX C = 1 IF A COVARIANCE MATRIX ESTIMATE IS DESIRED C C RDR = 0 FOR NO LEAVE-ONE-OUT DIAGNOSTIC INFORMATION. C = 1 TO HAVE ONE-STEP ESTIMATES OF F(X(I)) - F(X*) STORED IN RD, C WHERE X(I) MINIMIZES F (THE OBJECTIVE FUNCTION) WITH C COMPONENT I OF R REMOVED AND X* MINIMIZES THE FULL F. C = 2 FOR MORE DETAILED ONE-STEP LEAVE-ONE-OUT INFORMATION, AS C DICTATED BY THE IV COMPONENTS DESCRIBED BELOW. C C FOR RDR = 2, THE FOLLOWING COMPONENTS OF IV ARE RELEVANT... C C NFIX = IV(83) = NUMBER OF TRAILING NUISANCE PARAMETERS TO TREAT C AS FIXED WHEN COMPUTING DIAGNOSTIC VECTORS (0 .LE. NFIX .LE. C P - PS, SO X(I) IS KEPT FIXED FOR P - NFIX .LT. I .LE. P). C C LOO = IV(84) TELLS WHAT TO LEAVE OUT... C = 1 MEANS LEAVE OUT EACH COMPONENT OF R SEPARATELY, AND C = 2 MEANS LEAVE OUT CONTIGUOUS BLOCKS OF R COMPONENTS. C FOR LOO = 2, IV(85) IS THE STARTING SUBSCRIPT IN RHOI C OF AN ARRAY BS OF BLOCK SIZES, IV(86) IS THE STRIDE FOR BS, C AND IV(87) = NB IS THE NUMBER OF BLOCKS, SO THAT C BS(I) = RHOI(IV(85) + (I-1)*IV(86)), 1 .LE. I .LE. NB. C NOTE THAT IF ALL BLOCKS ARE THE SAME SIZE, THEN IT SUFFICES C TO SET RHOI(IV(85)) = BLOCKSIZE AND IV(86) = 0. C NOTE THAT LOO = 1 IS EQUIVALENT TO LOO = 2 WITH C RHOI(IV(85)) = 1, IV(86) = 0, IV(87) = N. C = 3,4 ARE SIMILAR TO LOO = 1,2, RESPECTIVELY, BUT LEAVING A C FRACTION OUT. IN THIS CASE, IV(88) IS THE STARTING C SUBSCRIPT IN RHOR OF AN ARRAY FLO OF FRACTIONS TO LEAVE OUT, C AND IV(89) IS THE STRIDE FOR FLO... C FLO(I) = RHOR(IV(88) + (I-1)*IV(89)), 1 .LE. I .LE. NB. C C XNOTI = IV(90) TELLS WHAT DIAGNOSTIC INFORMATION TO STORE... C = 0 MEANS JUST STORE ONE-STEP ESTIMATES OF F(X(I)) - F(X*) IN C RD(I), 1 .LE. I .LE. NB. C .GT. 0 MEANS ALSO STORE ONE-STEP ESTIMATES OF X(I) ESTIMATES C IN RHOR, STARTING AT RHOR(XNOTI)... C X(I)(J) = RHOR((I-1)*(P-NFIX) + J + XNOTI-1), C 1 .LE. I .LE. NB, 1 .LE. J .LE. P - NFIX. C C SOMETIMES ONE-STEP ESTIMATES OF X(I) DO NOT EXIST, BECAUSE THE C APPROXIMATE UPDATED HESSIAN MATRIX IS INDEFINITE. IN SUCH CASES, C THE CORRESPONDING RD COMPONENT IS SET TO -1, AND, IF XNOTI IS C POSITIVE, THE SOLUTION X IS RETURNED AS X(I). WHEN ONE-STEP ESTIMATES C OF X(I) DO EXIST, THE CORRESPONDING COMPONENT OF RD IS POSITIVE. C C SUMMARY OF RHOI COMPONENTS (FOR RDR = MOD(IV(RDREQ)/2, 3) = 2)... C C IV(83) = NFIX C IV(84) = LOO C IV(85) = START IN RHOI OF BS C IV(86) = STRIDE FOR BS C IV(87) = NB C IV(88) = START IN RHOR OF FLO C IV(89) = STRIDE FOR FLO C IV(90) = XNOTI (START IN RHOR OF X(I)). C C C *** COVARIANCE MATRIX ESTIMATE *** C C IF IV(RDREQ) INDICATES THAT A COVARIANCE MATRIX IS TO BE COMPUTED, C THEN IV(COVREQ) = IV(15) DETERMINES THE FORM OF THE COMPUTED C COVARIANCE MATRIX ESTIMATE AND, SIMULTANEOUSLY, THE FORM OF C APPROXIMATE HESSIAN MATRIX USED IN COMPUTING REGRESSION DIAGNOSTIC C INFORMATION. IN ALL CASES, SOME APPROXIMATE FINAL HESSIAN MATRIX C IS OBTAINED, AND ITS INVERSE IS THE COVARIANCE MATRIX ESTIMATE C (WHICH MAY HAVE TO BE SCALED APPROPRIATELY -- THAT IS UP TO YOU). C IF IV(COVREQ) IS AT MOST 2 IN ABSOLUTE VALUE, THEN THE FINAL C HESSIAN APPROXIMATION IS COMPUTED BY FINITE DIFFERENCES -- GRADIENT C DIFFERENCES IF IV(COVREQ) IS NONNEGATIVE, FUNCTION DIFFERENCES C OTHERWISE. IF (IV(COVREQ)) IS AT LEAST 3 IN ABSOLUTE VALUE, THEN THE C CURRENT GAUSS-NEWTON HESSIAN APPROXIMATION IS TAKEN AS THE FINAL C HESSIAN APPROXIMATION. FOR SOME PROBLEMS THIS SAVES TIME AND YIELDS C THE SAME OR NEARLY THE SAME HESSIAN APPROXIMATION AS DO FINITE C DIFFERENCES. FOR OTHER PROBLEMS, THE TWO KINDS OF HESSIAN C APPROXIMATIONS MAY GIVE DECIDEDLY DIFFERENT REGRESSION DIAGNOSTICS AND C COVARIANCE MATRIX ESTIMATES. C C SET IV(RDREQ) TO 2, 4 OR 6 TO REQUEST A CHOLESKY FACTOR FOR THE FINAL C APPROXIMATE HESSIAN MATRIX. IF IV(FDH) IS POSITIVE WHEN THIS ROUTINE C RETURNS, THEN THE CHOLESKY FACTOR IS STORED IN V, STARTING AT C V(IV(FDH)), IN THE SEQUENCE L(1,1), L(2,1), L(2,2), L(3,1), ... . C C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C+++++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++ C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL D7UP5, IVSET, G2LRD, N3RDP, D7TPR, Q7ADR, VSUM, 1 G7LIT, ITSUM, L7NVR, L7ITV, L7IVM, L7SRT, L7SQR, 2 L7SVX, L7SVN, L7TSQ, L7VML, O7PRD, V2AXY, V7CPY, 3 V7SCL, V7SCP REAL D7TPR, L7SVX, L7SVN, VSUM C C D7UP5... UPDATES SCALE VECTOR D. C IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C G2LRD.... COMPUTES REGRESSION DIAGNOSTIC. C N3RDP... PRINTS REGRESSION DIAGNOSTIC. C D7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. C Q7ADR.... ADDS ROWS TO QR FACTORIZATION. C VSUM..... RETURNS SUM OF ELEMENTS OF A VECTOR. C G7LIT.... PERFORMS BASIC MINIMIZATION ALGORITHM. C ITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X. C L7NVR... INVERTS COMPACTLY STORED TRIANGULAR MATRIX. C L7ITV... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. C L7IVM... APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. C L7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. C L7SQR... COMPUTES L*(L**T) FOR LOWER TRIANG. MATRIX L. C L7SVX... UNDERESTIMATES LARGEST SINGULAR VALUE OF TRIANG. MATRIX. C L7SVN... OVERESTIMATES SMALLEST SINGULAR VALUE OF TRIANG. MATRIX. C L7TSQ... COMPUTES (L**T)*L FOR LOWER TRIANG. MATRIX L. C L7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C O7PRD.... ADDS OUTER PRODUCT OF VECTORS TO A MATRIX. C V2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER. C V7CPY.... COPIES ONE VECTOR TO ANOTHER. C V7SCL... MULTIPLIES A VECTOR BY A SCALAR. C V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C C *** LOCAL VARIABLES *** C LOGICAL JUSTG, UPDATD, ZEROG INTEGER G1, HN1, I, II, IV1, J, J1, JTOL1, K, L, LH, 1 NEED1(2), NEED2(2), PMPS, PS1, PSLEN, QTR1, 2 RMAT1, STEP1, TEMP1, TEMP2, TEMP3, TEMP4, W, WI, Y1 REAL RHMAX, RHTOL, RHO1, RHO2, T C REAL ONE, ZERO C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, COVMAT, DINIT, DTYPE, DTINIT, D0INIT, F, 1 F0, FDH, G, H, HC, IPIVOT, IVNEED, JCN, JTOL, LMAT, 2 MODE, NEXTIV, NEXTV, NF0, NF1, NFCALL, NFCOV, NFGCAL, 3 NGCALL, NGCOV, PERM, QTR, RDREQ, REGD, RESTOR, 4 RMAT, RSPTOL, STEP, TOOBIG, VNEED, XNOTI, Y C C *** IV SUBSCRIPT VALUES *** C PARAMETER (CNVCOD=55, COVMAT=26, DTYPE=16, F0=13, FDH=74, G=28, 1 H=56, HC=71, IPIVOT=76, IVNEED=3, JCN=66, JTOL=59, 2 LMAT=42, MODE=35, NEXTIV=46, NEXTV=47, NFCALL=6, 3 NFCOV=52, NF0=68, NF1=69, NFGCAL=7, NGCALL=30, 4 NGCOV=53, PERM=58, QTR=77, RESTOR=9, RMAT=78, RDREQ=57, 5 REGD=67, STEP=40, TOOBIG=2, VNEED=4, XNOTI=90, Y=48) C C *** V SUBSCRIPT VALUES *** C PARAMETER (DINIT=38, DTINIT=39, D0INIT=40, F=10, RSPTOL=49) PARAMETER (ONE=1.E+0, ZERO=0.E+0) SAVE NEED1, NEED2 DATA NEED1(1)/1/, NEED1(2)/0/, NEED2(1)/2/, NEED2(2)/0/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C LH = P * (P+1) / 2 IF (IV(1) .EQ. 0) CALL IVSET(1, IV, LIV, LV, V) PS1 = PS + 1 IV1 = IV(1) IF (IV1 .GT. 2) GO TO 10 W = IV(Y) + P IV(RESTOR) = 0 I = IV1 + 2 IF (IV(TOOBIG) .EQ. 0) GO TO (120, 110, 110, 130), I V(F) = V(F0) IF (I .NE. 3) IV(1) = 2 GO TO 40 C C *** FRESH START OR RESTART -- CHECK INPUT INTEGERS *** C 10 IF (ND .LT. PS) GO TO 360 IF (PS .GT. P) GO TO 360 IF (PS .LE. 0) GO TO 360 IF (N .LE. 0) GO TO 360 IF (IV1 .EQ. 14) GO TO 30 IF (IV1 .GT. 16) GO TO 420 IF (IV1 .LT. 12) GO TO 40 IF (IV1 .EQ. 12) IV(1) = 13 IF (IV(1) .NE. 13) GO TO 20 IV(IVNEED) = IV(IVNEED) + P IV(VNEED) = IV(VNEED) + P*(P+13)/2 + 2*N + 4*PS C *** ADJUST IV(PERM) TO MAKE ROOM FOR IV INPUT COMPONENTS C *** NEEDED WHEN IV(RDREQ) IS 4 OR 5... I = XNOTI + 1 IF (IV(PERM) .LT. I) IV(PERM) = I C 20 CALL G7LIT(D, X, IV, LIV, LV, P, PS, V, X, X) IF (IV(1) .NE. 14) GO TO 999 C C *** STORAGE ALLOCATION *** C IV(IPIVOT) = IV(NEXTIV) IV(NEXTIV) = IV(IPIVOT) + P IV(Y) = IV(NEXTV) IV(G) = IV(Y) + P + N IV(RMAT) = IV(G) + P + 4*PS IV(QTR) = IV(RMAT) + LH IV(JTOL) = IV(QTR) + P + N IV(JCN) = IV(JTOL) + 2*P IV(NEXTV) = IV(JCN) + P IF (IV1 .EQ. 13) GO TO 999 C 30 JTOL1 = IV(JTOL) IF (V(DINIT) .GE. ZERO) CALL V7SCP(P, D, V(DINIT)) IF (V(DTINIT) .GT. ZERO) CALL V7SCP(P, V(JTOL1), V(DTINIT)) I = JTOL1 + P IF (V(D0INIT) .GT. ZERO) CALL V7SCP(P, V(I), V(D0INIT)) IV(NF0) = 0 IV(NF1) = 0 C 40 G1 = IV(G) Y1 = IV(Y) CALL G7LIT(D, V(G1), IV, LIV, LV, P, PS, V, X, V(Y1)) IF (IV(1) - 2) 50, 60, 380 C 50 V(F) = ZERO IF (IV(NF1) .EQ. 0) GO TO 999 IF (IV(RESTOR) .NE. 2) GO TO 999 IV(NF0) = IV(NF1) CALL V7CPY(N, RD, R) IV(REGD) = 0 GO TO 999 C 60 IF (IV(MODE) .GT. 0) GO TO 370 CALL V7SCP(P, V(G1), ZERO) RMAT1 = IABS(IV(RMAT)) QTR1 = IABS(IV(QTR)) CALL V7SCP(PS, V(QTR1), ZERO) IV(REGD) = 0 CALL V7SCP(PS, V(Y1), ZERO) CALL V7SCP(LH, V(RMAT1), ZERO) IF (IV(RESTOR) .NE. 3) GO TO 70 CALL V7CPY(N, R, RD) IV(NF1) = IV(NF0) 70 CALL RHO(NEED2, T, N, IV(NFGCAL), X(PS1), R, RD, RHOI, RHOR, V(W)) IF (IV(NFGCAL) .GT. 0) GO TO 90 80 IV(TOOBIG) = 1 GO TO 40 90 IF (IV(MODE) .LT. 0) GO TO 999 DO 100 I = 1, N 100 CALL V2AXY(PS, V(Y1), R(I), DR(1,I), V(Y1)) GO TO 999 C C *** COMPUTE F(X) *** C 110 I = IV(NFCALL) NEED1(2) = IV(NFGCAL) CALL RHO(NEED1, V(F), N, I, X(PS1), R, RD, RHOI, RHOR, V(W)) IV(NF1) = I IF (I .LE. 0) GO TO 80 GO TO 40 C C *** COMPUTE GRADIENT INFORMATION FOR FINITE-DIFFERENCE HESSIAN *** C 120 IV(1) = 2 JUSTG = .TRUE. I = IV(NFCALL) CALL RHO(NEED1, T, N, I, X(PS1), R, RD, RHOI, RHOR, V(W)) IF (I .LE. 0) GO TO 80 CALL RHO(NEED2, T, N, I, X(PS1), R, RD, RHOI, RHOR, V(W)) IF (I .LE. 0) GO TO 80 GO TO 250 C C *** PREPARE TO COMPUTE GRADIENT INFORMATION WHILE ITERATING *** C 130 JUSTG = .FALSE. G1 = IV(G) C C *** DECIDE WHETHER TO UPDATE D BELOW *** C I = IV(DTYPE) UPDATD = .FALSE. IF (I .LE. 0) GO TO 140 IF (I .EQ. 1 .OR. IV(MODE) .LT. 0) UPDATD = .TRUE. C C *** COMPUTE RMAT AND QTR *** C 140 QTR1 = IABS(IV(QTR)) RMAT1 = IABS(IV(RMAT)) IV(RMAT) = RMAT1 IV(HC) = 0 IV(NF0) = 0 IV(NF1) = 0 IF (IV(MODE) .LT. 0) GO TO 160 C C *** ADJUST Y *** C Y1 = IV(Y) WI = W STEP1 = IV(STEP) DO 150 I = 1, N T = V(WI) - RD(I) WI = WI + 1 IF (T .NE. ZERO) CALL V2AXY(PS, V(Y1), 1 T* D7TPR(PS,V(STEP1),DR(1,I)), DR(1,I), V(Y1)) 150 CONTINUE C C *** CHECK FOR NEGATIVE W COMPONENTS *** C 160 J1 = W + N - 1 DO 170 WI = W, J1 IF (V(WI) .LT. ZERO) GO TO 240 170 CONTINUE C C *** W IS NONNEGATIVE. COMPUTE QR FACTORIZATION *** C *** AND, IF NECESSARY, USE SEMINORMAL EQUATIONS *** C RHMAX = ZERO RHTOL = V(RSPTOL) TEMP1 = G1 + P ZEROG = .TRUE. WI = W DO 200 I = 1, N RHO1 = R(I) RHO2 = V(WI) WI = WI + 1 T = SQRT(RHO2) IF (RHMAX .LT. RHO2) RHMAX = RHO2 IF (RHO2 .GT. RHTOL*RHMAX) GO TO 180 C *** SEMINORMAL EQUATIONS *** CALL V2AXY(PS, V(G1), RHO1, DR(1,I), V(G1)) RHO1 = ZERO ZEROG = .FALSE. GO TO 190 180 RHO1 = RHO1 / T C *** QR ACCUMULATION *** 190 CALL V7SCL(PS, V(TEMP1), T, DR(1,I)) CALL Q7ADR(PS, V(QTR1), V(RMAT1), V(TEMP1), RHO1) 200 CONTINUE C C *** COMPUTE G FROM RMAT AND QTR *** C TEMP2 = TEMP1 + PS CALL L7VML(PS, V(TEMP1), V(RMAT1), V(QTR1)) IF (ZEROG) GO TO 220 IV(QTR) = -QTR1 IF ( L7SVX(PS, V(RMAT1), V(TEMP2), V(TEMP2)) * RHTOL .GE. 1 L7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2))) GO TO 230 CALL L7IVM(PS, V(TEMP2), V(RMAT1), V(G1)) C C *** SEMINORMAL EQUATIONS CORRECTION OF BJOERCK -- C *** ONE CYCLE OF ITERATIVE REFINEMENT... C TEMP3 = TEMP2 + PS TEMP4 = TEMP3 + PS CALL L7ITV(PS, V(TEMP3), V(RMAT1), V(TEMP2)) CALL V7SCP(PS, V(TEMP4), ZERO) RHMAX = ZERO WI = W DO 210 I = 1, N RHO2 = V(WI) WI = WI + 1 IF (RHMAX .LT. RHO2) RHMAX = RHO2 RHO1 = ZERO IF (RHO2 .LE. RHTOL*RHMAX) RHO1 = R(I) T = RHO1 - RHO2* D7TPR(PS, V(TEMP3), DR(1,I)) CALL V2AXY(PS, V(TEMP4), T, DR(1,I), V(TEMP4)) 210 CONTINUE CALL L7IVM(PS, V(TEMP3), V(RMAT1), V(TEMP4)) CALL V2AXY(PS, V(TEMP2), ONE, V(TEMP3), V(TEMP2)) CALL V2AXY(PS, V(QTR1), ONE, V(TEMP2), V(QTR1)) 220 IV(QTR) = QTR1 230 CALL V2AXY(PS, V(G1), ONE, V(TEMP1), V(G1)) IF (PS .GE. P) GO TO 350 GO TO 270 C C *** INDEFINITE GN HESSIAN... *** C 240 IV(RMAT) = -RMAT1 IV(HC) = RMAT1 CALL O7PRD(N, LH, PS, V(RMAT1), V(W), DR, DR) C C *** COMPUTE GRADIENT *** C 250 G1 = IV(G) CALL V7SCP(P, V(G1), ZERO) DO 260 I = 1, N 260 CALL V2AXY(PS, V(G1), R(I), DR(1,I), V(G1)) IF (PS .GE. P) GO TO 350 C C *** COMPUTE GRADIENT COMPONENTS OF NUISANCE PARAMETERS *** C 270 K = P - PS J1 = 1 G1 = G1 + PS DO 280 J = 1, K J1 = J1 + NN V(G1) = VSUM(N, R(J1)) G1 = G1 + 1 280 CONTINUE IF (JUSTG) GO TO 390 C C *** COMPUTE HESSIAN COMPONENTS OF NUISANCE PARAMETERS *** C I = PS*PS1/2 PSLEN = P*(P+1)/2 - I HN1 = RMAT1 + I CALL V7SCP(PSLEN, V(HN1), ZERO) PMPS = P - PS K = HN1 J1 = 1 DO 310 II = 1, PMPS J1 = J1 + NN J = J1 DO 290 I = 1, N CALL V2AXY(PS, V(K), RD(J), DR(1,I), V(K)) J = J + 1 290 CONTINUE K = K + PS DO 300 I = 1, II J1 = J1 + NN V(K) = VSUM(N, RD(J1)) K = K + 1 300 CONTINUE 310 CONTINUE IF (IV(RMAT) .LE. 0) GO TO 350 J = IV(LMAT) CALL V7CPY(PSLEN, V(J), V(HN1)) IF ( L7SVN(PS, V(RMAT1), V(TEMP2), V(TEMP2)) .LE. ZERO) GO TO 320 CALL L7SRT(PS1, P, V(RMAT1), V(RMAT1), I) IF (I .LE. 0) GO TO 330 C C *** HESSIAN IS NOT POSITIVE DEFINITE *** C 320 CALL L7SQR(PS, V(RMAT1), V(RMAT1)) CALL V7CPY(PSLEN, V(HN1), V(J)) IV(HC) = RMAT1 IV(RMAT) = -RMAT1 GO TO 350 C C *** NUISANCE PARS LEAVE HESSIAN POS. DEF. GET REST OF QTR *** C 330 J = QTR1 + PS G1 = IV(G) + PS DO 340 I = PS1, P T = D7TPR(I-1, V(HN1), V(QTR1)) HN1 = HN1 + I V(J) = (V(G1) - T) / V(HN1-1) J = J + 1 G1 = G1 + 1 340 CONTINUE 350 IF (JUSTG) GO TO 390 IF (UPDATD) CALL D7UP5(D, IV, LIV, LV, P, PS, V) GO TO 40 C C *** MISC. DETAILS *** C C *** BAD N, ND, OR P *** C 360 IV(1) = 66 GO TO 420 C C *** COVARIANCE OR INITIAL S COMPUTATION *** C 370 IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(NFGCAL) = IV(NFCALL) IV(1) = -1 GO TO 999 C C *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** C 380 IF (IV(COVMAT) .NE. 0) GO TO 410 IF (IV(REGD) .NE. 0) GO TO 410 C C *** SEE IF CHOLESKY FACTOR OF HESSIAN IS AVAILABLE *** C K = IV(FDH) IF (K .LE. 0) GO TO 400 IF (IV(RDREQ) .LE. 0) GO TO 410 C C *** COMPUTE REGRESSION DIAGNOSTICS AND DEFAULT COVARIANCE IF C DESIRED *** C IV(MODE) = P + 1 IV(NGCALL) = IV(NGCALL) + 1 IV(NGCOV) = IV(NGCOV) + 1 IV(CNVCOD) = IV(1) IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(NFGCAL) = IV(NFCALL) IV(1) = -1 GO TO 999 C 390 IF (IV(MODE) .LE. P) GO TO 40 C *** SAVE RD IN W FOR POSSIBLE USE IN OTHER DIAGNOSTICS *** CALL V7CPY(N, V(W), RD) C *** OVERWRITE RD WITH REGRESSION DIAGNOSTICS *** L = IV(LMAT) I = IV(JCN) STEP1 = IV(STEP) CALL G2LRD(DR, IV, V(L), LH, LIV, LV, ND, N, P, PS, R, RD, 1 RHOI, RHOR, V, V(STEP1), X, V(I)) IV(1) = IV(CNVCOD) IV(CNVCOD) = 0 IF (MOD(IV(RDREQ),2) .EQ. 0) GO TO 410 C C *** FINISH COVARIANCE COMPUTATION *** C I = IABS(IV(H)) IV(FDH) = 0 CALL L7NVR(P, V(I), V(L)) CALL L7TSQ(P, V(I), V(I)) IV(COVMAT) = I GO TO 410 C C *** COME HERE FOR INDEFINITE FINITE-DIFFERENCE HESSIAN *** C 400 IV(COVMAT) = K IV(REGD) = K C C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** C 410 G1 = IV(G) 420 CALL ITSUM(D, V(G1), IV, LIV, LV, P, V, X) IF (IV(1) .LE. 6 .AND. IV(RDREQ) .GT. 0) 1 CALL N3RDP(IV, LIV, LV, N, P, RD, RHOI, RHOR, V) C 999 RETURN C *** LAST LINE OF RGLG FOLLOWS *** END SUBROUTINE F7HES(D, G, IRT, IV, LIV, LV, P, V, X) C C *** COMPUTE FINITE-DIFFERENCE HESSIAN, STORE IT IN V STARTING C *** AT V(IV(FDH)) = V(-IV(H)). C C *** IF IV(COVREQ) .GE. 0 THEN F7HES USES GRADIENT DIFFERENCES, C *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN G7LIT. C C IRT VALUES... C 1 = COMPUTE FUNCTION VALUE, I.E., V(F). C 2 = COMPUTE G. C 3 = DONE. C C C *** PARAMETER DECLARATIONS *** C INTEGER IRT, LIV, LV, P INTEGER IV(LIV) REAL D(P), G(P), V(LV), X(P) C C *** LOCAL VARIABLES *** C INTEGER GSAVE1, HES, HMI, HPI, HPM, I, K, KIND, L, M, MM1, MM1O2, 1 PP1O2, STPI, STPM, STP0 REAL DEL, HALF, NEGPT5, ONE, TWO, ZERO C C *** EXTERNAL SUBROUTINES *** C EXTERNAL V7CPY C C V7CPY.... COPY ONE VECTOR TO ANOTHER. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER COVREQ, DELTA, DELTA0, DLTFDC, F, FDH, FX, H, KAGQT, MODE, 1 NFGCAL, SAVEI, SWITCH, TOOBIG, W, XMSAVE C PARAMETER (HALF=0.5E+0, NEGPT5=-0.5E+0, ONE=1.E+0, TWO=2.E+0, 1 ZERO=0.E+0) C PARAMETER (COVREQ=15, DELTA=52, DELTA0=44, DLTFDC=42, F=10, 1 FDH=74, FX=53, H=56, KAGQT=33, MODE=35, NFGCAL=7, 2 SAVEI=63, SWITCH=12, TOOBIG=2, W=65, XMSAVE=51) C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C IRT = 4 KIND = IV(COVREQ) M = IV(MODE) IF (M .GT. 0) GO TO 10 IV(H) = -IABS(IV(H)) IV(FDH) = 0 IV(KAGQT) = -1 V(FX) = V(F) 10 IF (M .GT. P) GO TO 999 IF (KIND .LT. 0) GO TO 110 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING BOTH FUNCTION AND C *** GRADIENT VALUES. C GSAVE1 = IV(W) + P IF (M .GT. 0) GO TO 20 C *** FIRST CALL ON F7HES. SET GSAVE = G, TAKE FIRST STEP *** CALL V7CPY(P, V(GSAVE1), G) IV(SWITCH) = IV(NFGCAL) GO TO 90 C 20 DEL = V(DELTA) X(M) = V(XMSAVE) IF (IV(TOOBIG) .EQ. 0) GO TO 40 C C *** HANDLE OVERSIZE V(DELTA) *** C IF (DEL*X(M) .GT. ZERO) GO TO 30 C *** WE ALREADY TRIED SHRINKING V(DELTA), SO QUIT *** IV(FDH) = -2 GO TO 220 C C *** TRY SHRINKING V(DELTA) *** 30 DEL = NEGPT5 * DEL GO TO 100 C 40 HES = -IV(H) C C *** SET G = (G - GSAVE)/DEL *** C DO 50 I = 1, P G(I) = (G(I) - V(GSAVE1)) / DEL GSAVE1 = GSAVE1 + 1 50 CONTINUE C C *** ADD G AS NEW COL. TO FINITE-DIFF. HESSIAN MATRIX *** C K = HES + M*(M-1)/2 L = K + M - 2 IF (M .EQ. 1) GO TO 70 C C *** SET H(I,M) = 0.5 * (H(I,M) + G(I)) FOR I = 1 TO M-1 *** C MM1 = M - 1 DO 60 I = 1, MM1 V(K) = HALF * (V(K) + G(I)) K = K + 1 60 CONTINUE C C *** ADD H(I,M) = G(I) FOR I = M TO P *** C 70 L = L + 1 DO 80 I = M, P V(L) = G(I) L = L + I 80 CONTINUE C 90 M = M + 1 IV(MODE) = M IF (M .GT. P) GO TO 210 C C *** CHOOSE NEXT FINITE-DIFFERENCE STEP, RETURN TO GET G THERE *** C DEL = V(DELTA0) * MAX(ONE/D(M), ABS(X(M))) IF (X(M) .LT. ZERO) DEL = -DEL V(XMSAVE) = X(M) 100 X(M) = X(M) + DEL V(DELTA) = DEL IRT = 2 GO TO 999 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN USING FUNCTION VALUES ONLY. C 110 STP0 = IV(W) + P - 1 MM1 = M - 1 MM1O2 = M*MM1/2 IF (M .GT. 0) GO TO 120 C *** FIRST CALL ON F7HES. *** IV(SAVEI) = 0 GO TO 200 C 120 I = IV(SAVEI) HES = -IV(H) IF (I .GT. 0) GO TO 180 IF (IV(TOOBIG) .EQ. 0) GO TO 140 C C *** HANDLE OVERSIZE STEP *** C STPM = STP0 + M DEL = V(STPM) IF (DEL*X(XMSAVE) .GT. ZERO) GO TO 130 C *** WE ALREADY TRIED SHRINKING THE STEP, SO QUIT *** IV(FDH) = -2 GO TO 220 C C *** TRY SHRINKING THE STEP *** 130 DEL = NEGPT5 * DEL X(M) = X(XMSAVE) + DEL V(STPM) = DEL IRT = 1 GO TO 999 C C *** SAVE F(X + STP(M)*E(M)) IN H(P,M) *** C 140 PP1O2 = P * (P-1) / 2 HPM = HES + PP1O2 + MM1 V(HPM) = V(F) C C *** START COMPUTING ROW M OF THE FINITE-DIFFERENCE HESSIAN H. *** C HMI = HES + MM1O2 IF (MM1 .EQ. 0) GO TO 160 HPI = HES + PP1O2 DO 150 I = 1, MM1 V(HMI) = V(FX) - (V(F) + V(HPI)) HMI = HMI + 1 HPI = HPI + 1 150 CONTINUE 160 V(HMI) = V(F) - TWO*V(FX) C C *** COMPUTE FUNCTION VALUES NEEDED TO COMPLETE ROW M OF H. *** C I = 1 C 170 IV(SAVEI) = I STPI = STP0 + I V(DELTA) = X(I) X(I) = X(I) + V(STPI) IF (I .EQ. M) X(I) = V(XMSAVE) - V(STPI) IRT = 1 GO TO 999 C 180 X(I) = V(DELTA) IF (IV(TOOBIG) .EQ. 0) GO TO 190 C *** PUNT IN THE EVENT OF AN OVERSIZE STEP *** IV(FDH) = -2 GO TO 220 C C *** FINISH COMPUTING H(M,I) *** C 190 STPI = STP0 + I HMI = HES + MM1O2 + I - 1 STPM = STP0 + M V(HMI) = (V(HMI) + V(F)) / (V(STPI)*V(STPM)) I = I + 1 IF (I .LE. M) GO TO 170 IV(SAVEI) = 0 X(M) = V(XMSAVE) C 200 M = M + 1 IV(MODE) = M IF (M .GT. P) GO TO 210 C C *** PREPARE TO COMPUTE ROW M OF THE FINITE-DIFFERENCE HESSIAN H. C *** COMPUTE M-TH STEP SIZE STP(M), THEN RETURN TO OBTAIN C *** F(X + STP(M)*E(M)), WHERE E(M) = M-TH STD. UNIT VECTOR. C DEL = V(DLTFDC) * MAX(ONE/D(M), ABS(X(M))) IF (X(M) .LT. ZERO) DEL = -DEL V(XMSAVE) = X(M) X(M) = X(M) + DEL STPM = STP0 + M V(STPM) = DEL IRT = 1 GO TO 999 C C *** RESTORE V(F), ETC. *** C 210 IV(FDH) = HES 220 V(F) = V(FX) IRT = 3 IF (KIND .LT. 0) GO TO 999 IV(NFGCAL) = IV(SWITCH) GSAVE1 = IV(W) + P CALL V7CPY(P, G, V(GSAVE1)) GO TO 999 C 999 RETURN C *** LAST LINE OF F7HES FOLLOWS *** END SUBROUTINE G2LRD(DR, IV, L, LH, LIV, LV, ND, N, P, PS, R, RD, 1 RHOI, RHOR, V, W, X, Z) C C *** COMPUTE REGRESSION DIAGNOSTIC FOR RGLG *** C C *** PARAMETERS *** C INTEGER LH, LIV, LV, ND, N, P, PS INTEGER IV(LIV), RHOI(*) REAL DR(ND,P), L(LH), R(N), RD(N), RHOR(*), V(LV), 1 W(P), X(P), Z(P) C C *** CODED BY DAVID M. GAY (SPRING 1986, SUMMER 1991) *** C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL D7TPR, L7ITV, L7IVM, L7SRT, L7SQR, S7LVM, 1 V2AXY, V7CPY, V7SCP REAL D7TPR C C D7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. C L7ITV... MULTIPLIES INVERSE TRANSPOSE OF LOWER TRIANGLE TIMES VECTOR. C L7IVM... APPLY INVERSE OF COMPACT LOWER TRIANG. MATRIX. C L7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. C L7SQR... COMPUTES L*(L**T) FOR LOWER TRIANG. MATRIX L. C S7LVM... MULTIPLIES COMPACTLY STORED SYM. MATRIX TIMES VECTOR. C V2AXY.... ADDS A MULTIPLE OF ONE VECTOR TO ANOTHER. C V7CPY.... COPIES ONE VECTOR TO ANOTHER. C V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C C *** LOCAL VARIABLES *** C LOGICAL USEFLO INTEGER BS1, BSINC, FLO1, FLOINC, H1, HPS1, I, 1 J, J1, K, KI, KI1, KID, L1, LE, LL, LOO1, N1, 2 PMPS, PP1O2, PS1, PX, RDR, XNI, ZAP1, ZAPLEN REAL FRAC, HI, RI, S, T, T1 C C *** CONSTANTS *** C REAL HALF, NEGONE, ONE, ZERO C C C *** IV SUBSCRIPTS *** C INTEGER BS, BSSTR, COVREQ, FDH, FLO, FLOSTR, LOO, NB, NFIX, 1 RDREQ, REGD, XNOTI PARAMETER (BS=85, BSSTR=86, COVREQ=15, FDH=74, FLO=88, FLOSTR=89, 1 LOO=84, NB=87, NFIX=83, RDREQ=57, REGD=67, XNOTI=90) PARAMETER (HALF=0.5E+0, NEGONE=-1.E+0, ONE=1.E+0, ZERO=0.E+0) C C++++++++++++++++++++++++++++++++ BODY +++++++++++++++++++++++++++++++ C I = IV(RDREQ) RDR = MOD(I/2, 3) IF (RDR .EQ. 0) GO TO 999 H1 = IV(FDH) USEFLO = .FALSE. PX = P N1 = N FRAC = ONE XNI = 0 IF (RDR .EQ. 1) GO TO 120 LOO1 = IV(LOO) IF (LOO1 .LE. 0 .OR. LOO1 .GT. 6) THEN IV(REGD) = -1 GO TO 999 ENDIF IF (LOO1 .GT. 3) THEN USEFLO = .TRUE. FLO1 = IV(FLO) FLOINC = IV(FLOSTR) LOO1 = LOO1 - 3 ENDIF XNI = IV(XNOTI) PX = P - IV(NFIX) IF (PX .LT. PS .OR. PX .GT. P) THEN IV(REGD) = -2 GO TO 999 ENDIF IF (LOO1 .EQ. 1) GO TO 120 N1 = IV(NB) IF (N1 .LE. 0 .OR. N1 .GT. N) THEN IV(REGD) = -3 GO TO 999 ENDIF BS1 = IV(BS) BSINC = IV(BSSTR) IF (H1 .LE. 0) GO TO 190 IF (IABS(IV(COVREQ)) .GE. 3) CALL L7SQR(P, V(H1), L) PP1O2 = PX*(PX+1)/2 PS1 = PS + 1 ZAP1 = PS*(PS1)/2 + 1 LE = 0 DO 100 I = 1, N1 IF (USEFLO) THEN FRAC = RHOR(FLO1) FLO1 = FLO1 + FLOINC ENDIF L1 = LE + 1 IF (L1 .GT. N) GO TO 110 LE = LE + RHOI(BS1) IF (LE .GT. N) LE = N BS1 = BS1 + BSINC CALL V7CPY(PP1O2, L, V(H1)) IF (PS .GE. PX) GO TO 50 K = ZAP1 KI = L1 DO 40 J = PS1, P KI = KI + N KI1 = KI DO 10 LL = L1, LE CALL V2AXY(PS, L(K), -FRAC*RD(KI1), DR(1,LL), L(K)) KI1 = KI1 + 1 10 CONTINUE K = K + PS DO 30 J1 = PS1, J KI = KI + N KI1 = KI T = ZERO DO 20 LL = L1, LE T = T + RD(KI1) KI1 = KI1 + 1 20 CONTINUE L(K) = L(K) - FRAC*T K = K + 1 30 CONTINUE 40 CONTINUE 50 DO 70 LL = L1, LE T = -FRAC*RD(LL) K = 1 DO 60 J = 1, PS CALL V2AXY(J, L(K), T*DR(J,LL), DR(1,LL), L(K)) K = K + J 60 CONTINUE 70 CONTINUE CALL L7SRT(1, PX, L, L, J) IF (J .EQ. 0) THEN CALL V7SCP(PX, W, ZERO) DO 90 LL = L1, LE CALL V2AXY(PS, W, R(LL), DR(1,LL), W) IF (PS1 .GT. PX) GO TO 90 K = L1 DO 80 J = PS1, P K = K + N W(J) = W(J) + R(K) 80 CONTINUE 90 CONTINUE CALL L7IVM(PX, W, L, W) CALL L7ITV(PX, W, L, W) CALL S7LVM(PX, Z, V(H1), W) RD(I) = HALF * FRAC * D7TPR(PX, W, Z) IF (XNI .GT. 0) THEN CALL V2AXY(PX, RHOR(XNI), FRAC, W, X) XNI = XNI + PX ENDIF ELSE RD(I) = NEGONE IF (XNI .GT. 0) THEN CALL V7CPY(PX, RHOR(XNI), X) XNI = XNI + PX ENDIF ENDIF 100 CONTINUE 110 IV(REGD) = 1 C *** RESTORE L *** CALL L7SRT(1, P, L, V(H1), J) GO TO 999 C 120 IF (H1 .LE. 0) GO TO 190 IF (IABS(IV(COVREQ)) .GE. 3) CALL L7SQR(P, V(H1), L) IF (PS .GE. PX) GO TO 170 PS1 = PS + 1 PMPS = PX - PS ZAP1 = PS*(PS1)/2 ZAPLEN = PX*(PX+1)/2 - ZAP1 HPS1 = H1 + ZAP1 ZAP1 = ZAP1 + 1 DO 160 I = 1, N IF (USEFLO) THEN FRAC = RHOR(FLO1) FLO1 = FLO1 + FLOINC ENDIF CALL V7CPY(ZAPLEN, L(ZAP1), V(HPS1)) CALL V7SCP(PS, W, ZERO) K = ZAP1 KI = I KID = KI DO 140 J = PS1, PX KI = KI + N CALL V2AXY(PS, L(K), -FRAC*RD(KI), DR(1,I), L(K)) K = K + PS KID = KID + N W(J) = FRAC*R(KID) DO 130 J1 = PS1, J KI = KI + N L(K) = L(K) - FRAC*RD(KI) K = K + 1 130 CONTINUE 140 CONTINUE CALL L7SRT(PS1, PX, L, L, J) IF (J .NE. 0) GO TO 150 CALL V7CPY(PS, Z, DR(1,I)) CALL V7SCP(PMPS, Z(PS1), ZERO) CALL L7IVM(PX, Z, L, Z) HI = D7TPR(PX, Z, Z) CALL L7IVM(PX, W, L, W) RI = FRAC*R(I) C *** FIRST PS ELEMENTS OF W VANISH *** T = D7TPR(PMPS, W(PS1), Z(PS1)) S = FRAC*RD(I) T1 = ONE - S*HI IF (T1 .LE. ZERO) GO TO 150 CALL V2AXY(PX, W, (RI + S*T)/T1, Z, W) CALL L7ITV(PX, W, L, W) CALL S7LVM(PX, Z, V(H1), W) RD(I) = HALF * D7TPR(PX, W, Z) IF (XNI .GT. 0) THEN CALL V2AXY(PX, RHOR(XNI), ONE, W, X) XNI = XNI + PX ENDIF GO TO 160 150 RD(I) = NEGONE IF (XNI .GT. 0) THEN CALL V7CPY(PX, RHOR(XNI), X) XNI = XNI + PX ENDIF 160 CONTINUE C C *** RESTORE L *** C CALL V7CPY(ZAPLEN, L(ZAP1), V(HPS1)) CALL L7SRT(PS1, PX, L, L, J) GO TO 200 C 170 DO 180 I = 1, N IF (USEFLO) THEN FRAC = RHOR(FLO1) FLO1 = FLO1 + FLOINC ENDIF CALL L7IVM(PX, Z, L, DR(1,I)) S = D7TPR(PX, Z, Z) T = ONE - FRAC*RD(I) * S IF (T .LE. ZERO) THEN RD(I) = NEGONE IF (XNI .GT. 0) THEN CALL V7CPY(PX, RHOR(XNI), X) XNI = XNI + PX ENDIF ELSE RD(I) = HALF * FRAC * (R(I)/T)**2 * S IF (XNI .GT. 0) THEN CALL L7ITV(PX, Z, L, Z) CALL V2AXY(PX, RHOR(XNI), FRAC*R(I)/T, Z, X) XNI = XNI + PX ENDIF ENDIF 180 CONTINUE GO TO 200 C 190 CALL V7SCP(N1, RD, NEGONE) 200 IV(REGD) = 1 C 999 RETURN C *** LAST LINE OF G2LRD FOLLOWS *** END SUBROUTINE G7LIT(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 C *** PARAMETER DECLARATIONS *** C INTEGER LIV, LV, P, PS INTEGER IV(LIV) REAL D(P), G(P), V(LV), X(P), Y(P) C C-------------------------- PARAMETER USAGE -------------------------- C C D.... SCALE VECTOR. C IV... INTEGER VALUE ARRAY. C LIV.. LENGTH OF IV. MUST BE AT LEAST 82. 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 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 G7LIT 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 IN V, STARTING AT IV(HC)), AND G7LIT BUILDS AN C APPROXIMATION, S, TO THE SECOND-ORDER TERM. THE CALLER ALSO C PROVIDES THE FUNCTION VALUE, GRADIENT, AND PART OF THE YIELD C VECTOR USED IN UPDATING S. G7LIT DECIDES DYNAMICALLY WHETHER OR C NOT TO USE S WHEN CHOOSING THE NEXT STEP TO TRY... THE HESSIAN C APPROXIMATION USED IS EITHER HC ALONE (GAUSS-NEWTON MODEL) OR C HC + S (AUGMENTED MODEL). C 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, G7LIT 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 G7LIT 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 NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER C (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS C NOT NEEDED). MOREOVER, COMPARED WITH NL2SOL, IV(1) MAY HAVE THE C TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, C AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). THE VALUES IV(D), C IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM NL2SOL (AND C NL2SNO), ARE NOT REFERENCED BY G7LIT OR THE SUBROUTINES IT CALLS. C C WHEN G7LIT IS FIRST CALLED, I.E., WHEN G7LIT IS CALLED WITH C IV(1) = 0 OR 12, V(F), G, AND HC NEED NOT BE INITIALIZED. TO C OBTAIN THESE STARTING VALUES, G7LIT 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 G7LIT WILL MAKE A C 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 G7LIT 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 G7LIT 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 G7LIT 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 G7LIT 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(TOOBIG) TO 1, IN WHICH CASE G7LIT WILL RETURN WITH C IV(1) = 15. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH C SUPPORTED IN PART BY D.O.E. GRANT EX-76-A-01-2295 TO MIT/CCREMS. C C (SEE NL2SOL FOR REFERENCES.) C C+++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C INTEGER DUMMY, DIG1, G01, H1, HC1, I, IPIV1, J, K, L, LMAT1, 1 LSTGST, PP1O2, QTR1, RMAT1, RSTRST, STEP1, STPMOD, S1, 2 TEMP1, TEMP2, W1, X01 REAL E, STTSST, T, T1 C C *** CONSTANTS *** C REAL HALF, NEGONE, ONE, ONEP2, ZERO C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C LOGICAL STOPX REAL D7TPR, L7SVX, L7SVN, RLDST, R7MDC, V2NRM EXTERNAL A7SST, D7TPR, F7HES, G7QTS, ITSUM, L7MST, L7SRT, 1 L7SQR, L7SVX, L7SVN, L7TVM, L7VML, PARCK, RLDST, 2 R7MDC, S7LUP, S7LVM, STOPX, V2AXY, V7CPY, V7SCP, 3 V2NRM C C A7SST.... ASSESSES CANDIDATE STEP. C D7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. C F7HES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE). C G7QTS.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). C ITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. C L7MST... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). C L7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. C L7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L. C L7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C L7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX. C L7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX. C L7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C PARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS. C RLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. C R7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS. C S7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI- C ANGLE OF A SYMMETRIC MATRIX. C STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. C V2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. C V7CPY.... COPIES ONE VECTOR TO ANOTHER. C V7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C V2NRM... RETURNS THE 2-NORM OF A VECTOR. C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, COSMIN, COVMAT, COVREQ, DGNORM, DIG, DSTNRM, F, 1 FDH, FDIF, FUZZ, F0, GTSTEP, H, HC, IERR, INCFAC, INITS, 2 IPIVOT, IRC, KAGQT, KALM, LMAT, LMAX0, LMAXS, MODE, MODEL, 3 MXFCAL, MXITER, NEXTV, NFCALL, NFGCAL, NFCOV, NGCOV, 4 NGCALL, NITER, NVSAVE, PHMXFC, PREDUC, QTR, RADFAC, 5 RADINC, RADIUS, RAD0, RCOND, RDREQ, REGD, RELDX, RESTOR, 6 RMAT, S, SIZE, STEP, STGLIM, STLSTG, STPPAR, SUSED, 7 SWITCH, TOOBIG, TUNER4, TUNER5, VNEED, VSAVE, W, WSCALE, 8 XIRC, X0 C C *** IV SUBSCRIPT VALUES *** C PARAMETER (CNVCOD=55, COVMAT=26, COVREQ=15, DIG=37, FDH=74, H=56, 1 HC=71, IERR=75, INITS=25, IPIVOT=76, IRC=29, KAGQT=33, 2 KALM=34, LMAT=42, MODE=35, MODEL=5, MXFCAL=17, 3 MXITER=18, NEXTV=47, NFCALL=6, NFGCAL=7, NFCOV=52, 4 NGCOV=53, NGCALL=30, NITER=31, QTR=77, RADINC=8, 5 RDREQ=57, REGD=67, RESTOR=9, RMAT=78, S=62, STEP=40, 6 STGLIM=11, STLSTG=41, SUSED=64, SWITCH=12, TOOBIG=2, 7 VNEED=4, VSAVE=60, W=65, XIRC=13, X0=43) C C *** V SUBSCRIPT VALUES *** C 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, RCOND=53, RELDX=17, SIZE=55, STPPAR=5, 4 TUNER4=29, TUNER5=30, WSCALE=56) C C PARAMETER (HALF=0.5E+0, NEGONE=-1.E+0, ONE=1.E+0, ONEP2=1.2E+0, 1 ZERO=0.E+0) C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C I = IV(1) IF (I .EQ. 1) GO TO 40 IF (I .EQ. 2) GO TO 50 C IF (I .EQ. 12 .OR. I .EQ. 13) 1 IV(VNEED) = IV(VNEED) + P*(3*P + 19)/2 + 7 CALL PARCK(1, D, IV, LIV, LV, P, V) I = IV(1) - 2 IF (I .GT. 12) GO TO 999 GO TO (290, 290, 290, 290, 290, 290, 170, 120, 170, 10, 10, 20), I C C *** STORAGE ALLOCATION *** C 10 PP1O2 = P * (P + 1) / 2 IV(S) = IV(LMAT) + PP1O2 IV(X0) = IV(S) + PP1O2 IV(STEP) = IV(X0) + P IV(STLSTG) = IV(STEP) + P IV(DIG) = IV(STLSTG) + P IV(W) = IV(DIG) + P IV(H) = IV(W) + 4*P + 7 IV(NEXTV) = IV(H) + PP1O2 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(STGLIM) = 2 IV(TOOBIG) = 0 IV(CNVCOD) = 0 IV(COVMAT) = 0 IV(NFCOV) = 0 IV(NGCOV) = 0 IV(RADINC) = 0 IV(RESTOR) = 0 IV(FDH) = 0 V(RAD0) = ZERO V(STPPAR) = ZERO V(RADIUS) = V(LMAX0) / (ONE + V(PHMXFC)) C C *** SET INITIAL MODEL AND S MATRIX *** C IV(MODEL) = 1 IF (IV(S) .LT. 0) GO TO 999 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) IV(1) = 1 J = IV(IPIVOT) IF (J .LE. 0) GO TO 999 DO 30 I = 1, P IV(J) = I J = J + 1 30 CONTINUE GO TO 999 C C *** NEW FUNCTION VALUE *** C 40 IF (IV(MODE) .EQ. 0) GO TO 290 IF (IV(MODE) .GT. 0) GO TO 520 C IV(1) = 2 IF (IV(TOOBIG) .EQ. 0) GO TO 999 IV(1) = 63 GO TO 999 C C *** NEW GRADIENT *** C 50 IV(KALM) = -1 IV(KAGQT) = -1 IV(FDH) = 0 IF (IV(MODE) .GT. 0) GO TO 520 C C *** MAKE SURE GRADIENT COULD BE COMPUTED *** C IF (IV(TOOBIG) .EQ. 0) GO TO 60 IV(1) = 65 GO TO 999 60 IF (IV(HC) .LE. 0 .AND. IV(RMAT) .LE. 0) GO TO 610 C C *** COMPUTE D**-1 * GRADIENT *** C DIG1 = IV(DIG) K = DIG1 DO 70 I = 1, P V(K) = G(I) / D(I) K = K + 1 70 CONTINUE V(DGNORM) = V2NRM(P, V(DIG1)) C IF (IV(CNVCOD) .NE. 0) GO TO 510 IF (IV(MODE) .EQ. 0) GO TO 440 IV(MODE) = 0 V(F0) = V(F) IF (IV(INITS) .LE. 2) GO TO 100 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 530 C C *** COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S *** C 80 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 90 CALL V2AXY(PP1O2, V(S1), NEGONE, V(HC1), V(H1)) GO TO 100 90 RMAT1 = IV(RMAT) CALL L7SQR(PS, V(S1), V(RMAT1)) CALL V2AXY(PP1O2, V(S1), NEGONE, V(S1), V(H1)) 100 IV(1) = 2 C C C----------------------------- MAIN LOOP ----------------------------- C C C *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** C 110 CALL ITSUM(D, G, IV, LIV, LV, P, V, X) 120 K = IV(NITER) IF (K .LT. IV(MXITER)) GO TO 130 IV(1) = 10 GO TO 999 130 IV(NITER) = K + 1 C C *** UPDATE RADIUS *** C IF (K .EQ. 0) GO TO 150 STEP1 = IV(STEP) DO 140 I = 1, P V(STEP1) = D(I) * V(STEP1) STEP1 = STEP1 + 1 140 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 150 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 160 IF (.NOT. STOPX(DUMMY)) GO TO 180 IV(1) = 11 GO TO 190 C C *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. C 170 IF (V(F) .GE. V(F0)) GO TO 180 V(RADFAC) = ONE K = IV(NITER) GO TO 130 C 180 IF (IV(NFCALL) .LT. IV(MXFCAL) + IV(NFCOV)) GO TO 200 IV(1) = 9 190 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 430 C C. . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . . C 200 STEP1 = IV(STEP) W1 = IV(W) H1 = IV(H) T1 = ONE IF (IV(MODEL) .EQ. 2) GO TO 210 T1 = ZERO C C *** COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE... C RMAT1 = IV(RMAT) IF (RMAT1 .LE. 0) GO TO 210 QTR1 = IV(QTR) IF (QTR1 .LE. 0) GO TO 210 IPIV1 = IV(IPIVOT) CALL L7MST(D, G, IV(IERR), IV(IPIV1), IV(KALM), P, V(QTR1), 1 V(RMAT1), V(STEP1), V, V(W1)) 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 260 C 210 IF (H1 .GT. 0) GO TO 250 C C *** SET H TO D**-1 * (HC + T1*S) * D**-1. *** C H1 = -H1 IV(H) = H1 IV(FDH) = 0 J = IV(HC) IF (J .GT. 0) GO TO 220 J = H1 RMAT1 = IV(RMAT) CALL L7SQR(P, V(H1), V(RMAT1)) 220 S1 = IV(S) DO 240 I = 1, P T = ONE / D(I) DO 230 K = 1, I V(H1) = T * (V(J) + T1*V(S1)) / D(K) J = J + 1 H1 = H1 + 1 S1 = S1 + 1 230 CONTINUE 240 CONTINUE H1 = IV(H) IV(KAGQT) = -1 C C *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP *** C 250 DIG1 = IV(DIG) LMAT1 = IV(LMAT) CALL G7QTS(D, V(DIG1), V(H1), IV(KAGQT), V(LMAT1), P, V(STEP1), 1 V, V(W1)) IF (IV(KALM) .GT. 0) IV(KALM) = 0 C 260 IF (IV(IRC) .NE. 6) GO TO 270 IF (IV(RESTOR) .NE. 2) GO TO 290 RSTRST = 2 GO TO 300 C C *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** C 270 IV(TOOBIG) = 0 IF (V(DSTNRM) .LE. ZERO) GO TO 290 IF (IV(IRC) .NE. 5) GO TO 280 IF (V(RADFAC) .LE. ONE) GO TO 280 IF (V(PREDUC) .GT. ONEP2 * V(FDIF)) GO TO 280 IF (IV(RESTOR) .NE. 2) GO TO 290 RSTRST = 0 GO TO 300 C C *** COMPUTE F(X0 + STEP) *** C 280 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 999 C C. . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . . C 290 RSTRST = 3 300 X01 = IV(X0) V(RELDX) = RLDST(P, D, X, V(X01)) CALL A7SST(IV, LIV, LV, V) STEP1 = IV(STEP) LSTGST = IV(STLSTG) I = IV(RESTOR) + 1 GO TO (340, 310, 320, 330), I 310 CALL V7CPY(P, X, V(X01)) GO TO 340 320 CALL V7CPY(P, V(LSTGST), V(STEP1)) GO TO 340 330 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 340 IF (IV(SWITCH) .EQ. 0) GO TO 350 IV(H) = -IABS(IV(H)) IV(SUSED) = IV(SUSED) + 2 L = IV(VSAVE) CALL V7CPY(NVSAVE, V, V(L)) 350 L = IV(IRC) - 4 STPMOD = IV(MODEL) IF (L .GT. 0) GO TO (370,380,390,390,390,390,390,390,500,440), 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 360 C C *** SWITCH MODELS *** C IV(MODEL) = 3 - IV(MODEL) IF (-2 .LT. L) GO TO 400 IV(H) = -IABS(IV(H)) IV(SUSED) = IV(SUSED) + 2 L = IV(VSAVE) CALL V7CPY(NVSAVE, V(L), V) GO TO 160 C 360 IF (-3 .LT. L) GO TO 400 C C *** RECOMPUTE STEP WITH NEW RADIUS *** C 370 V(RADIUS) = V(RADFAC) * V(DSTNRM) GO TO 160 C C *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST C 380 V(RADIUS) = V(LMAXS) GO TO 200 C C *** CONVERGENCE OR FALSE CONVERGENCE *** C 390 IV(CNVCOD) = L IF (V(F) .GE. V(F0)) GO TO 510 IF (IV(XIRC) .EQ. 14) GO TO 510 IV(XIRC) = 14 C C. . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . . C 400 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 430 STEP1 = IV(STEP) TEMP1 = IV(STLSTG) TEMP2 = IV(W) C C *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** C HC1 = IV(HC) IF (HC1 .LE. 0) GO TO 410 CALL S7LVM(P, V(TEMP1), V(HC1), V(STEP1)) GO TO 420 410 RMAT1 = IV(RMAT) CALL L7TVM(P, V(TEMP1), V(RMAT1), V(STEP1)) CALL L7VML(P, V(TEMP1), V(RMAT1), V(TEMP1)) C 420 IF (STPMOD .EQ. 1) GO TO 430 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 430 IV(NGCALL) = IV(NGCALL) + 1 G01 = IV(W) CALL V7CPY(P, V(G01), G) IV(1) = 2 IV(TOOBIG) = 0 GO TO 999 C C *** INITIALIZATIONS -- G0 = G - G0, ETC. *** C 440 G01 = IV(W) CALL V2AXY(P, V(G01), NEGONE, V(G01), G) STEP1 = IV(STEP) TEMP1 = IV(STLSTG) TEMP2 = IV(W) IF (IV(IRC) .NE. 3) GO TO 470 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 450 I = 1, P V(K) = (V(K) - V(L)) / D(I) K = K + 1 L = L + 1 450 CONTINUE C C *** DO GRADIENT TESTS *** C IF ( V2NRM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4)) GO TO 460 IF ( D7TPR(P, G, V(STEP1)) 1 .GE. V(GTSTEP) * V(TUNER5)) GO TO 470 460 V(RADFAC) = V(INCFAC) C C *** COMPUTE Y VECTOR NEEDED FOR UPDATING S *** C 470 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 480 CALL S7LVM(PS, V(G01), V(HC1), V(STEP1)) GO TO 490 C 480 RMAT1 = IV(RMAT) CALL L7TVM(PS, V(G01), V(RMAT1), V(STEP1)) CALL L7VML(PS, V(G01), V(RMAT1), V(G01)) C 490 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 110 C C. . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . . C C *** BAD PARAMETERS TO ASSESS *** C 500 IV(1) = 64 GO TO 999 C C C *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** C 510 IF (IV(RDREQ) .EQ. 0) GO TO 600 IF (IV(FDH) .NE. 0) GO TO 600 IF (IV(CNVCOD) .GE. 7) GO TO 600 IF (IV(REGD) .GT. 0) GO TO 600 IF (IV(COVMAT) .GT. 0) GO TO 600 IF (IABS(IV(COVREQ)) .GE. 3) GO TO 560 IF (IV(RESTOR) .EQ. 0) IV(RESTOR) = 2 GO TO 530 C C *** COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE *** C 520 IV(RESTOR) = 0 530 CALL F7HES(D, G, I, IV, LIV, LV, P, V, X) GO TO (540, 550, 580), I 540 IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(1) = 1 GO TO 999 C 550 IV(NGCOV) = IV(NGCOV) + 1 IV(NGCALL) = IV(NGCALL) + 1 IV(NFGCAL) = IV(NFCALL) + IV(NGCOV) IV(1) = 2 GO TO 999 C 560 H1 = IABS(IV(H)) IV(H) = -H1 PP1O2 = P * (P + 1) / 2 RMAT1 = IV(RMAT) IF (RMAT1 .LE. 0) GO TO 570 LMAT1 = IV(LMAT) CALL V7CPY(PP1O2, V(LMAT1), V(RMAT1)) V(RCOND) = ZERO GO TO 590 570 HC1 = IV(HC) IV(FDH) = H1 CALL V7CPY(P*(P+1)/2, V(H1), V(HC1)) C C *** COMPUTE CHOLESKY FACTOR OF FINITE-DIFFERENCE HESSIAN C *** FOR USE IN CALLER*S COVARIANCE CALCULATION... C 580 LMAT1 = IV(LMAT) H1 = IV(FDH) IF (H1 .LE. 0) GO TO 600 IF (IV(CNVCOD) .EQ. 70) GO TO 80 CALL L7SRT(1, P, V(LMAT1), V(H1), I) IV(FDH) = -1 V(RCOND) = ZERO IF (I .NE. 0) GO TO 600 C 590 IV(FDH) = -1 STEP1 = IV(STEP) T = L7SVN(P, V(LMAT1), V(STEP1), V(STEP1)) IF (T .LE. ZERO) GO TO 600 T = T / L7SVX(P, V(LMAT1), V(STEP1), V(STEP1)) IF (T .GT. R7MDC(4)) IV(FDH) = H1 V(RCOND) = T C 600 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 610 IV(1) = 1400 C 999 RETURN C C *** LAST LINE OF G7LIT FOLLOWS *** END SUBROUTINE L7NVR(N, LIN, L) C C *** COMPUTE LIN = L**-1, BOTH N X N LOWER TRIANG. STORED *** C *** COMPACTLY BY ROWS. LIN AND L MAY SHARE THE SAME STORAGE. *** C C *** PARAMETERS *** C INTEGER N REAL L(1), LIN(1) C DIMENSION L(N*(N+1)/2), LIN(N*(N+1)/2) C C *** LOCAL VARIABLES *** C INTEGER I, II, IM1, JJ, J0, J1, K, K0, NP1 REAL ONE, T, ZERO PARAMETER (ONE=1.E+0, ZERO=0.E+0) C C *** BODY *** C NP1 = N + 1 J0 = N*(NP1)/2 DO 30 II = 1, N I = NP1 - II LIN(J0) = ONE/L(J0) IF (I .LE. 1) GO TO 999 J1 = J0 IM1 = I - 1 DO 20 JJ = 1, IM1 T = ZERO J0 = J1 K0 = J1 - JJ DO 10 K = 1, JJ T = T - L(K0)*LIN(J0) J0 = J0 - 1 K0 = K0 + K - I 10 CONTINUE LIN(J0) = T/L(K0) 20 CONTINUE J0 = J0 - 1 30 CONTINUE 999 RETURN C *** LAST LINE OF L7NVR FOLLOWS *** END SUBROUTINE L7TSQ(N, A, L) C C *** SET A TO LOWER TRIANGLE OF (L**T) * L *** C C *** L = N X N LOWER TRIANG. MATRIX STORED ROWWISE. *** C *** A IS ALSO STORED ROWWISE AND MAY SHARE STORAGE WITH L. *** C INTEGER N REAL A(1), L(1) C DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) C INTEGER I, II, IIM1, I1, J, K, M REAL LII, LJ C II = 0 DO 50 I = 1, N I1 = II + 1 II = II + I M = 1 IF (I .EQ. 1) GO TO 30 IIM1 = II - 1 DO 20 J = I1, IIM1 LJ = L(J) DO 10 K = I1, J A(M) = A(M) + LJ*L(K) M = M + 1 10 CONTINUE 20 CONTINUE 30 LII = L(II) DO 40 J = I1, II 40 A(J) = LII * L(J) 50 CONTINUE C 999 RETURN C *** LAST LINE OF L7TSQ FOLLOWS *** END SUBROUTINE N3RDP(IV, LIV, LV, N, P, RD, RHOI, RHOR, V) C C *** PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 *** C INTEGER LIV, LV, N, P INTEGER IV(LIV), RHOI(*) REAL RD(N), RHOR(*), V(LV) C C *** NOTE -- V IS PASSED FOR POSSIBLE USE BY REVISED VERSIONS OF C *** THIS ROUTINE. C INTEGER COV1, I, I1, I2, IEND, II, J, K, K1, NI, PU, PX, PX1, XNI C C *** IV AND V SUBSCRIPTS *** C INTEGER BS, BSSTR, COVMAT, COVPRT, COVREQ, LOO, NB, NEEDHD, NFCOV, 1 NFIX, NGCOV, PRUNIT, RDREQ, REGD, RCOND, STATPR, XNOTI C PARAMETER (BS=85, BSSTR=86, COVMAT=26, COVPRT=14, COVREQ=15, 1 LOO=84, NB=87, NEEDHD=36, NFCOV=52, NFIX=83, NGCOV=53, 2 PRUNIT=21, RDREQ=57, REGD=67, RCOND=53, STATPR=23, 3 XNOTI=90) C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C PU = IV(PRUNIT) IF (PU .EQ. 0) GO TO 999 IF (IV(STATPR) .EQ. 0) GO TO 30 IF (IV(NFCOV) .GT. 0) WRITE(PU,10) IV(NFCOV) 10 FORMAT(/1X,I4,50H EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOST 1ICS.) IF (IV(NGCOV) .GT. 0) WRITE(PU,20) IV(NGCOV) 20 FORMAT(1X,I4,50H EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTI 1CS.) IF (IV(NFCOV) .GT. 0 .OR. IV(NGCOV) .GT. 0) IV(NEEDHD) = 1 C 30 IF (IV(COVPRT) .LE. 0) GO TO 999 COV1 = IV(COVMAT) IF (IV(REGD) .LE. 0 .AND. COV1 .LE. 0) GO TO 70 IV(NEEDHD) = 1 IF (IABS(IV(COVREQ)) .GT. 2) GO TO 50 C WRITE(PU,40) V(RCOND) 40 FORMAT(/53H SQRT(RECIPROCAL CONDITION OF F.D. HESSIAN) = AT MOST, 1 G10.2) GO TO 70 C 50 WRITE(PU,60) V(RCOND) 60 FORMAT(/54H SQRT(RECIPROCAL CONDITION OF (J**T)*RHO"*J) = AT MOST, 1 G10.2) C 70 IF (MOD(IV(COVPRT),2) .EQ. 0) GO TO 210 IV(NEEDHD) = 1 IF (COV1) 80,110,130 80 IF (-1 .EQ. COV1) WRITE(PU,90) 90 FORMAT(/43H ++++++ INDEFINITE COVARIANCE MATRIX ++++++) IF (-2 .EQ. COV1) WRITE(PU,100) 100 FORMAT(/52H ++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++) GO TO 999 C 110 WRITE(PU,120 ) 120 FORMAT(/45H ++++++ COVARIANCE MATRIX NOT COMPUTED ++++++) GO TO 210 C 130 IF (IABS(IV(COVREQ)) .LT. 3) GO TO 150 WRITE(PU,140) 140 FORMAT(/35H COVARIANCE = (J**T * RHO" * J)**-1/) GO TO 170 150 WRITE(PU,160) 160 FORMAT(/56H COVARIANCE = H**-1, WHERE H = FINITE-DIFFERENCE HESSIA 1N/) 170 II = COV1 - 1 DO 180 I = 1, P I1 = II + 1 I2 = II + MIN(I, 5) II = II + I WRITE(PU,190) I, (V(J), J = I1, I2) IF (I .LE. 5) GO TO 180 I2 = I2 + 1 WRITE(PU,200) (V(J), J = I2, II) 180 CONTINUE 190 FORMAT(4H ROW,I3,2X,5G12.3) 200 FORMAT(9X,5G12.3) 210 IF (IV(COVPRT) .LT. 2) GO TO 999 I = IV(REGD) + 4 GO TO (230, 250, 270, 290, 310), I WRITE(PU,220) IV(REGD) 220 FORMAT(/18H BUG... IV(REGD) =,I10) GO TO 999 230 WRITE(PU,240) NB, IV(NB) 240 FORMAT(/17H BAD IV(NB) = IV(,I2,3H) =,I10) GO TO 999 250 WRITE(PU,260) NFIX, IV(NFIX) 260 FORMAT(/19H BAD IV(NFIX) = IV(,I2,3H) =,I10) GO TO 999 270 WRITE(PU,280) LOO, IV(LOO) 280 FORMAT(/18H BAD IV(LOO) = IV(,I2,3H) =,I10) GO TO 999 290 WRITE(PU,300) 300 FORMAT(/42H REGRESSION DIAGNOSTIC VECTOR NOT COMPUTED) GO TO 999 310 IV(NEEDHD) = 1 XNI = 0 I = MOD(IV(RDREQ)/2, 3) + 1 GO TO (999, 330, 320), I 320 XNI = IV(XNOTI) PX = P - IV(NFIX) PX1 = PX - 1 IF (IV(LOO) .GT. 1) GO TO 400 330 WRITE(PU,340) 340 FORMAT (74H REGRESSION DIAGNOSTIC = 0.5 * G(I)**T * H(I)**-1 * H * 1 H(I)**-1 * G(I)...) IF (XNI .LE. 0) GO TO 380 WRITE(PU, 350) 350 FORMAT(29H I RD(I) X(I)) DO 360 I = 1, N WRITE(PU, 370) I, RD(I), (RHOR(J), J = XNI, XNI+PX1) XNI = XNI + PX 360 CONTINUE 370 FORMAT(1X,I5,G13.3,4G15.6/(19X,4G15.6)) GO TO 999 C 380 WRITE(PU,390) RD 390 FORMAT(6G12.3) GO TO 999 C 400 WRITE(PU,410) 410 FORMAT(/77H BLOCK REGRESSION DIAGNOSTIC = 0.5 * G(I)**T * H(I)**-1 1 * H * H(I)**-1 * G(I)) NI = IV(NB) K = IV(BS) K1 = IV(BSSTR) IEND = 0 IF (XNI .GT. 0) GO TO 450 WRITE(PU,420) 420 FORMAT(28H BLOCK FIRST LAST RD(I)) DO 440 I = 1, NI I1 = IEND + 1 IF (I1 .GT. N) GO TO 999 IEND = IEND + RHOI(K) K = K + K1 IF (IEND .GT. N) IEND = N WRITE(PU,430) I, I1, IEND, RD(I) 430 FORMAT(I6,I7,I6,G12.3) 440 CONTINUE GO TO 999 C 450 WRITE(PU, 460) 460 FORMAT(41H BLOCK FIRST LAST RD(I) X(I)) DO 480 I = 1, NI I1 = IEND + 1 IF (I1 .GT. N) GO TO 999 IEND = IEND + RHOI(K) K = K + K1 IF (IEND .GT. N) IEND = N WRITE(PU,470) I, I1, IEND, RD(I), (RHOR(J), J = XNI, XNI+PX1) 470 FORMAT(I6,I7,I6,G12.3,3G15.6/(31X,3G15.6)) XNI = XNI + PX 480 CONTINUE C 999 RETURN C *** LAST LINE OF N3RDP FOLLOWS *** END