1 SUBROUTINE DNLAGU(N, P, X, DCALCR, DCALCJ, IV, LIV, LV, V) 2c Copyright (c) 1996 California Institute of Technology, Pasadena, CA. 3c ALL RIGHTS RESERVED. 4c Based on Government Sponsored Research NAS7-03001. 5c>> 1996-04-27 DNLAGU Krogh Changes to get desired C prototypes. 6c>> 1994-10-20 DNLAGU Krogh Changes to use M77CON 7c>> 1990-07-02 DNLAGU C. L. Lawson, JPL 8c>> 1990-01-31 C. L. Lawson, JPL 9C 10C *** VERSION OF NL2SOL THAT CALLS DRN2G *** 11C 12C *** PARAMETERS *** 13C 14 INTEGER N, P, LIV, LV 15 INTEGER IV(LIV) 16 DOUBLE PRECISION X(P), V(LV) 17C/ 18 EXTERNAL DCALCR, DCALCJ 19C 20C *** PARAMETER USAGE *** 21C 22C N....... TOTAL NUMBER OF RESIDUALS. 23C P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. 24C X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, 25C OUTPUT = BEST VALUE FOUND). 26C DCALCR... SUBROUTINE FOR COMPUTING RESIDUAL VECTOR. 27C DCALCJ... SUBROUTINE FOR COMPUTING JACOBIAN MATRIX = MATRIX OF FIRST 28C PARTIALS OF THE RESIDUAL VECTOR. 29C IV...... INTEGER VALUES ARRAY. 30C LIV..... LENGTH OF IV (SEE DISCUSSION BELOW). 31C LV...... LENGTH OF V (SEE DISCUSSION BELOW). 32C V....... FLOATING-POINT VALUES ARRAY. 33C 34C 35C *** DISCUSSION *** 36C 37C NOTE... NL2SOL (MENTIONED BELOW) IS A CODE FOR SOLVING 38C NONLINEAR LEAST-SQUARES PROBLEMS. IT IS DESCRIBED IN 39C ACM TRANS. MATH. SOFTWARE, VOL. 9, PP. 369-383 (AN ADAPTIVE 40C NONLINEAR LEAST-SQUARES ALGORITHM, BY J.E. DENNIS, D.M. GAY, 41C AND R.E. WELSCH). 42C 43C LIV GIVES THE LENGTH OF IV. IT MUST BE AT LEAST 82+P. IF NOT, 44C THEN DNLAGU RETURNS WITH IV(1) = 15. WHEN DNLAGU RETURNS, THE 45C MINIMUM ACCEPTABLE VALUE OF LIV IS STORED IN IV(LASTIV) = IV(44), 46C (PROVIDED THAT LIV .GE. 44). 47C 48C LV GIVES THE LENGTH OF V. THE MINIMUM VALUE FOR LV IS 49C LV0 = 105 + P*(N + 2*P + 17) + 2*N. IF LV IS SMALLER THAN THIS, 50C THEN DNLAGU RETURNS WITH IV(1) = 16. WHEN DNLAGU RETURNS, THE 51C MINIMUM ACCEPTABLE VALUE OF LV IS STORED IN IV(LASTV) = IV(45) 52C (PROVIDED LIV .GE. 45). 53C 54C RETURN CODES AND CONVERGENCE TOLERANCES ARE THE SAME AS FOR 55C NL2SOL, WITH SOME SMALL EXTENSIONS... IV(1) = 15 MEANS LIV WAS 56C TOO SMALL. IV(1) = 16 MEANS LV WAS TOO SMALL. 57C 58C THERE ARE TWO NEW V INPUT COMPONENTS... V(LMAXS) = V(36) AND 59C V(SCTOL) = V(37) SERVE WHERE V(LMAX0) AND V(RFCTOL) FORMERLY DID 60C IN THE SINGULAR CONVERGENCE TEST -- SEE THE NL2SOL DOCUMENTATION. 61C 62C *** DEFAULT VALUES *** 63C 64C DEFAULT VALUES ARE PROVIDED BY SUBROUTINE DIVSET, RATHER THAN 65C DFAULT. THE CALLING SEQUENCE IS... 66C CALL DIVSET(1, IV, LIV, LV, V) 67C THE FIRST PARAMETER IS AN INTEGER 1. IF LIV AND LV ARE LARGE 68C ENOUGH FOR DIVSET, THEN DIVSET SETS IV(1) TO 12. OTHERWISE IT 69C SETS IV(1) TO 15 OR 16. CALLING DNLAGU WITH IV(1) = 0 CAUSES ALL 70C DEFAULT VALUES TO BE USED FOR THE INPUT COMPONENTS OF IV AND V. 71C IF YOU FIRST CALL DIVSET, THEN SET IV(1) TO 13 AND CALL DNLAGU, 72C THEN STORAGE ALLOCATION ONLY WILL BE PERFORMED. IN PARTICULAR, 73C IV(D) = IV(27), IV(J) = IV(70), AND IV(R) = IV(61) WILL BE SET 74C TO THE FIRST SUBSCRIPT IN V OF THE SCALE VECTOR, THE JACOBIAN 75C MATRIX, AND THE RESIDUAL VECTOR RESPECTIVELY, PROVIDED LIV AND LV 76C ARE LARGE ENOUGH. IF SO, THEN DNLAGU RETURNS WITH IV(1) = 14. 77C WHEN CALLED WITH IV(1) = 14, DNLAGU ASSUMES THAT STORAGE HAS 78C BEEN ALLOCATED, AND IT BEGINS THE MINIMIZATION ALGORITHM. 79C 80C *** SCALE VECTOR *** 81C 82C ONE DIFFERENCE WITH NL2SOL IS THAT THE SCALE VECTOR D IS 83C STORED IN V, STARTING AT A DIFFERENT SUBSCRIPT. THE STARTING 84C SUBSCRIPT VALUE IS STILL STORED IN IV(D) = IV(27). THE 85C DISCUSSION OF DEFAULT VALUES ABOVE TELLS HOW TO HAVE IV(D) SET 86C BEFORE THE ALGORITHM IS STARTED. 87C 88C *** REGRESSION DIAGNOSTICS *** 89C 90C IF IV(RDREQ) SO DICTATES, THEN ESTIMATES ARE COMPUTED OF THE 91C INFLUENCE EACH RESIDUAL COMPONENT HAS ON THE FINAL PARAMETER 92C ESTIMATE X. THE GENERAL IDEA IS THAT ONE MAY WISH TO EXAMINE 93C RESIDUAL COMPONENTS (AND THE DATA BEHIND THEM) FOR WHICH THE 94C INFLUENCE ESTIMATE IS SIGNIFICANTLY LARGER THAN MOST OF THE OTHER 95C INFLUENCE ESTIMATES. THESE ESTIMATES, HEREAFTER CALLED 96C REGRESSION DIAGNOSTICS, ARE ONLY COMPUTED IF IV(RDREQ) = 2 OR 3. 97C IN THIS CASE, FOR I = 1(1)N, 98C SQRT( G(I)**T * H(I)**-1 * G(I) ) 99C IS COMPUTED AND STORED IN V, STARTING AT V(IV(REGD)), WHERE 100C RDREQ = 57 AND REGD = 67. HERE G(I) STANDS FOR THE GRADIENT 101C RESULTING WHEN THE I-TH OBSERVATION IS DELETED AND H(I) STANDS 102C FOR AN APPROXIMATION TO THE CORRESPONDING HESSIAN AT X, THE SOLU- 103C TION CORRESPONDING TO ALL OBSERVATIONS. (THIS APPROXIMATION IS 104C OBTAINED BY SUBTRACTING THE FIRST-ORDER CONTRIBUTION OF THE I-TH 105C OBSERVATION TO THE HESSIAN FROM A FINITE-DIFFERENCE HESSIAN 106C APPROXIMATION. IF H IS INDEFINITE, THEN IV(REGD) IS SET TO -1. 107C IF H(I) IS INDEFINITE, THEN -1 IS RETURNED AS THE DIAGNOSTIC FOR 108C OBSERVATION I. IF NO DIAGNOSTICS ARE COMPUTED, PERHAPS BECAUSE 109C OF A FAILURE TO CONVERGE, THEN IV(REGD) = 0 IS RETURNED.) 110C PRINTING OF THE REGRESSION DIAGNOSTICS IS CONTROLLED BY 111C IV(COVPRT) = IV(14)... IF IV(COVPRT) = 3, THEN BOTH THE 112C COVARIANCE MATRIX AND THE REGRESSION DIAGNOSTICS ARE PRINTED. 113C IV(COVPRT) = 2 CAUSES ONLY THE REGRESSION DIAGNOSTICS TO BE 114C PRINTED, IV(COVPRT) = 1 CAUSES ONLY THE COVARIANCE MATRIX TO BE 115C PRINTED, AND IV(COVPRT) = 0 CAUSES NEITHER TO BE PRINTED. 116C 117C RDREQ = 57 AND REGD = 67. 118C 119C *** GENERAL *** 120C 121C CODED BY DAVID M. GAY. 122C 123C+++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++ 124C 125C *** EXTERNAL SUBROUTINES *** 126C 127 EXTERNAL DIVSET, DRN2G, DN2RDP 128c 129c--D replaces "?": ?NLAGU, ?RN2G, ?IVSET, ?N2RDP, ?CALCR, ?CALCJ 130c 131C DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. 132C DRN2G ... CARRIES OUT OPTIMIZATION ITERATIONS. 133C DN2RDP... PRINTS REGRESSION DIAGNOSTICS. 134C 135C *** NO INTRINSIC FUNCTIONS *** 136C 137C *** LOCAL VARIABLES *** 138C 139 INTEGER D1, DR1, IV1, N1, N2, NF, R1, RD1 140C 141C *** IV COMPONENTS *** 142C 143 INTEGER D, J, NEXTV, NFCALL, NFGCAL, R, REGD, REGD0, TOOBIG, VNEED 144 PARAMETER (D=27, J=70, NEXTV=47, NFCALL=6, NFGCAL=7, R=61, 145 1 REGD=67, REGD0=82, TOOBIG=2, VNEED=4) 146c -------------------------------- BODY ------------------------------ 147C 148 IF (IV(1) .EQ. 0) CALL DIVSET(1, IV, LIV, LV, V) 149 IV1 = IV(1) 150 IF (IV1 .EQ. 14) GO TO 10 151 IF (IV1 .GT. 2 .AND. IV1 .LT. 12) GO TO 10 152 IF (IV1 .EQ. 12) IV(1) = 13 153 IF (IV(1) .EQ. 13) IV(VNEED) = IV(VNEED) + P + N*(P+2) 154 CALL DRN2G(X, V, IV, LIV, LV, N, N, N1, N2, P, V, V, V, X) 155 IF (IV(1) .NE. 14) GO TO 999 156C 157C *** STORAGE ALLOCATION *** 158C 159 IV(D) = IV(NEXTV) 160 IV(R) = IV(D) + P 161 IV(REGD0) = IV(R) + N 162 IV(J) = IV(REGD0) + N 163 IV(NEXTV) = IV(J) + N*P 164 IF (IV1 .EQ. 13) GO TO 999 165C 166 10 D1 = IV(D) 167 DR1 = IV(J) 168 R1 = IV(R) 169 RD1 = IV(REGD0) 170C 171 20 CALL DRN2G(V(D1), V(DR1), IV, LIV, LV, N, N, N1, N2, P, V(R1), 172 1 V(RD1), V, X) 173 IF (IV(1)-2) 30, 50, 60 174C 175C *** NEW FUNCTION VALUE (R VALUE) NEEDED *** 176C 177 30 NF = IV(NFCALL) 178C%% (*dcalcr)( n, p, x, &nf, &V[r1] ); 179 CALL DCALCR(N, P, X, NF, V(R1)) 180 IF (NF .GT. 0) GO TO 40 181 IV(TOOBIG) = 1 182 GO TO 20 183 40 IF (IV(1) .GT. 0) GO TO 20 184C 185C *** COMPUTE DR = GRADIENT OF R COMPONENTS *** 186C 187 50 CONTINUE 188C%% (*dcalcj)( n, p, x, &Iv[NFGCAL], &V[dr1] ); 189 CALL DCALCJ(N, P, X, IV(NFGCAL), V(DR1)) 190 IF (IV(NFGCAL) .EQ. 0) IV(TOOBIG) = 1 191 GO TO 20 192C 193C *** INDICATE WHETHER THE REGRESSION DIAGNOSTIC ARRAY WAS COMPUTED 194C *** AND PRINT IT IF SO REQUESTED... 195C 196 60 IF (IV(REGD) .GT. 0) IV(REGD) = RD1 197 CALL DN2RDP(IV, LIV, N, V(RD1)) 198C 199 999 RETURN 200C 201C *** LAST CARD OF DNLAGU FOLLOWS *** 202 END 203