1 /* dn2g.f -- translated by f2c (version 19970211).
2 You must link the resulting object file with the libraries:
3 -lf2c -lm (in that order)
4 */
5
6 #include "f2c.h"
7
8 /* Table of constant values */
9
10 static integer c__1 = 1;
11
dn2g_(n,p,x,calcr,calcj,iv,liv,lv,v,ui,ur,uf)12 /* Subroutine */ int dn2g_(n, p, x, calcr, calcj, iv, liv, lv, v, ui, ur, uf)
13 integer *n, *p;
14 doublereal *x;
15 /* Subroutine */ int (*calcr) (), (*calcj) ();
16 integer *iv, *liv, *lv;
17 doublereal *v;
18 integer *ui;
19 doublereal *ur;
20 /* Subroutine */ int (*uf) ();
21 {
22 /* System generated locals */
23 integer i__1;
24
25 /* Local variables */
26 extern /* Subroutine */ int drn2g_();
27 static integer d1, n1, n2, r1;
28 extern /* Subroutine */ int dn2rdp_();
29 static integer nf;
30 extern /* Subroutine */ int divset_();
31 static integer dr1, rd1, iv1;
32
33
34 /* *** VERSION OF NL2SOL THAT CALLS DRN2G *** */
35
36 /* *** PARAMETERS *** */
37
38 /* /6 */
39 /* INTEGER IV(LIV), UI(1) */
40 /* DOUBLE PRECISION X(P), V(LV), UR(1) */
41 /* /7 */
42 /* / */
43
44 /* *** PARAMETER USAGE *** */
45
46 /* N....... TOTAL NUMBER OF RESIDUALS. */
47 /* P....... NUMBER OF PARAMETERS (COMPONENTS OF X) BEING ESTIMATED. */
48 /* X....... PARAMETER VECTOR BEING ESTIMATED (INPUT = INITIAL GUESS, */
49 /* OUTPUT = BEST VALUE FOUND). */
50 /* CALCR... SUBROUTINE FOR COMPUTING RESIDUAL VECTOR. */
51 /* CALCJ... SUBROUTINE FOR COMPUTING JACOBIAN MATRIX = MATRIX OF FIRST */
52 /* PARTIALS OF THE RESIDUAL VECTOR. */
53 /* IV...... INTEGER VALUES ARRAY. */
54 /* LIV..... LENGTH OF IV (SEE DISCUSSION BELOW). */
55 /* LV...... LENGTH OF V (SEE DISCUSSION BELOW). */
56 /* V....... FLOATING-POINT VALUES ARRAY. */
57 /* UI...... PASSED UNCHANGED TO CALCR AND CALCJ. */
58 /* UR...... PASSED UNCHANGED TO CALCR AND CALCJ. */
59 /* UF...... PASSED UNCHANGED TO CALCR AND CALCJ. */
60
61
62 /* *** DISCUSSION *** */
63
64 /* NOTE... NL2SOL (MENTIONED BELOW) IS A CODE FOR SOLVING */
65 /* NONLINEAR LEAST-SQUARES PROBLEMS. IT IS DESCRIBED IN */
66 /* ACM TRANS. MATH. SOFTWARE, VOL. 9, PP. 369-383 (AN ADAPTIVE */
67 /* NONLINEAR LEAST-SQUARES ALGORITHM, BY J.E. DENNIS, D.M. GAY, */
68 /* AND R.E. WELSCH). */
69
70 /* LIV GIVES THE LENGTH OF IV. IT MUST BE AT LEAST 82+P. IF NOT,
71 */
72 /* THEN DN2G RETURNS WITH IV(1) = 15. WHEN DN2G RETURNS, THE */
73 /* MINIMUM ACCEPTABLE VALUE OF LIV IS STORED IN IV(LASTIV) = IV(44),
74 */
75 /* (PROVIDED THAT LIV .GE. 44). */
76
77 /* LV GIVES THE LENGTH OF V. THE MINIMUM VALUE FOR LV IS */
78 /* LV0 = 105 + P*(N + 2*P + 17) + 2*N. IF LV IS SMALLER THAN THIS, */
79 /* THEN DN2G RETURNS WITH IV(1) = 16. WHEN DN2G RETURNS, THE */
80 /* MINIMUM ACCEPTABLE VALUE OF LV IS STORED IN IV(LASTV) = IV(45) */
81 /* (PROVIDED LIV .GE. 45). */
82
83 /* RETURN CODES AND CONVERGENCE TOLERANCES ARE THE SAME AS FOR */
84 /* NL2SOL, WITH SOME SMALL EXTENSIONS... IV(1) = 15 MEANS LIV WAS */
85 /* TOO SMALL. IV(1) = 16 MEANS LV WAS TOO SMALL. */
86
87 /* THERE ARE TWO NEW V INPUT COMPONENTS... V(LMAXS) = V(36) AND */
88 /* V(SCTOL) = V(37) SERVE WHERE V(LMAX0) AND V(RFCTOL) FORMERLY DID */
89 /* IN THE SINGULAR CONVERGENCE TEST -- SEE THE NL2SOL DOCUMENTATION.
90 */
91
92 /* *** DEFAULT VALUES *** */
93
94 /* DEFAULT VALUES ARE PROVIDED BY SUBROUTINE DIVSET, RATHER THAN */
95 /* DFAULT. THE CALLING SEQUENCE IS... */
96 /* CALL DIVSET(1, IV, LIV, LV, V) */
97 /* THE FIRST PARAMETER IS AN INTEGER 1. IF LIV AND LV ARE LARGE */
98 /* ENOUGH FOR DIVSET, THEN DIVSET SETS IV(1) TO 12. OTHERWISE IT */
99 /* SETS IV(1) TO 15 OR 16. CALLING DN2G WITH IV(1) = 0 CAUSES ALL
100 */
101 /* DEFAULT VALUES TO BE USED FOR THE INPUT COMPONENTS OF IV AND V. */
102 /* IF YOU FIRST CALL DIVSET, THEN SET IV(1) TO 13 AND CALL DN2G,
103 */
104 /* THEN STORAGE ALLOCATION ONLY WILL BE PERFORMED. IN PARTICULAR, */
105 /* IV(D) = IV(27), IV(J) = IV(70), AND IV(R) = IV(61) WILL BE SET */
106 /* TO THE FIRST SUBSCRIPT IN V OF THE SCALE VECTOR, THE JACOBIAN */
107 /* MATRIX, AND THE RESIDUAL VECTOR RESPECTIVELY, PROVIDED LIV AND LV
108 */
109 /* ARE LARGE ENOUGH. IF SO, THEN DN2G RETURNS WITH IV(1) = 14. */
110 /* WHEN CALLED WITH IV(1) = 14, DN2G ASSUMES THAT STORAGE HAS */
111 /* BEEN ALLOCATED, AND IT BEGINS THE MINIMIZATION ALGORITHM. */
112
113 /* *** SCALE VECTOR *** */
114
115 /* ONE DIFFERENCE WITH NL2SOL IS THAT THE SCALE VECTOR D IS */
116 /* STORED IN V, STARTING AT A DIFFERENT SUBSCRIPT. THE STARTING */
117 /* SUBSCRIPT VALUE IS STILL STORED IN IV(D) = IV(27). THE */
118 /* DISCUSSION OF DEFAULT VALUES ABOVE TELLS HOW TO HAVE IV(D) SET */
119 /* BEFORE THE ALGORITHM IS STARTED. */
120
121 /* *** REGRESSION DIAGNOSTICS *** */
122
123 /* IF IV(RDREQ) SO DICTATES, THEN ESTIMATES ARE COMPUTED OF THE */
124 /* INFLUENCE EACH RESIDUAL COMPONENT HAS ON THE FINAL PARAMETER */
125 /* ESTIMATE X. THE GENERAL IDEA IS THAT ONE MAY WISH TO EXAMINE */
126 /* RESIDUAL COMPONENTS (AND THE DATA BEHIND THEM) FOR WHICH THE */
127 /* INFLUENCE ESTIMATE IS SIGNIFICANTLY LARGER THAN MOST OF THE OTHER
128 */
129 /* INFLUENCE ESTIMATES. THESE ESTIMATES, HEREAFTER CALLED */
130 /* REGRESSION DIAGNOSTICS, ARE ONLY COMPUTED IF IV(RDREQ) = 2 OR 3. */
131 /* IN THIS CASE, FOR I = 1(1)N, */
132 /* SQRT( G(I)**T * H(I)**-1 * G(I) ) */
133 /* IS COMPUTED AND STORED IN V, STARTING AT V(IV(REGD)), WHERE */
134 /* RDREQ = 57 AND REGD = 67. HERE G(I) STANDS FOR THE GRADIENT */
135 /* RESULTING WHEN THE I-TH OBSERVATION IS DELETED AND H(I) STANDS */
136 /* FOR AN APPROXIMATION TO THE CORRESPONDING HESSIAN AT X, THE SOLU-
137 */
138 /* TION CORRESPONDING TO ALL OBSERVATIONS. (THIS APPROXIMATION IS */
139 /* OBTAINED BY SUBTRACTING THE FIRST-ORDER CONTRIBUTION OF THE I-TH */
140 /* OBSERVATION TO THE HESSIAN FROM A FINITE-DIFFERENCE HESSIAN */
141 /* APPROXIMATION. IF H IS INDEFINITE, THEN IV(REGD) IS SET TO -1. */
142 /* IF H(I) IS INDEFINITE, THEN -1 IS RETURNED AS THE DIAGNOSTIC FOR */
143 /* OBSERVATION I. IF NO DIAGNOSTICS ARE COMPUTED, PERHAPS BECAUSE */
144 /* OF A FAILURE TO CONVERGE, THEN IV(REGD) = 0 IS RETURNED.) */
145 /* PRINTING OF THE REGRESSION DIAGNOSTICS IS CONTROLLED BY */
146 /* IV(COVPRT) = IV(14)... IF IV(COVPRT) = 3, THEN BOTH THE */
147 /* COVARIANCE MATRIX AND THE REGRESSION DIAGNOSTICS ARE PRINTED. */
148 /* IV(COVPRT) = 2 CAUSES ONLY THE REGRESSION DIAGNOSTICS TO BE */
149 /* PRINTED, IV(COVPRT) = 1 CAUSES ONLY THE COVARIANCE MATRIX TO BE */
150 /* PRINTED, AND IV(COVPRT) = 0 CAUSES NEITHER TO BE PRINTED. */
151
152 /* RDREQ = 57 AND REGD = 67. */
153
154 /* *** GENERAL *** */
155
156 /* CODED BY DAVID M. GAY. */
157
158 /* +++++++++++++++++++++++++++ DECLARATIONS +++++++++++++++++++++++++++
159 */
160
161 /* *** EXTERNAL SUBROUTINES *** */
162
163 /* DIVSET.... PROVIDES DEFAULT IV AND V INPUT COMPONENTS. */
164 /* DRN2G... CARRIES OUT OPTIMIZATION ITERATIONS. */
165 /* DN2RDP... PRINTS REGRESSION DIAGNOSTICS. */
166
167 /* *** NO INTRINSIC FUNCTIONS *** */
168
169 /* *** LOCAL VARIABLES *** */
170
171
172 /* *** IV COMPONENTS *** */
173
174 /* /6 */
175 /* DATA D/27/, J/70/, NEXTV/47/, NFCALL/6/, NFGCAL/7/, R/61/, */
176 /* 1 REGD/67/, REGD0/82/, TOOBIG/2/, VNEED/4/ */
177 /* /7 */
178 /* / */
179 /* --------------------------------- BODY ------------------------------
180 */
181
182 /* Parameter adjustments */
183 --x;
184 --iv;
185 --v;
186 --ui;
187 --ur;
188
189 /* Function Body */
190 if (iv[1] == 0) {
191 divset_(&c__1, &iv[1], liv, lv, &v[1]);
192 }
193 iv1 = iv[1];
194 if (iv1 == 14) {
195 goto L10;
196 }
197 if (iv1 > 2 && iv1 < 12) {
198 goto L10;
199 }
200 if (iv1 == 12) {
201 iv[1] = 13;
202 }
203 if (iv[1] == 13) {
204 iv[4] = iv[4] + *p + *n * (*p + 2);
205 }
206 drn2g_(&x[1], &v[1], &iv[1], liv, lv, n, n, &n1, &n2, p, &v[1], &v[1], &v[
207 1], &x[1]);
208 if (iv[1] != 14) {
209 goto L999;
210 }
211
212 /* *** STORAGE ALLOCATION *** */
213
214 iv[27] = iv[47];
215 iv[61] = iv[27] + *p;
216 iv[82] = iv[61] + *n;
217 iv[70] = iv[82] + *n;
218 iv[47] = iv[70] + *n * *p;
219 if (iv1 == 13) {
220 goto L999;
221 }
222
223 L10:
224 d1 = iv[27];
225 dr1 = iv[70];
226 r1 = iv[61];
227 rd1 = iv[82];
228
229 L20:
230 drn2g_(&v[d1], &v[dr1], &iv[1], liv, lv, n, n, &n1, &n2, p, &v[r1], &v[
231 rd1], &v[1], &x[1]);
232 if ((i__1 = iv[1] - 2) < 0) {
233 goto L30;
234 } else if (i__1 == 0) {
235 goto L50;
236 } else {
237 goto L60;
238 }
239
240 /* *** NEW FUNCTION VALUE (R VALUE) NEEDED *** */
241
242 L30:
243 nf = iv[6];
244 (*calcr)(n, p, &x[1], &nf, &v[r1], &ui[1], &ur[1], uf);
245 if (nf > 0) {
246 goto L40;
247 }
248 iv[2] = 1;
249 goto L20;
250 L40:
251 if (iv[1] > 0) {
252 goto L20;
253 }
254
255 /* *** COMPUTE DR = GRADIENT OF R COMPONENTS *** */
256
257 L50:
258 (*calcj)(n, p, &x[1], &iv[7], &v[dr1], &ui[1], &ur[1], uf);
259 if (iv[7] == 0) {
260 iv[2] = 1;
261 }
262 goto L20;
263
264 /* *** INDICATE WHETHER THE REGRESSION DIAGNOSTIC ARRAY WAS COMPUTED */
265 /* *** AND PRINT IT IF SO REQUESTED... */
266
267 L60:
268 if (iv[67] > 0) {
269 iv[67] = rd1;
270 }
271 dn2rdp_(&iv[1], liv, lv, n, &v[rd1], &v[1]);
272
273 L999:
274 return 0;
275
276 /* *** LAST LINE OF DN2G FOLLOWS *** */
277 } /* dn2g_ */
278
279