SUBROUTINE RN2G(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) REAL 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. RN2G 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 RN2G IS CALLED WITH IV(1) = 0 OR 12. RN2G ALSO ALLOWS C R AND DR TO BE SUPPLIED ROW-WISE -- JUST SET ND = 1 AND CALL C RN2G ONCE FOR EACH ROW WHEN PROVIDING RESIDUALS AND JACOBIANS. C ANOTHER NEW FEATURE IS THAT CALLING RN2G 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 RN2G 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 RN2G USES IV AND V IN THE SAME WAY AS NL2SOL, WITH A SMALL C NUMBER OF OBVIOUS CHANGES. ONE DIFFERENCE BETWEEN RN2G AND C NL2ITR IS THAT INITIAL FUNCTION AND GRADIENT INFORMATION NEED NOT C BE SUPPLIED IN THE VERY FIRST CALL ON RN2G, THE ONE WITH C IV(1) = 0 OR 12. ANOTHER DIFFERENCE IS THAT RN2G 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... RN2G 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 RN2G 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 RN2G. EACH TIME RN2G RETURNS WITH IV(1) = 1, N1 WILL HAVE C BEEN SET TO THE INDEX OF THE NEXT RESIDUAL THAT RN2G 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 RN2G 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. RN2G 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 RN2G(..., R, ...) C 10 CONTINUE C C THE SITUATION IS SIMILAR WHEN GRADIENT INFORMATION IS C REQUIRED, I.E., WHEN RN2G RETURNS WITH IV(1) = 2, -1, OR -2. C NOTE THAT RN2G OVERWRITES R, BUT THAT IN THE SPECIAL CASE OF C N1 = 1 AND N2 = N ON PREVIOUS CALLS, RN2G 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 N2G. 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 REAL D7TPR, V2NRM EXTERNAL C7VFN, IVSET, D7TPR, D7UPD, G7LIT, ITSUM, L7VML, 1 N2CVP, N2LRD, Q7APL, Q7RAD, V7CPY, V7SCP, V2NRM C C C7VFN... FINISHES COVARIANCE COMPUTATION. C IVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. C D7TPR... COMPUTES INNER PRODUCT OF TWO VECTORS. C D7UPD... UPDATES SCALE VECTOR D. C G7LIT.... PERFORMS BASIC MINIMIZATION ALGORITHM. C ITSUM.... PRINTS ITERATION SUMMARY, INFO ABOUT INITIAL AND FINAL X. C L7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. C N2CVP... PRINTS COVARIANCE MATRIX. C N2LRD... COMPUTES REGRESSION DIAGNOSTICS. C Q7APL... APPLIES QR TRANSFORMATIONS STORED BY Q7RAD. C Q7RAD.... ADDS A NEW BLOCK OF ROWS TO QR DECOMPOSITION. C V7CPY.... COPIES ONE VECTOR TO ANOTHER. C V7SCP... 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 REAL T C REAL 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.5E+0/, ZERO/0.E+0/ C/7 PARAMETER (HALF=0.5E+0, ZERO=0.E+0) C/ C C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++ C LH = P * (P+1) / 2 IF (IV(1) .EQ. 0) CALL IVSET(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 G7LIT(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 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 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 G7LIT(D, V(G1), IV, LIV, LV, P, P, V, X, V(Y1)) IF (IV(1) .NE. 1) GO TO 220 V(F) = ZERO CALL V7SCP(P, V(G1), ZERO) IV(1) = -1 QTR1 = IV(QTR) CALL V7SCP(P, V(QTR1), ZERO) IV(REGD) = 0 RMAT1 = IV(RMAT) GO TO 100 C 40 G1 = IV(G) Y1 = IV(Y) CALL G7LIT(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 V7CPY(N, RD, R) IV(REGD) = 0 GO TO 260 C 60 CALL V7SCP(P, V(G1), ZERO) IF (IV(MODE) .GT. 0) GO TO 230 RMAT1 = IV(RMAT) QTR1 = IV(QTR) CALL V7SCP(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 V7CPY(N, R, RD) GO TO 80 70 CALL V7CPY(N, RD, R) 80 CALL Q7APL(ND, N, P, DR, RD, 0) CALL L7VML(P, V(Y1), V(RMAT1), RD) GO TO 110 C 90 IV(1) = -2 IF (IV(MODE) .LT. 0) IV(1) = -1 100 CALL V7SCP(P, V(Y1), ZERO) 110 CALL V7SCP(LH, V(RMAT1), ZERO) GO TO 260 C C *** COMPUTE F(X) *** C 120 T = V2NRM(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) + D7TPR(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) + D7TPR(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 = V2NRM(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 D7UPD(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 Q7RAD(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 L7VML(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 G7LIT(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 V7SCP(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 N2LRD(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 G7LIT) FOR SCRATCH. C RMAT1 = IV(RMAT) CALL V7SCP(LH, V(RMAT1), ZERO) CALL Q7RAD(NN, ND, P, R, .FALSE., V(RMAT1), DR, R) IV(NF1) = 0 C C *** FINISH COMPUTING COVARIANCE *** C 250 L = IV(LMAT) CALL C7VFN(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 ITSUM(D, V(G1), IV, LIV, LV, P, V, X) IF (IV(1) .LE. 6 .AND. IV(RDREQ) .GT. 0) 1 CALL N2CVP(IV, LIV, LV, P, V) C 999 RETURN C *** LAST LINE OF RN2G FOLLOWS *** END