SUBROUTINE DRN2G(D, DR, IV, LIV, LV, N, ND, N1, N2, P, R, 1 RD, V, X) c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. c ALL RIGHTS RESERVED. c Based on Government Sponsored Research NAS7-03001. c File: DRN2G.for Ten subrs used by the c David Gay & Linda Kaufman nonlinear LS package. c Needed for versions that do not allow Bounded variables. c DRN2G is called by DNLAFU, DNLAGU, & DRNSG. c C>> 2015-07-09 DRN2G Krogh Introduced TP to avoid divide by 0, C>> 2000-01-07 DRN2G Krogh Moved COV1 = IV(COMAT) up in DN2CVP. C>> 1998-10-29 DRN2G Krogh Moved external statement up for mangle. c>> 1996-07-09 DRN2G Krogh Changes for conversion to C. c>> 1995-01-26 DRN2G Krogh Moved formats up for C conversion. c>> 1994-11-02 DRN2G Krogh Changes to use M77CON c>> 1993-03-10 DRN2G CLL Moved stmt NN = ... to follow IF (IV1 ... c>> 1992-04-27 CLL Comment out unreferenced stmt labels. c>> 1992-04-13 CLL Change from Hollerith to '...' syntax in formats. c>> 1990-06-29 CLL Changes to formats in DN2CVP. c>> 1990-06-12 CLL Revised DRN2G & DG7LIT from DMG 4/19/90 c>> 1990-03-30 CLL JPL c>> 1990-03-14 CLL JPL c>> 1990-06-12 CLL c>> 1990-04-23 CLL (Recent revision by DMG) *** from netlib, Thu Apr 19 11:58:57 EDT 1990 *** c--D replaces "?": ?RN2G,?C7VFN,?D7TPR,?D7UPD,?G7LIT,?ITSUM,?IVSET c--& ?L7VML,?N2CVP,?N2LRD,?Q7APL,?Q7RAD,?V2NRM,?V7CPY,?V7SCP, ?N2G c--& ?A7SST,?F7HES,?G7QTS,?L7MST,?L7SQR,?L7SRT,?L7SVN,?L7SVX,?L7TVM c--& ?PARCK,?R7MDC,?RLDST,?S7LUP,?S7LVM,?V2AXY,?L7ITV,?L7IVM,?O7PRD c--& ?L7NVR,?L7TSQ,?V7SCL,?N2RDP,?NLAFU,?NLAGU,?RNSG C C *** REVISED ITERATION DRIVER FOR NL2SOL (VERSION 2.3) *** C INTEGER LIV, LV, N, ND, N1, N2, P INTEGER IV(LIV) DOUBLE PRECISION D(P), DR(ND,P), R(ND), RD(ND), V(LV), X(P) 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 + 82. C LV....... LENGTH OF V... LV MUST BE AT LEAST 105 + P*(2*P+16). C N........ TOTAL NUMBER OF RESIDUALS. C ND....... MAX. NO. OF RESIDUALS PASSED ON ONE CALL. C N1....... LOWEST ROW INDEX FOR RESIDUALS SUPPLIED THIS TIME. C N2....... HIGHEST ROW INDEX FOR RESIDUALS SUPPLIED THIS TIME. C P........ NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. C R........ RESIDUALS. C RD....... RD(I) = SQRT(G(I)**T * H(I)**-1 * G(I)) ON OUTPUT WHEN C IV(RDREQ) IS NONZERO. DRN2G 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. IF ND .GE. N, THEN RD IS ALSO USED AS C TEMPORARY STORAGE. C V........ FLOATING-POINT VALUES ARRAY. C X........ PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, C OUTPUT = BEST VALUE FOUND). C C *** DISCUSSION *** C C NOTE... NL2SOL AND NL2ITR (MENTIONED BELOW) ARE DESCRIBED IN C ACM TRANS. MATH. SOFTWARE, VOL. 7, PP. 369-383 (AN ADAPTIVE C NONLINEAR LEAST-SQUARES ALGORITHM, BY J.E. DENNIS, D.M. GAY, C AND R.E. WELSCH). C C THIS ROUTINE CARRIES OUT ITERATIONS FOR SOLVING NONLINEAR C LEAST SQUARES PROBLEMS. WHEN ND = N, IT IS SIMILAR TO NL2ITR C (WITH J = DR), EXCEPT THAT R(X) AND DR(X) NEED NOT BE INITIALIZED C WHEN DRN2G IS CALLED WITH IV(1) = 0 OR 12. DRN2G ALSO ALLOWS C R AND DR TO BE SUPPLIED ROW-WISE -- JUST SET ND = 1 AND CALL C DRN2G ONCE FOR EACH ROW WHEN PROVIDING RESIDUALS AND JACOBIANS. C ANOTHER NEW FEATURE IS THAT CALLING DRN2G WITH IV(1) = 13 C CAUSES STORAGE ALLOCATION ONLY TO BE PERFORMED -- ON RETURN, SUCH C COMPONENTS AS IV(G) (THE FIRST SUBSCRIPT IN G OF THE GRADIENT) C AND IV(S) (THE FIRST SUBSCRIPT IN V OF THE S LOWER TRIANGLE OF C THE S MATRIX) WILL HAVE BEEN SET (UNLESS LIV OR LV IS TOO SMALL), C AND IV(1) WILL HAVE BEEN SET TO 14. CALLING DRN2G WITH IV(1) = 14 C CAUSES EXECUTION OF THE ALGORITHM TO BEGIN UNDER THE ASSUMPTION C THAT STORAGE HAS BEEN ALLOCATED. C C *** SUPPLYING R AND DR *** C C DRN2G USES IV AND V IN THE SAME WAY AS NL2SOL, WITH A SMALL C NUMBER OF OBVIOUS CHANGES. ONE DIFFERENCE BETWEEN DRN2G AND C NL2ITR IS THAT INITIAL FUNCTION AND GRADIENT INFORMATION NEED NOT C BE SUPPLIED IN THE VERY FIRST CALL ON DRN2G, THE ONE WITH C IV(1) = 0 OR 12. ANOTHER DIFFERENCE IS THAT DRN2G RETURNS WITH C IV(1) = -2 WHEN IT WANTS ANOTHER LOOK AT THE OLD JACOBIAN MATRIX C AND THE CURRENT RESIDUAL -- THE ONE CORRESPONDING TO X AND C IV(NFGCAL). IT THEN RETURNS WITH IV(1) = -3 WHEN IT WANTS TO SEE C BOTH THE NEW RESIDUAL AND THE NEW JACOBIAN MATRIX AT ONCE. NOTE C THAT IV(NFGCAL) = IV(7) CONTAINS THE VALUE THAT IV(NFCALL) = IV(6) C HAD WHEN THE CURRENT RESIDUAL WAS EVALUATED. ALSO NOTE THAT THE C VALUE OF X CORRESPONDING TO THE OLD JACOBIAN MATRIX IS STORED IN C V, STARTING AT V(IV(X0)) = V(IV(43)). C ANOTHER NEW RETURN... DRN2G IV(1) = -1 WHEN IT WANTS BOTH THE C RESIDUAL AND THE JACOBIAN TO BE EVALUATED AT X. C A NEW RESIDUAL VECTOR MUST BE SUPPLIED WHEN DRN2G RETURNS WITH C IV(1) = 1 OR -1. THIS TAKES THE FORM OF VALUES OF R(I,X) PASSED C IN R(I-N1+1), I = N1(1)N2. YOU MAY PASS ALL THESE VALUES AT ONCE C (I.E., N1 = 1 AND N2 = N) OR IN PIECES BY MAKING SEVERAL CALLS ON C DRN2G. EACH TIME DRN2G RETURNS WITH IV(1) = 1, N1 WILL HAVE C BEEN SET TO THE INDEX OF THE NEXT RESIDUAL THAT DRN2G EXPECTS TO C SEE, AND N2 WILL BE SET TO THE INDEX OF THE HIGHEST RESIDUAL THAT C COULD BE GIVEN ON THE NEXT CALL, I.E., N2 = N1 + ND - 1. (THUS C WHEN DRN2G FIRST RETURNS WITH IV(1) = 1 FOR A NEW X, IT WILL C HAVE SET N1 TO 1 AND N2 TO MIN(ND,N).) THE CALLER MAY PROVIDE C FEWER THAN N2-N1+1 RESIDUALS ON THE NEXT CALL BY SETTING N2 TO C A SMALLER VALUE. DRN2G ASSUMES IT HAS SEEN ALL THE RESIDUALS C FOR THE CURRENT X WHEN IT IS CALLED WITH N2 .GE. N. C EXAMPLE... SUPPOSE N = 80 AND THAT R IS TO BE PASSED IN 8 C BLOCKS OF SIZE 10. THE FOLLOWING CODE WOULD DO THE JOB. C C N = 80 C ND = 10 C ... C DO 10 K = 1, 8 C *** COMPUTE R(I,X) FOR I = 10*K-9 TO 10*K *** C *** AND STORE THEM IN R(1),...,R(10) *** C CALL DRN2G(..., R, ...) C 10 CONTINUE C C THE SITUATION IS SIMILAR WHEN GRADIENT INFORMATION IS C REQUIRED, I.E., WHEN DRN2G RETURNS WITH IV(1) = 2, -1, OR -2. C NOTE THAT DRN2G OVERWRITES R, BUT THAT IN THE SPECIAL CASE OF C N1 = 1 AND N2 = N ON PREVIOUS CALLS, DRN2G NEVER RETURNS WITH C IV(1) = -2. IT SHOULD BE CLEAR THAT THE PARTIAL DERIVATIVE OF C R(I,X) WITH RESPECT TO X(L) IS TO BE STORED IN DR(I-N1+1,L), C L = 1(1)P, I = N1(1)N2. IT IS ESSENTIAL THAT R(I) AND DR(I,L) C ALL CORRESPOND TO THE SAME RESIDUALS WHEN IV(1) = -1 OR -2. C C *** COVARIANCE MATRIX *** C C IV(RDREQ) = IV(57) TELLS WHETHER TO COMPUTE A COVARIANCE C MATRIX AND/OR REGRESSION DIAGNOSTICS... 0 MEANS NEITHER, C 1 MEANS COVARIANCE MATRIX ONLY, 2 MEANS REG. DIAGNOSTICS ONLY, C 3 MEANS BOTH. AS WITH NL2SOL, IV(COVREQ) = IV(15) TELLS WHAT C HESSIAN APPROXIMATION TO USE IN THIS COMPUTING. C C *** REGRESSION DIAGNOSTICS *** C C SEE THE COMMENTS IN SUBROUTINE DN2G. C C *** GENERAL *** C C CODED BY DAVID M. GAY. C C ++++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++ C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL DC7VFN,DIVSET, DD7TPR,DD7UPD,DG7LIT,DITSUM,DL7VML, 1 DN2CVP, DN2LRD, DQ7APL,DQ7RAD,DV7CPY, DV7SCP, DV2NRM DOUBLE PRECISION DD7TPR, DV2NRM c ------------------------------------------------------------------ C DC7VFN... FINISHES COVARIANCE COMPUTATION. C DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C DD7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. C DD7UPD... UPDATES SCALE VECTOR D. C DG7LIT.... PERFORMS BASIC MINIMIZATION ALGORITHM. C DITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X. C DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C DN2CVP... PRINTS COVARIANCE MATRIX. C DN2LRD... COMPUTES REGRESSION DIAGNOSTICS. C DQ7APL... APPLIES QR TRANSFORMATIONS STORED BY DQ7RAD. C DQ7RAD.... ADDS A NEW BLOCK OF ROWS TO QR DECOMPOSITION. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C C *** LOCAL VARIABLES *** C INTEGER G1, GI, I, IV1, IVMODE, JTOL1, K, L, LH, NN, QTR1, 1 RMAT1, YI, Y1 DOUBLE PRECISION T C DOUBLE PRECISION HALF, ZERO C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, COVMAT, COVREQ, DINIT, DTYPE, DTINIT, D0INIT, F, 1 FDH, G, H, IPIVOT, IVNEED, JCN, JTOL, LMAT, MODE, 2 NEXTIV, NEXTV, NF0, NF00, NF1, NFCALL, NFCOV, NFGCAL, 3 NGCALL, NGCOV, QTR, RDREQ, REGD, RESTOR, RLIMIT, RMAT, 4 TOOBIG, VNEED, Y C C *** IV SUBSCRIPT VALUES *** C PARAMETER (CNVCOD=55, COVMAT=26, COVREQ=15, DTYPE=16, FDH=74, 1 G=28, H=56, IPIVOT=76, IVNEED=3, JCN=66, JTOL=59, 2 LMAT=42, MODE=35, NEXTIV=46, NEXTV=47, NFCALL=6, 3 NFCOV=52, NF0=68, NF00=81, NF1=69, NFGCAL=7, NGCALL=30, 4 NGCOV=53, QTR=77, RESTOR=9, RMAT=78, RDREQ=57, REGD=67, 5 TOOBIG=2, VNEED=4, Y=48) C C *** V SUBSCRIPT VALUES *** PARAMETER (DINIT=38, DTINIT=39, D0INIT=40, F=10, RLIMIT=46) PARAMETER (HALF=0.5D+0, ZERO=0.D+0) C C ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C LH = P * (P+1) / 2 IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) IV1 = IV(1) IF (IV1 .GT. 2) GO TO 10 NN = N2 - N1 + 1 IV(RESTOR) = 0 I = IV1 + 4 IF (IV(TOOBIG) .EQ. 0) GO TO (150, 130, 150, 120, 120, 150), I IF (I .NE. 5) IV(1) = 2 GO TO 40 C C *** FRESH START OR RESTART -- CHECK INPUT INTEGERS *** C 10 IF (ND .LE. 0) GO TO 210 IF (P .LE. 0) GO TO 210 IF (N .LE. 0) GO TO 210 IF (IV1 .EQ. 14) GO TO 30 IF (IV1 .GT. 16) GO TO 300 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 20 CALL DG7LIT(D, X, IV, LIV, LV, P, P, 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 IV(JCN) = IV(G) + P IV(RMAT) = IV(JCN) + P IV(QTR) = IV(RMAT) + LH IV(JTOL) = IV(QTR) + P IV(NEXTV) = IV(JTOL) + 2*P IF (IV1 .EQ. 13) GO TO 999 C 30 JTOL1 = IV(JTOL) IF (V(DINIT) .GE. ZERO) CALL DV7SCP(P, D, V(DINIT)) IF (V(DTINIT) .GT. ZERO) CALL DV7SCP(P, V(JTOL1), V(DTINIT)) I = JTOL1 + P IF (V(D0INIT) .GT. ZERO) CALL DV7SCP(P, V(I), V(D0INIT)) IV(NF0) = 0 IV(NF1) = 0 IF (ND .GE. N) GO TO 40 C C *** SPECIAL CASE HANDLING OF FIRST FUNCTION AND GRADIENT EVALUATION C *** -- ASK FOR BOTH RESIDUAL AND JACOBIAN AT ONCE C G1 = IV(G) Y1 = IV(Y) CALL DG7LIT(D, V(G1), IV, LIV, LV, P, P, V, X, V(Y1)) IF (IV(1) .NE. 1) GO TO 220 V(F) = ZERO CALL DV7SCP(P, V(G1), ZERO) IV(1) = -1 QTR1 = IV(QTR) CALL DV7SCP(P, V(QTR1), ZERO) IV(REGD) = 0 RMAT1 = IV(RMAT) GO TO 100 C 40 G1 = IV(G) Y1 = IV(Y) CALL DG7LIT(D, V(G1), IV, LIV, LV, P, P, V, X, V(Y1)) IF (IV(1) - 2) 50, 60, 220 C 50 V(F) = ZERO IF (IV(NF1) .EQ. 0) GO TO 260 IF (IV(RESTOR) .NE. 2) GO TO 260 IV(NF0) = IV(NF1) CALL DV7CPY(N, RD, R) IV(REGD) = 0 GO TO 260 C 60 CALL DV7SCP(P, V(G1), ZERO) IF (IV(MODE) .GT. 0) GO TO 230 RMAT1 = IV(RMAT) QTR1 = IV(QTR) CALL DV7SCP(P, V(QTR1), ZERO) IV(REGD) = 0 IF (ND .LT. N) GO TO 90 IF (N1 .NE. 1) GO TO 90 IF (IV(MODE) .LT. 0) GO TO 100 IF (IV(NF1) .EQ. IV(NFGCAL)) GO TO 70 IF (IV(NF0) .NE. IV(NFGCAL)) GO TO 90 CALL DV7CPY(N, R, RD) GO TO 80 70 CALL DV7CPY(N, RD, R) 80 CALL DQ7APL(ND, N, P, DR, RD, 0) CALL DL7VML(P, V(Y1), V(RMAT1), RD) GO TO 110 C 90 IV(1) = -2 IF (IV(MODE) .LT. 0) IV(1) = -1 100 CALL DV7SCP(P, V(Y1), ZERO) 110 CALL DV7SCP(LH, V(RMAT1), ZERO) GO TO 260 C C *** COMPUTE F(X) *** C 120 T = DV2NRM(NN, R) IF (T .GT. V(RLIMIT)) GO TO 200 V(F) = V(F) + HALF * T**2 IF (N2 .LT. N) GO TO 270 IF (N1 .EQ. 1) IV(NF1) = IV(NFCALL) GO TO 40 C C *** COMPUTE Y *** C 130 Y1 = IV(Y) YI = Y1 DO 140 L = 1, P V(YI) = V(YI) + DD7TPR(NN, DR(1,L), R) YI = YI + 1 140 CONTINUE IF (N2 .LT. N) GO TO 270 IV(1) = 2 IF (N1 .GT. 1) IV(1) = -3 GO TO 260 C C *** COMPUTE GRADIENT INFORMATION *** C 150 IF (IV(MODE) .GT. P) GO TO 240 G1 = IV(G) IVMODE = IV(MODE) IF (IVMODE .LT. 0) GO TO 170 IF (IVMODE .EQ. 0) GO TO 180 IV(1) = 2 C C *** COMPUTE GRADIENT ONLY (FOR USE IN COVARIANCE COMPUTATION) *** C GI = G1 DO 160 L = 1, P V(GI) = V(GI) + DD7TPR(NN, R, DR(1,L)) GI = GI + 1 160 CONTINUE GO TO 190 C C *** COMPUTE INITIAL FUNCTION VALUE WHEN ND .LT. N *** C 170 IF (N .LE. ND) GO TO 180 T = DV2NRM(NN, R) IF (T .GT. V(RLIMIT)) GO TO 200 V(F) = V(F) + HALF * T**2 C C *** UPDATE D IF DESIRED *** C 180 IF (IV(DTYPE) .GT. 0) 1 CALL DD7UPD(D, DR, IV, LIV, LV, N, ND, NN, N2, P, V) C C *** COMPUTE RMAT AND QTR *** C QTR1 = IV(QTR) RMAT1 = IV(RMAT) CALL DQ7RAD(NN, ND, P, V(QTR1), .TRUE., V(RMAT1), DR, R) IV(NF1) = 0 C 190 IF (N2 .LT. N) GO TO 270 IF (IVMODE .GT. 0) GO TO 40 IV(NF00) = IV(NFGCAL) C C *** COMPUTE G FROM RMAT AND QTR *** C CALL DL7VML(P, V(G1), V(RMAT1), V(QTR1)) IV(1) = 2 IF (IVMODE .EQ. 0) GO TO 40 IF (N .LE. ND) GO TO 40 C C *** FINISH SPECIAL CASE HANDLING OF FIRST FUNCTION AND GRADIENT C Y1 = IV(Y) IV(1) = 1 CALL DG7LIT(D, V(G1), IV, LIV, LV, P, P, V, X, V(Y1)) IF (IV(1) .NE. 2) GO TO 220 GO TO 40 C C *** MISC. DETAILS *** C C *** X IS OUT OF RANGE (OVERSIZE STEP) *** C 200 IV(TOOBIG) = 1 GO TO 40 C C *** BAD N, ND, OR P *** C 210 IV(1) = 66 GO TO 300 C C *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** C 220 IF (IV(COVMAT) .NE. 0) GO TO 290 IF (IV(REGD) .NE. 0) GO TO 290 C C *** SEE IF CHOLESKY FACTOR OF HESSIAN IS AVAILABLE *** C K = IV(FDH) IF (K .LE. 0) GO TO 280 IF (IV(RDREQ) .LE. 0) GO TO 290 C C *** COMPUTE REGRESSION DIAGNOSTICS AND DEFAULT COVARIANCE IF C DESIRED *** C I = 0 IF (mod(IV(RDREQ),4) .GE. 2) I = 1 IF (mod(IV(RDREQ),2) .EQ. 1 .AND. abs(IV(COVREQ)) .LE. 1) I = I+2 IF (I .EQ. 0) GO TO 250 IV(MODE) = P + I IV(NGCALL) = IV(NGCALL) + 1 IV(NGCOV) = IV(NGCOV) + 1 IV(CNVCOD) = IV(1) IF (I .LT. 2) GO TO 230 L = abs(IV(H)) CALL DV7SCP(LH, V(L), ZERO) 230 IV(NFCOV) = IV(NFCOV) + 1 IV(NFCALL) = IV(NFCALL) + 1 IV(NFGCAL) = IV(NFCALL) IV(1) = -1 GO TO 260 C 240 L = IV(LMAT) CALL DN2LRD(DR, IV, V(L), LH, LIV, LV, ND, NN, P, R, RD, V) IF (N2 .LT. N) GO TO 270 IF (N1 .GT. 1) GO TO 250 C C *** ENSURE WE CAN RESTART -- AND MAKE RETURN STATE OF DR C *** INDEPENDENT OF WHETHER REGRESSION DIAGNOSTICS ARE COMPUTED. C *** USE STEP VECTOR (ALLOCATED BY DG7LIT) FOR SCRATCH. C RMAT1 = IV(RMAT) CALL DV7SCP(LH, V(RMAT1), ZERO) CALL DQ7RAD(NN, ND, P, R, .FALSE., V(RMAT1), DR, R) IV(NF1) = 0 C C *** FINISH COMPUTING COVARIANCE *** C 250 L = IV(LMAT) CALL DC7VFN(IV, V(L), LH, LIV, LV, N, P, V) GO TO 290 C C *** RETURN FOR MORE FUNCTION OR GRADIENT INFORMATION *** C 260 N2 = 0 270 N1 = N2 + 1 N2 = N2 + ND IF (N2 .GT. N) N2 = N GO TO 999 C C *** COME HERE FOR INDEFINITE FINITE-DIFFERENCE HESSIAN *** C 280 IV(COVMAT) = K IV(REGD) = K C C *** PRINT SUMMARY OF FINAL ITERATION AND OTHER REQUESTED ITEMS *** C 290 G1 = IV(G) 300 CALL DITSUM(D, V(G1), IV, LIV, LV, P, V, X) IF (IV(1) .LE. 6 .AND. IV(RDREQ) .GT. 0) 1 CALL DN2CVP(IV, LIV, LV, P, V) C 999 RETURN C *** LAST LINE OF DRN2G FOLLOWS *** END c ================================================================== SUBROUTINE DG7LIT(D, GG, IV, LIV, LV, P, PS, V, X, YY) c>> 1990-06-12 CLL @ JPL c>> 1990-04-23 CLL (Recent revision by DMG) *** from netlib, Mon Apr 23 20:37:24 EDT 1990 *** c>> 1990-02-20 CLL @ JPL 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) DOUBLE PRECISION D(P), GG(P), V(LV), X(P), YY(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 GG... 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 YY... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE). C C *** DISCUSSION *** C C DG7LIT 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 DG7LIT 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. DG7LIT 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,DG7LIT 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 YY, NAMELY THE SUM OVER I OF C RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING DG7LIT WITH IV(1) = 2 AND C IV(MODE) = 0 (WHERE MODE = 38). GG THEN CONTANS THE OTHER PART, C SO THAT THE DESIRED YIELD VECTOR IS GG - YY. IF PS .LT. P, THEN C THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF C GRAD(R(I,X)), STEP, AND YY. 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 DG7LIT OR THE SUBROUTINES IT CALLS. C C WHEN DG7LIT IS FIRST CALLED, I.E., WHEN DG7LIT IS CALLED WITH C IV(1) = 0 OR 12, V(F), GG, AND HC NEED NOT BE INITIALIZED. TO C OBTAIN THESE STARTING VALUES,DG7LIT 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 YY MUST ALSO BE SUPPLIED. (NOTE THAT YY IS USED FOR SCRATCH -- c ITS 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 DG7LIT WILL MAKE A C NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE. C WHEN IV(MODE) IS POSITIVE, YY 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 DG7LIT 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 DG7LIT 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 GG, HC, AND YY THE NEXT TIME DG7LIT RETURNS WITH C IV(1) = 2. SEE MLPIT FOR AN EXAMPLE OF THIS. C IV(1) = 2 MEANS THE CALLER SHOULD SET GG 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 DG7LIT 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 GG OR HC C CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET C IV(TOOBIG) TO 1, IN WHICH CASE DG7LIT 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 ------------------------------------------------------------------ c References to the function STOPX have been commented out of this c subroutine. If one wishes to be able to terminate this package c gracefully using a keybord "Break" key, one can provide a STOPX c function that returns .true. if the Break key has been pressed c since the last call to STOPX, and otherwise returns .false., and c then uncomment the references to STOPX in this subr. c -- CLL 6/12/90 C ++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++ C C *** LOCAL VARIABLES *** C c integer DUMMY INTEGER DIG1, G01, H1, HC1, I, IPIV1, J, K, L, LMAT1, 1 LSTGST, PP1O2, QTR1, RMAT1, RSTRST, STEP1, STPMOD, S1, 2 TEMP1, TEMP2, W1, X01 DOUBLE PRECISION E, STTSST, T, T1, TP C C *** CONSTANTS *** C DOUBLE PRECISION HALF, NEGONE, ONE, ONEP2, ZERO C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C c external STOPX c LOGICAL STOPX EXTERNAL DA7SST, DD7TPR,DF7HES,DG7QTS,DITSUM, DL7MST,DL7SRT, 1 DL7SQR, DL7SVX, DL7SVN, DL7TVM,DL7VML,DPARCK, DRLDST, 2 DR7MDC, DS7LUP, DS7LVM, DV2AXY,DV7CPY, DV7SCP, 3 DV2NRM DOUBLE PRECISION DD7TPR, DL7SVX, DL7SVN, DRLDST, DR7MDC, DV2NRM c ------------------------------------------------------------------ C DA7SST.... ASSESSES CANDIDATE STEP. C DD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. C DF7HES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE). C DG7QTS.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). C DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. C DL7MST... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). C DL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX. C DL7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L. C DL7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C DL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX. C DL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX. C DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C DPARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS. C DRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. C DR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS. C DS7LUP... 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 Call to STOPX commented out. -- CLL 6/12/90 C DV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. C DV7CPY.... COPIES ONE VECTOR TO ANOTHER. C DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. C DV2NRM... 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) PARAMETER (HALF=0.5D+0, NEGONE=-1.D+0, ONE=1.D+0, ONEP2=1.2D+0, 1 ZERO=0.D+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 DPARCK(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 DV7SCP(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) = GG(I) / D(I) K = K + 1 70 CONTINUE V(DGNORM) = DV2NRM(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 DV2AXY(PP1O2, V(S1), NEGONE, V(HC1), V(H1)) GO TO 100 90 RMAT1 = IV(RMAT) CALL DL7SQR(PS, V(S1), V(RMAT1)) CALL DV2AXY(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 DITSUM(D, GG, 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) * DV2NRM(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) = -abs(IV(H)) IV(SUSED) = IV(MODEL) C C *** COPY X TO X0 *** C CALL DV7CPY(P, V(X01), X) C C *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** C 160 continue c if (STOPX(DUMMY)) then c IV(1) = 11 c GO TO 190 c else go to 180 c endif 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 c 190 continue 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 DL7MST(D, GG, 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) = -abs(H1) C *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO C *** MARK INVALID THE INFORMATION DG7QTS 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 DL7SQR(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 DG7QTS(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 DV2AXY(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) = DRLDST(P, D, X, V(X01)) CALL DA7SST(IV, LIV, LV, V) STEP1 = IV(STEP) LSTGST = IV(STLSTG) I = IV(RESTOR) + 1 GO TO (340, 310, 320, 330), I 310 CALL DV7CPY(P, X, V(X01)) GO TO 340 320 CALL DV7CPY(P, V(LSTGST), V(STEP1)) GO TO 340 330 CALL DV7CPY(P, V(STEP1), V(LSTGST)) CALL DV2AXY(P, X, ONE, V(STEP1), V(X01)) V(RELDX) = DRLDST(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) = -abs(IV(H)) IV(SUSED) = IV(SUSED) + 2 L = IV(VSAVE) CALL DV7CPY(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 DS7LVM(PS, YY, V(S1), V(STEP1)) STTSST = HALF * DD7TPR(PS, V(STEP1), YY) 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) = -abs(IV(H)) IV(SUSED) = IV(SUSED) + 2 L = IV(VSAVE) CALL DV7CPY(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 DS7LVM(P, V(TEMP1), V(HC1), V(STEP1)) GO TO 420 410 RMAT1 = IV(RMAT) CALL DL7TVM(P, V(TEMP1), V(RMAT1), V(STEP1)) CALL DL7VML(P, V(TEMP1), V(RMAT1), V(TEMP1)) C 420 IF (STPMOD .EQ. 1) GO TO 430 S1 = IV(S) CALL DS7LVM(PS, V(TEMP2), V(S1), V(STEP1)) CALL DV2AXY(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 DV7CPY(P, V(G01), GG) IV(1) = 2 IV(TOOBIG) = 0 GO TO 999 C C *** INITIALIZATIONS -- G0 = GG - G0, ETC. *** C 440 G01 = IV(W) CALL DV2AXY(P, V(G01), NEGONE, V(G01), GG) 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 (DV2NRM(P, V(TEMP1)) .LE. V(DGNORM) * V(TUNER4)) GO TO 460 IF (DD7TPR(P, GG, V(STEP1)) 1 .GE. V(GTSTEP) * V(TUNER5)) GO TO 470 460 V(RADFAC) = V(INCFAC) C C *** COMPUTE YY VECTOR NEEDED FOR UPDATING S *** C 470 CALL DV2AXY(PS, YY, NEGONE, YY, GG) C C *** DETERMINE SIZING FACTOR V(SIZE) *** C C *** SET TEMP1 = S * STEP *** S1 = IV(S) CALL DS7LVM(PS, V(TEMP1), V(S1), V(STEP1)) C T1 = abs(DD7TPR(PS, V(STEP1), V(TEMP1))) T = abs(DD7TPR(PS, V(STEP1), YY)) 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 DS7LVM(PS, V(G01), V(HC1), V(STEP1)) GO TO 490 C 480 RMAT1 = IV(RMAT) CALL DL7TVM(PS, V(G01), V(RMAT1), V(STEP1)) CALL DL7VML(PS, V(G01), V(RMAT1), V(G01)) C 490 CALL DV2AXY(PS, V(G01), ONE, YY, V(G01)) C C *** UPDATE S *** C CALL DS7LUP(V(S1), V(COSMIN), PS, V(SIZE), V(STEP1), V(TEMP1), 1 V(TEMP2), V(G01), V(WSCALE), YY) 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 (abs(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 DF7HES(D, GG, 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 = abs(IV(H)) IV(H) = -H1 PP1O2 = P * (P + 1) / 2 RMAT1 = IV(RMAT) IF (RMAT1 .LE. 0) GO TO 570 LMAT1 = IV(LMAT) CALL DV7CPY(PP1O2, V(LMAT1), V(RMAT1)) V(RCOND) = ZERO GO TO 590 570 HC1 = IV(HC) IV(FDH) = H1 CALL DV7CPY(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 DL7SRT(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 = DL7SVN(P, V(LMAT1), V(STEP1), V(STEP1)) IF (T .LE. ZERO) GO TO 600 TP = DL7SVX(P, V(LMAT1), V(STEP1), V(STEP1)) IF (TP .NE. ZERO) then T = T / TP IF (T .GT. DR7MDC(4)) IV(FDH) = H1 V(RCOND) = T END IF 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 DG7LIT FOLLOWS *** END c ================================================================== SUBROUTINE DN2LRD(DR, IV, L, LH, LIV, LV, ND, NN, P, R, RD, V) C C *** COMPUTE REGRESSION DIAGNOSTIC AND DEFAULT COVARIANCE MATRIX FOR C DRN2G *** C C *** PARAMETERS *** C INTEGER LH, LIV, LV, ND, NN, P INTEGER IV(LIV) DOUBLE PRECISION DR(ND,P), L(LH), R(NN), RD(NN), V(LV) C C *** CODED BY DAVID M. GAY (WINTER 1982, FALL 1983) *** C C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C EXTERNAL DD7TPR, DL7ITV, DL7IVM,DO7PRD, DV7SCP DOUBLE PRECISION DD7TPR C C *** LOCAL VARIABLES *** C INTEGER COV, I, J, M, STEP1 DOUBLE PRECISION A, S, T C C *** CONSTANTS *** C DOUBLE PRECISION NEGONE, ONE, ONEV(1), ZERO C C *** IV SUBSCRIPTS *** C INTEGER D, H, MODE, RDREQ, STEP PARAMETER (D=27, H=56, MODE=35, RDREQ=57, STEP=40) PARAMETER (NEGONE=-1.D+0, ONE=1.D+0, ZERO=0.D+0) DATA ONEV(1)/1.D+0/ C C +++++++++++++++++++++++++++++++ BODY +++++++++++++++++++++++++++++++ C STEP1 = IV(STEP) I = IV(RDREQ) IF (I .LE. 0) GO TO 999 IF (MOD(I,4) .LT. 2) GO TO 30 CALL DV7SCP(NN, RD, NEGONE) DO 20 I = 1, NN A = R(I)**2 M = STEP1 DO 10 J = 1, P V(M) = DR(I,J) M = M + 1 10 CONTINUE CALL DL7IVM(P, V(STEP1), L, V(STEP1)) S = DD7TPR(P, V(STEP1), V(STEP1)) T = ONE - S IF (T .LE. ZERO) GO TO 20 A = A * S / T RD(I) = sqrt(A) 20 CONTINUE C 30 IF (IV(MODE) - P .LT. 2) GO TO 999 C C *** COMPUTE DEFAULT COVARIANCE MATRIX *** C COV = abs(IV(H)) DO 50 I = 1, NN M = STEP1 DO 40 J = 1, P V(M) = DR(I,J) M = M + 1 40 CONTINUE CALL DL7IVM(P, V(STEP1), L, V(STEP1)) CALL DL7ITV(P, V(STEP1), L, V(STEP1)) CALL DO7PRD(1, LH, P, V(COV), ONEV, V(STEP1), V(STEP1)) 50 CONTINUE C 999 RETURN C *** LAST CARD OF DN2LRD FOLLOWS *** END SUBROUTINE DC7VFN(IV, L, LH, LIV, LV, N, P, V) C C *** FINISH COVARIANCE COMPUTATION FOR DRN2G, DRNSG *** C INTEGER LH, LIV, LV, N, P INTEGER IV(LIV) DOUBLE PRECISION L(LH), V(LV) C EXTERNAL DL7NVR, DL7TSQ, DV7SCL C C *** LOCAL VARIABLES *** C INTEGER COV, I DOUBLE PRECISION HALF C C *** SUBSCRIPTS FOR IV AND V *** C INTEGER CNVCOD, COVMAT, F, FDH, H, MODE, RDREQ, REGD C PARAMETER (CNVCOD=55, COVMAT=26, F=10, FDH=74, H=56, MODE=35, 1 RDREQ=57, REGD=67) DATA HALF/0.5D+0/ C C *** BODY *** C IV(1) = IV(CNVCOD) I = IV(MODE) - P IV(MODE) = 0 IV(CNVCOD) = 0 IF (IV(FDH) .LE. 0) GO TO 999 IF ((I-2)**2 .EQ. 1) IV(REGD) = 1 IF (MOD(IV(RDREQ),2) .NE. 1) GO TO 999 C C *** FINISH COMPUTING COVARIANCE MATRIX = INVERSE OF F.D. HESSIAN. C COV = abs(IV(H)) IV(FDH) = 0 C IF (IV(COVMAT) .NE. 0) GO TO 999 IF (I .GE. 2) GO TO 10 CALL DL7NVR(P, V(COV), L) CALL DL7TSQ(P, V(COV), V(COV)) C 10 CALL DV7SCL(LH, V(COV), V(F)/(HALF * dble(max(1,N-P))), V(COV)) IV(COVMAT) = COV C 999 RETURN C *** LAST LINE OF DC7VFN FOLLOWS *** END SUBROUTINE DF7HES(D, GG, 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 DF7HES USES GRADIENT DIFFERENCES, C *** OTHERWISE FUNCTION DIFFERENCES. STORAGE IN V IS AS IN DG7LIT. 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) DOUBLE PRECISION D(P), GG(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 DOUBLE PRECISION DEL, HALF, NEGPT5, ONE, TWO, ZERO C C *** EXTERNAL SUBROUTINES *** C EXTERNAL DV7CPY C C DV7CPY.... 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 PARAMETER (HALF=0.5D+0, NEGPT5=-0.5D+0, ONE=1.D+0, TWO=2.D+0, 1 ZERO=0.D+0) 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) = -abs(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 DF7HES. SET GSAVE = G, TAKE FIRST STEP *** CALL DV7CPY(P, V(GSAVE1), GG) 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 GG = (GG - GSAVE)/DEL *** C DO 50 I = 1, P GG(I) = (GG(I) - V(GSAVE1)) / DEL GSAVE1 = GSAVE1 + 1 50 CONTINUE C C *** ADD GG 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) + GG(I)) FOR I = 1 TO M-1 *** C MM1 = M - 1 DO 60 I = 1, MM1 V(K) = HALF * (V(K) + GG(I)) K = K + 1 60 CONTINUE C C *** ADD H(I,M) = GG(I) FOR I = M TO P *** C 70 L = L + 1 DO 80 I = M, P V(L) = GG(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 GG 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 DF7HES. *** 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 DV7CPY(P, GG, V(GSAVE1)) GO TO 999 C 999 RETURN C *** LAST CARD OF DF7HES FOLLOWS *** END SUBROUTINE DN2CVP(IV, LIV, LV, P, V) C C *** PRINT COVARIANCE MATRIX FOR DRN2G *** C 6/27/90 CLL changed 'SCALE' to 'VARFAC' in output labels. c ------------------------------------------------------------------ c%% long int j, k; INTEGER J INTEGER LIV, LV, P INTEGER IV(LIV) DOUBLE PRECISION V(LV) C C *** LOCAL VARIABLES *** C INTEGER COV1, I, II, I1, PU DOUBLE PRECISION T C C *** IV SUBSCRIPTS *** C INTEGER COVMAT, COVPRT, COVREQ, NEEDHD, NFCOV, NGCOV, PRUNIT, 1 RCOND, REGD, STATPR C PARAMETER (COVMAT=26, COVPRT=14, COVREQ=15, NEEDHD=36, NFCOV=52, 1 NGCOV=53, PRUNIT=21, RCOND=53, REGD=67, STATPR=23) C *** BODY *** C 10 FORMAT(/1X,I4, * ' EXTRA FUNC. EVALS FOR COVARIANCE AND DIAGNOSTICS.') 20 FORMAT(1X,I4, * ' EXTRA GRAD. EVALS FOR COVARIANCE AND DIAGNOSTICS.') 40 FORMAT(/' RECIPROCAL CONDITION OF F.D. HESSIAN = AT MOST',g10.2) 60 FORMAT(/' RECIPROCAL CONDITION OF (J**T)*J = AT LEAST',g10.2) 90 FORMAT(/' ++++++ INDEFINITE COVARIANCE MATRIX ++++++') 100 FORMAT(/' ++++++ OVERSIZE STEPS IN COMPUTING COVARIANCE +++++') 120 FORMAT(/' ++++++ COVARIANCE MATRIX NOT COMPUTED ++++++') c c++(~.C.) Default UNITNO='(PU,' c++(.C.) Default UNITNO='(*,' c++ Replace "(PU," = UNITNO c IF (IV(1) .GT. 8) GO TO 999 PU = IV(PRUNIT) IF (PU .EQ. 0) GO TO 999 COV1 = IV(COVMAT) IF (-2 .EQ. COV1) WRITE(PU,100) IF (IV(STATPR) .EQ. 0) GO TO 30 IF (IV(NFCOV) .GT. 0) WRITE(PU,10) IV(NFCOV) IF (IV(NGCOV) .GT. 0) WRITE(PU,20) IV(NGCOV) C 30 IF (IV(COVPRT) .LE. 0) GO TO 999 IF (IV(REGD) .LE. 0 .AND. COV1 .LE. 0) GO TO 70 IV(NEEDHD) = 1 T = V(RCOND)**2 IF (abs(IV(COVREQ)) .GT. 2) GO TO 50 C WRITE(PU,40) T GO TO 70 C 50 WRITE(PU,60) T C 70 IF (MOD(IV(COVPRT),2) .EQ. 0) GO TO 999 IV(NEEDHD) = 1 IF (COV1) 80,110,130 80 IF (-1 .EQ. COV1) WRITE(PU,90) GO TO 999 C 110 WRITE(PU,120) GO TO 999 C 130 I = abs(IV(COVREQ)) IF (I .LE. 1) WRITE(PU,'(/ * '' COVARIANCE = VARFAC * H**-1 * (J**T * J) * H**-1''/ * '' WHERE H = F.D. HESSIAN''/1x)') IF (I .EQ. 2) WRITE(PU,'(1x/'' COVARIANCE = VARFAC * H**-1,'', * '' WHERE H = FINITE-DIFFERENCE HESSIAN''/1x)') IF (I.GT.2) WRITE(PU,'(/'' COVARIANCE = VARFAC * J**T * J''/1x)') II = COV1 - 1 DO 170 I = 1, P I1 = II + 1 II = II + I C%% printf( " ROW%3ld ", i ); C%% for (j = i1; j <= ii; j+=5){ C%% for (k = j; k <= (j <= ii - 5 ? j+4 : ii); k++) C%% printf( "%12.3g", V[k] ); C%% printf( "\n"); C%% if (j <= ii - 5) printf( " ");} WRITE(PU,'('' ROW'',I3,2X,5g12.3/(9X,5g12.3))')I,(V(J),J=I1,II) 170 CONTINUE C 999 RETURN C *** LAST CARD OF DN2CVP FOLLOWS *** END SUBROUTINE DN2RDP(IV, LIV, N, RD) C C *** PRINT REGRESSION DIAGNOSTICS FOR MLPSL AND NL2S1 *** C c ------------------------------------------------------------------ c++ Code for .C. is inactive c%% long int j,k; c++ End INTEGER LIV, N INTEGER IV(LIV) DOUBLE PRECISION RD(N) INTEGER PU C C *** IV SUBSCRIPTS *** C INTEGER COVPRT, NEEDHD, PRUNIT, RDREQ, REGD C C DATA COVPRT/14/, NEEDHD/36/, PRUNIT/21/, RDREQ/57/, REGD/67/ PARAMETER (COVPRT=14, NEEDHD=36, PRUNIT=21, RDREQ=57, REGD=67) C C ++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C PU = IV(PRUNIT) IF (PU .EQ. 0) GO TO 999 IF (IV(COVPRT) .LT. 2) GO TO 999 IF (IV(REGD) .LE. 0) GO TO 999 IV(NEEDHD) = 1 WRITE(PU,'('' REGRESSION DIAGNOSTIC = SQRT(G(I)**T * H(I)**-1 *'', * ''G(I))...''/)') C%% for(j=0; j < n; j+=6){ C%% for (k = j; k < (j < n - 6 ? j+6 : n); k++) C%% printf( "%12.3g", rd[k] ); C%% printf( "\n" );} WRITE(PU,'(6g12.3)') RD C 999 RETURN C *** LAST CARD OF DN2RDP FOLLOWS *** END SUBROUTINE DO7PRD(L, LS, PP, SS, WW, YY, Z) C C *** FOR I = 1..L, SET SS = SS + WW(I)*YY(.,I)*(Z(.,I)**T), I.E., C *** ADD WW(I) TIMES THE OUTER PRODUCT OF Y(.,I) AND Z(.,I). C c ------------------------------------------------------------------ INTEGER L, LS, PP DOUBLE PRECISION SS(LS), WW(L), YY(PP,L), Z(PP,L) C DIMENSION SS(PP*(PP+1)/2) C INTEGER I, J, K, M DOUBLE PRECISION WK, YI, ZERO parameter(ZERO = 0.0D0) C DO 30 K = 1, L WK = WW(K) IF (WK .EQ. ZERO) GO TO 30 M = 1 DO 20 I = 1, PP YI = WK * YY(I,K) DO 10 J = 1, I SS(M) = SS(M) + YI*Z(J,K) M = M + 1 10 CONTINUE 20 CONTINUE 30 CONTINUE C RETURN C *** LAST CARD OF DO7PRD FOLLOWS *** END SUBROUTINE DL7NVR(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 c ------------------------------------------------------------------ INTEGER N DOUBLE PRECISION L(*), LIN(*) 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 DOUBLE PRECISION ONE, T, ZERO PARAMETER (ONE=1.D+0, ZERO=0.D+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 CARD OF DL7NVR FOLLOWS *** END SUBROUTINE DL7TSQ(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 c ------------------------------------------------------------------ INTEGER N DOUBLE PRECISION A(*), L(*) C DIMENSION A(N*(N+1)/2), L(N*(N+1)/2) C INTEGER I, II, IIM1, I1, J, K, M DOUBLE PRECISION 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 RETURN C *** LAST CARD OF DL7TSQ FOLLOWS *** END