SUBROUTINE DRN2G(D, DR, IV, LIV, LV, N, ND, N1, N2, P, R, 1 RD, V, X) 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 *** INTRINSIC FUNCTIONS *** C/+ INTEGER IABS, MOD C/ C *** EXTERNAL FUNCTIONS AND SUBROUTINES *** C DOUBLE PRECISION DD7TPR, DV2NRM EXTERNAL DC7VFN,DIVSET, DD7TPR,DD7UPD,DG7LIT,DITSUM,DL7VML, 1 DN2CVP, DN2LRD, DQ7APL,DQ7RAD,DV7CPY, DV7SCP, 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 C/6 C DATA CNVCOD/55/, COVMAT/26/, COVREQ/15/, DTYPE/16/, FDH/74/, C 1 G/28/, H/56/, IPIVOT/76/, IVNEED/3/, JCN/66/, JTOL/59/, C 2 LMAT/42/, MODE/35/, NEXTIV/46/, NEXTV/47/, NFCALL/6/, C 3 NFCOV/52/, NF0/68/, NF00/81/, NF1/69/, NFGCAL/7/, NGCALL/30/, C 4 NGCOV/53/, QTR/77/, RESTOR/9/, RMAT/78/, RDREQ/57/, REGD/67/, C 5 TOOBIG/2/, VNEED/4/, Y/48/ C/7 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 C *** V SUBSCRIPT VALUES *** C C/6 C DATA DINIT/38/, DTINIT/39/, D0INIT/40/, F/10/, RLIMIT/46/ C/7 PARAMETER (DINIT=38, DTINIT=39, D0INIT=40, F=10, RLIMIT=46) C/ C/6 C DATA HALF/0.5D+0/, ZERO/0.D+0/ C/7 PARAMETER (HALF=0.5D+0, ZERO=0.D+0) C/ 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. IABS(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 = IABS(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