1 /* dg7lit.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 static doublereal c_b11 = 0.;
12 static doublereal c_b23 = -1.;
13 static doublereal c_b43 = 1.;
14 static integer c__9 = 9;
15 static integer c__4 = 4;
16
dg7lit_(d__,g,iv,liv,lv,p,ps,v,x,y)17 /* Subroutine */ int dg7lit_(d__, g, iv, liv, lv, p, ps, v, x, y)
18 doublereal *d__, *g;
19 integer *iv, *liv, *lv, *p, *ps;
20 doublereal *v, *x, *y;
21 {
22 /* System generated locals */
23 integer i__1, i__2;
24 doublereal d__1;
25
26 /* Local variables */
27 static integer pp1o2, lmat1, rmat1, temp1, temp2, ipiv1, step1;
28 static doublereal e;
29 static integer i__, j, k, l;
30 static doublereal t;
31 static integer h1, dummy, s1;
32 static doublereal t1;
33 static integer w1;
34 extern logical stopx_();
35 extern doublereal dr7mdc_();
36 extern /* Subroutine */ int df7hes_();
37 extern doublereal dd7tpr_();
38 extern /* Subroutine */ int da7sst_(), dl7vml_(), dv7scp_();
39 static integer g01;
40 extern doublereal dv2nrm_();
41 extern /* Subroutine */ int dg7qts_(), dl7mst_(), dl7sqr_();
42 extern doublereal dl7svn_();
43 extern /* Subroutine */ int dl7tvm_(), dl7srt_(), ds7lup_(), ds7lvm_(),
44 dv2axy_(), dv7cpy_();
45 extern doublereal dl7svx_();
46 static integer x01;
47 extern /* Subroutine */ int dparck_();
48 static integer hc1;
49 extern doublereal drldst_();
50 extern /* Subroutine */ int ditsum_();
51 static integer stpmod, lstgst, rstrst;
52 static doublereal sttsst;
53 static integer dig1, qtr1;
54
55
56 /* *** CARRY OUT NL2SOL-LIKE ITERATIONS FOR GENERALIZED LINEAR *** */
57 /* *** REGRESSION PROBLEMS (AND OTHERS OF SIMILAR STRUCTURE) *** */
58
59 /* *** PARAMETER DECLARATIONS *** */
60
61
62 /* -------------------------- PARAMETER USAGE --------------------------
63 */
64
65 /* D.... SCALE VECTOR. */
66 /* IV... INTEGER VALUE ARRAY. */
67 /* LIV.. LENGTH OF IV. MUST BE AT LEAST 82. */
68 /* LH... LENGTH OF H = P*(P+1)/2. */
69 /* LV... LENGTH OF V. MUST BE AT LEAST P*(3*P + 19)/2 + 7. */
70 /* G.... GRADIENT AT X (WHEN IV(1) = 2). */
71 /* P.... NUMBER OF PARAMETERS (COMPONENTS IN X). */
72 /* PS... NUMBER OF NONZERO ROWS AND COLUMNS IN S. */
73 /* V.... FLOATING-POINT VALUE ARRAY. */
74 /* X.... PARAMETER VECTOR. */
75 /* Y.... PART OF YIELD VECTOR (WHEN IV(1)= 2, SCRATCH OTHERWISE). */
76
77 /* *** DISCUSSION *** */
78
79 /* DG7LIT PERFORMS NL2SOL-LIKE ITERATIONS FOR A VARIETY OF */
80 /* REGRESSION PROBLEMS THAT ARE SIMILAR TO NONLINEAR LEAST-SQUARES */
81 /* IN THAT THE HESSIAN IS THE SUM OF TWO TERMS, A READILY-COMPUTED */
82 /* FIRST-ORDER TERM AND A SECOND-ORDER TERM. THE CALLER SUPPLIES */
83 /* THE FIRST-ORDER TERM OF THE HESSIAN IN HC (LOWER TRIANGLE, STORED
84 */
85 /* COMPACTLY BY ROWS IN V, STARTING AT IV(HC)), AND DG7LIT BUILDS AN
86 */
87 /* APPROXIMATION, S, TO THE SECOND-ORDER TERM. THE CALLER ALSO */
88 /* PROVIDES THE FUNCTION VALUE, GRADIENT, AND PART OF THE YIELD */
89 /* VECTOR USED IN UPDATING S. DG7LIT DECIDES DYNAMICALLY WHETHER OR */
90 /* NOT TO USE S WHEN CHOOSING THE NEXT STEP TO TRY... THE HESSIAN */
91 /* APPROXIMATION USED IS EITHER HC ALONE (GAUSS-NEWTON MODEL) OR */
92 /* HC + S (AUGMENTED MODEL). */
93
94 /* IF PS .LT. P, THEN ROWS AND COLUMNS PS+1...P OF S ARE KEPT */
95 /* CONSTANT. THEY WILL BE ZERO UNLESS THE CALLER SETS IV(INITS) TO */
96 /* 1 OR 2 AND SUPPLIES NONZERO VALUES FOR THEM, OR THE CALLER SETS */
97 /* IV(INITS) TO 3 OR 4 AND THE FINITE-DIFFERENCE INITIAL S THEN */
98 /* COMPUTED HAS NONZERO VALUES IN THESE ROWS. */
99
100 /* IF IV(INITS) IS 3 OR 4, THEN THE INITIAL S IS COMPUTED BY */
101 /* FINITE DIFFERENCES. 3 MEANS USE FUNCTION DIFFERENCES, 4 MEANS */
102 /* USE GRADIENT DIFFERENCES. FINITE DIFFERENCING IS DONE THE SAME */
103 /* WAY AS IN COMPUTING A COVARIANCE MATRIX (WITH IV(COVREQ) = -1, -2,
104 */
105 /* 1, OR 2). */
106
107 /* FOR UPDATING S,DG7LIT ASSUMES THAT THE GRADIENT HAS THE FORM */
108 /* OF A SUM OVER I OF RHO(I,X)*GRAD(R(I,X)), WHERE GRAD DENOTES THE */
109 /* GRADIENT WITH RESPECT TO X. THE TRUE SECOND-ORDER TERM THEN IS */
110 /* THE SUM OVER I OF RHO(I,X)*HESSIAN(R(I,X)). IF X = X0 + STEP, */
111 /* THEN WE WISH TO UPDATE S SO THAT S*STEP IS THE SUM OVER I OF */
112 /* RHO(I,X)*(GRAD(R(I,X)) - GRAD(R(I,X0))). THE CALLER MUST SUPPLY */
113 /* PART OF THIS IN Y, NAMELY THE SUM OVER I OF */
114 /* RHO(I,X)*GRAD(R(I,X0)), WHEN CALLING DG7LIT WITH IV(1) = 2 AND */
115 /* IV(MODE) = 0 (WHERE MODE = 38). G THEN CONTANS THE OTHER PART, */
116 /* SO THAT THE DESIRED YIELD VECTOR IS G - Y. IF PS .LT. P, THEN */
117 /* THE ABOVE DISCUSSION APPLIES ONLY TO THE FIRST PS COMPONENTS OF */
118 /* GRAD(R(I,X)), STEP, AND Y. */
119
120 /* PARAMETERS IV, P, V, AND X ARE THE SAME AS THE CORRESPONDING */
121 /* ONES TO NL2SOL (WHICH SEE), EXCEPT THAT V CAN BE SHORTER */
122 /* (SINCE THE PART OF V THAT NL2SOL USES FOR STORING D, J, AND R IS */
123 /* NOT NEEDED). MOREOVER, COMPARED WITH NL2SOL, IV(1) MAY HAVE THE */
124 /* TWO ADDITIONAL OUTPUT VALUES 1 AND 2, WHICH ARE EXPLAINED BELOW, */
125 /* AS IS THE USE OF IV(TOOBIG) AND IV(NFGCAL). THE VALUES IV(D), */
126 /* IV(J), AND IV(R), WHICH ARE OUTPUT VALUES FROM NL2SOL (AND */
127 /* NL2SNO), ARE NOT REFERENCED BY DG7LIT OR THE SUBROUTINES IT CALLS.
128 */
129
130 /* WHEN DG7LIT IS FIRST CALLED, I.E., WHEN DG7LIT IS CALLED WITH */
131 /* IV(1) = 0 OR 12, V(F), G, AND HC NEED NOT BE INITIALIZED. TO */
132 /* OBTAIN THESE STARTING VALUES,DG7LIT RETURNS FIRST WITH IV(1) = 1,
133 */
134 /* THEN WITH IV(1) = 2, WITH IV(MODE) = -1 IN BOTH CASES. ON */
135 /* SUBSEQUENT RETURNS WITH IV(1) = 2, IV(MODE) = 0 IMPLIES THAT */
136 /* Y MUST ALSO BE SUPPLIED. (NOTE THAT Y IS USED FOR SCRATCH -- ITS
137 */
138 /* INPUT CONTENTS ARE LOST. BY CONTRAST, HC IS NEVER CHANGED.) */
139 /* ONCE CONVERGENCE HAS BEEN OBTAINED, IV(RDREQ) AND IV(COVREQ) MAY */
140 /* IMPLY THAT A FINITE-DIFFERENCE HESSIAN SHOULD BE COMPUTED FOR USE
141 */
142 /* IN COMPUTING A COVARIANCE MATRIX. IN THIS CASE DG7LIT WILL MAKE A
143 */
144 /* NUMBER OF RETURNS WITH IV(1) = 1 OR 2 AND IV(MODE) POSITIVE. */
145 /* WHEN IV(MODE) IS POSITIVE, Y SHOULD NOT BE CHANGED. */
146
147 /* IV(1) = 1 MEANS THE CALLER SHOULD SET V(F) (I.E., V(10)) TO F(X), THE
148 */
149 /* FUNCTION VALUE AT X, AND CALL DG7LIT AGAIN, HAVING CHANGED
150 */
151 /* NONE OF THE OTHER PARAMETERS. AN EXCEPTION OCCURS IF F(X)
152 */
153 /* CANNOT BE EVALUATED (E.G. IF OVERFLOW WOULD OCCUR), WHICH
154 */
155 /* MAY HAPPEN BECAUSE OF AN OVERSIZED STEP. IN THIS CASE */
156 /* THE CALLER SHOULD SET IV(TOOBIG) = IV(2) TO 1, WHICH WILL
157 */
158 /* CAUSE DG7LIT TO IGNORE V(F) AND TRY A SMALLER STEP. NOTE
159 */
160 /* THAT THE CURRENT FUNCTION EVALUATION COUNT IS AVAILABLE */
161 /* IN IV(NFCALL) = IV(6). THIS MAY BE USED TO IDENTIFY */
162 /* WHICH COPY OF SAVED INFORMATION SHOULD BE USED IN COM- */
163 /* PUTING G, HC, AND Y THE NEXT TIME DG7LIT RETURNS WITH */
164 /* IV(1) = 2. SEE MLPIT FOR AN EXAMPLE OF THIS. */
165 /* IV(1) = 2 MEANS THE CALLER SHOULD SET G TO G(X), THE GRADIENT OF F AT
166 */
167 /* X. THE CALLER SHOULD ALSO SET HC TO THE GAUSS-NEWTON */
168 /* HESSIAN AT X. IF IV(MODE) = 0, THEN THE CALLER SHOULD */
169 /* ALSO COMPUTE THE PART OF THE YIELD VECTOR DESCRIBED ABOVE.
170 */
171 /* THE CALLER SHOULD THEN CALL DG7LIT AGAIN (WITH IV(1) = 2).
172 */
173 /* THE CALLER MAY ALSO CHANGE D AT THIS TIME, BUT SHOULD NOT
174 */
175 /* CHANGE X. NOTE THAT IV(NFGCAL) = IV(7) CONTAINS THE */
176 /* VALUE THAT IV(NFCALL) HAD DURING THE RETURN WITH */
177 /* IV(1) = 1 IN WHICH X HAD THE SAME VALUE AS IT NOW HAS. */
178 /* IV(NFGCAL) IS EITHER IV(NFCALL) OR IV(NFCALL) - 1. MLPIT
179 */
180 /* IS AN EXAMPLE WHERE THIS INFORMATION IS USED. IF G OR HC
181 */
182 /* CANNOT BE EVALUATED AT X, THEN THE CALLER MAY SET */
183 /* IV(TOOBIG) TO 1, IN WHICH CASE DG7LIT WILL RETURN WITH */
184 /* IV(1) = 15. */
185
186 /* *** GENERAL *** */
187
188 /* CODED BY DAVID M. GAY. */
189 /* THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH */
190 /* SUPPORTED IN PART BY D.O.E. GRANT EX-76-A-01-2295 TO MIT/CCREMS. */
191
192 /* (SEE NL2SOL FOR REFERENCES.) */
193
194 /* +++++++++++++++++++++++++++ DECLARATIONS ++++++++++++++++++++++++++++
195 */
196
197 /* *** LOCAL VARIABLES *** */
198
199
200 /* *** CONSTANTS *** */
201
202
203 /* *** EXTERNAL FUNCTIONS AND SUBROUTINES *** */
204
205
206 /* DA7SST.... ASSESSES CANDIDATE STEP. */
207 /* DD7TPR... RETURNS INNER PRODUCT OF TWO VECTORS. */
208 /* DF7HES.... COMPUTE FINITE-DIFFERENCE HESSIAN (FOR COVARIANCE). */
209 /* DG7QTS.... COMPUTES GOLDFELD-QUANDT-TROTTER STEP (AUGMENTED MODEL). */
210 /* DITSUM.... PRINTS ITERATION SUMMARY AND INFO ON INITIAL AND FINAL X. */
211 /* DL7MST... COMPUTES LEVENBERG-MARQUARDT STEP (GAUSS-NEWTON MODEL). */
212 /* DL7SRT.... COMPUTES CHOLESKY FACTOR OF (LOWER TRIANG. OF) SYM. MATRIX.
213 */
214 /* DL7SQR... COMPUTES L * L**T FROM LOWER TRIANGULAR MATRIX L. */
215 /* DL7TVM... COMPUTES L**T * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX.
216 */
217 /* DL7SVX... ESTIMATES LARGEST SING. VALUE OF LOWER TRIANG. MATRIX. */
218 /* DL7SVN... ESTIMATES SMALLEST SING. VALUE OF LOWER TRIANG. MATRIX. */
219 /* DL7VML.... COMPUTES L * V, V = VECTOR, L = LOWER TRIANGULAR MATRIX. */
220 /* DPARCK.... CHECK VALIDITY OF IV AND V INPUT COMPONENTS. */
221 /* DRLDST... COMPUTES V(RELDX) = RELATIVE STEP SIZE. */
222 /* DR7MDC... RETURNS MACHINE-DEPENDENT CONSTANTS. */
223 /* DS7LUP... PERFORMS QUASI-NEWTON UPDATE ON COMPACTLY STORED LOWER TRI-
224 */
225 /* ANGLE OF A SYMMETRIC MATRIX. */
226 /* STOPX.... RETURNS .TRUE. IF THE BREAK KEY HAS BEEN PRESSED. */
227 /* DV2AXY.... COMPUTES SCALAR TIMES ONE VECTOR PLUS ANOTHER. */
228 /* DV7CPY.... COPIES ONE VECTOR TO ANOTHER. */
229 /* DV7SCP... SETS ALL ELEMENTS OF A VECTOR TO A SCALAR. */
230 /* DV2NRM... RETURNS THE 2-NORM OF A VECTOR. */
231
232 /* *** SUBSCRIPTS FOR IV AND V *** */
233
234
235 /* *** IV SUBSCRIPT VALUES *** */
236
237 /* /6 */
238 /* DATA CNVCOD/55/, COVMAT/26/, COVREQ/15/, DIG/37/, FDH/74/, H/56/,
239 */
240 /* 1 HC/71/, IERR/75/, INITS/25/, IPIVOT/76/, IRC/29/, KAGQT/33/,
241 */
242 /* 2 KALM/34/, LMAT/42/, MODE/35/, MODEL/5/, MXFCAL/17/, */
243 /* 3 MXITER/18/, NEXTV/47/, NFCALL/6/, NFGCAL/7/, NFCOV/52/, */
244 /* 4 NGCOV/53/, NGCALL/30/, NITER/31/, QTR/77/, RADINC/8/, */
245 /* 5 RDREQ/57/, REGD/67/, RESTOR/9/, RMAT/78/, S/62/, STEP/40/, */
246 /* 6 STGLIM/11/, STLSTG/41/, SUSED/64/, SWITCH/12/, TOOBIG/2/, */
247 /* 7 VNEED/4/, VSAVE/60/, W/65/, XIRC/13/, X0/43/ */
248 /* /7 */
249 /* / */
250
251 /* *** V SUBSCRIPT VALUES *** */
252
253 /* /6 */
254 /* DATA COSMIN/47/, DGNORM/1/, DSTNRM/2/, F/10/, FDIF/11/, FUZZ/45/,
255 */
256 /* 1 F0/13/, GTSTEP/4/, INCFAC/23/, LMAX0/35/, LMAXS/36/, */
257 /* 2 NVSAVE/9/, PHMXFC/21/, PREDUC/7/, RADFAC/16/, RADIUS/8/, */
258 /* 3 RAD0/9/, RCOND/53/, RELDX/17/, SIZE/55/, STPPAR/5/, */
259 /* 4 TUNER4/29/, TUNER5/30/, WSCALE/56/ */
260 /* /7 */
261 /* / */
262
263
264 /* /6 */
265 /* DATA HALF/0.5D+0/, NEGONE/-1.D+0/, ONE/1.D+0/, ONEP2/1.2D+0/, */
266 /* 1 ZERO/0.D+0/ */
267 /* /7 */
268 /* / */
269
270 /* +++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
271 */
272
273 /* Parameter adjustments */
274 --iv;
275 --v;
276 --y;
277 --x;
278 --g;
279 --d__;
280
281 /* Function Body */
282 i__ = iv[1];
283 if (i__ == 1) {
284 goto L40;
285 }
286 if (i__ == 2) {
287 goto L50;
288 }
289
290 if (i__ == 12 || i__ == 13) {
291 iv[4] = iv[4] + *p * (*p * 3 + 19) / 2 + 7;
292 }
293 dparck_(&c__1, &d__[1], &iv[1], liv, lv, p, &v[1]);
294 i__ = iv[1] - 2;
295 if (i__ > 12) {
296 goto L999;
297 }
298 switch ((int)i__) {
299 case 1: goto L290;
300 case 2: goto L290;
301 case 3: goto L290;
302 case 4: goto L290;
303 case 5: goto L290;
304 case 6: goto L290;
305 case 7: goto L170;
306 case 8: goto L120;
307 case 9: goto L170;
308 case 10: goto L10;
309 case 11: goto L10;
310 case 12: goto L20;
311 }
312
313 /* *** STORAGE ALLOCATION *** */
314
315 L10:
316 pp1o2 = *p * (*p + 1) / 2;
317 iv[62] = iv[42] + pp1o2;
318 iv[43] = iv[62] + pp1o2;
319 iv[40] = iv[43] + *p;
320 iv[41] = iv[40] + *p;
321 iv[37] = iv[41] + *p;
322 iv[65] = iv[37] + *p;
323 iv[56] = iv[65] + (*p << 2) + 7;
324 iv[47] = iv[56] + pp1o2;
325 if (iv[1] != 13) {
326 goto L20;
327 }
328 iv[1] = 14;
329 goto L999;
330
331 /* *** INITIALIZATION *** */
332
333 L20:
334 iv[31] = 0;
335 iv[6] = 1;
336 iv[30] = 1;
337 iv[7] = 1;
338 iv[35] = -1;
339 iv[11] = 2;
340 iv[2] = 0;
341 iv[55] = 0;
342 iv[26] = 0;
343 iv[52] = 0;
344 iv[53] = 0;
345 iv[8] = 0;
346 iv[9] = 0;
347 iv[74] = 0;
348 v[9] = 0.;
349 v[5] = 0.;
350 v[8] = v[35] / (v[21] + 1.);
351
352 /* *** SET INITIAL MODEL AND S MATRIX *** */
353
354 iv[5] = 1;
355 if (iv[62] < 0) {
356 goto L999;
357 }
358 if (iv[25] > 1) {
359 iv[5] = 2;
360 }
361 s1 = iv[62];
362 if (iv[25] == 0 || iv[25] > 2) {
363 i__1 = *p * (*p + 1) / 2;
364 dv7scp_(&i__1, &v[s1], &c_b11);
365 }
366 iv[1] = 1;
367 j = iv[76];
368 if (j <= 0) {
369 goto L999;
370 }
371 i__1 = *p;
372 for (i__ = 1; i__ <= i__1; ++i__) {
373 iv[j] = i__;
374 ++j;
375 /* L30: */
376 }
377 goto L999;
378
379 /* *** NEW FUNCTION VALUE *** */
380
381 L40:
382 if (iv[35] == 0) {
383 goto L290;
384 }
385 if (iv[35] > 0) {
386 goto L520;
387 }
388
389 iv[1] = 2;
390 if (iv[2] == 0) {
391 goto L999;
392 }
393 iv[1] = 63;
394 goto L999;
395
396 /* *** NEW GRADIENT *** */
397
398 L50:
399 iv[34] = -1;
400 iv[33] = -1;
401 iv[74] = 0;
402 if (iv[35] > 0) {
403 goto L520;
404 }
405
406 /* *** MAKE SURE GRADIENT COULD BE COMPUTED *** */
407
408 if (iv[2] == 0) {
409 goto L60;
410 }
411 iv[1] = 65;
412 goto L999;
413 L60:
414 if (iv[71] <= 0 && iv[78] <= 0) {
415 goto L610;
416 }
417
418 /* *** COMPUTE D**-1 * GRADIENT *** */
419
420 dig1 = iv[37];
421 k = dig1;
422 i__1 = *p;
423 for (i__ = 1; i__ <= i__1; ++i__) {
424 v[k] = g[i__] / d__[i__];
425 ++k;
426 /* L70: */
427 }
428 v[1] = dv2nrm_(p, &v[dig1]);
429
430 if (iv[55] != 0) {
431 goto L510;
432 }
433 if (iv[35] == 0) {
434 goto L440;
435 }
436 iv[35] = 0;
437 v[13] = v[10];
438 if (iv[25] <= 2) {
439 goto L100;
440 }
441
442 /* *** ARRANGE FOR FINITE-DIFFERENCE INITIAL S *** */
443
444 iv[13] = iv[15];
445 iv[15] = -1;
446 if (iv[25] > 3) {
447 iv[15] = 1;
448 }
449 iv[55] = 70;
450 goto L530;
451
452 /* *** COME TO NEXT STMT AFTER COMPUTING F.D. HESSIAN FOR INIT. S ***
453 */
454
455 L80:
456 iv[55] = 0;
457 iv[35] = 0;
458 iv[52] = 0;
459 iv[53] = 0;
460 iv[15] = iv[13];
461 s1 = iv[62];
462 pp1o2 = *ps * (*ps + 1) / 2;
463 hc1 = iv[71];
464 if (hc1 <= 0) {
465 goto L90;
466 }
467 dv2axy_(&pp1o2, &v[s1], &c_b23, &v[hc1], &v[h1]);
468 goto L100;
469 L90:
470 rmat1 = iv[78];
471 dl7sqr_(ps, &v[s1], &v[rmat1]);
472 dv2axy_(&pp1o2, &v[s1], &c_b23, &v[s1], &v[h1]);
473 L100:
474 iv[1] = 2;
475
476
477 /* ----------------------------- MAIN LOOP -----------------------------
478 */
479
480
481 /* *** PRINT ITERATION SUMMARY, CHECK ITERATION LIMIT *** */
482
483 L110:
484 ditsum_(&d__[1], &g[1], &iv[1], liv, lv, p, &v[1], &x[1]);
485 L120:
486 k = iv[31];
487 if (k < iv[18]) {
488 goto L130;
489 }
490 iv[1] = 10;
491 goto L999;
492 L130:
493 iv[31] = k + 1;
494
495 /* *** UPDATE RADIUS *** */
496
497 if (k == 0) {
498 goto L150;
499 }
500 step1 = iv[40];
501 i__1 = *p;
502 for (i__ = 1; i__ <= i__1; ++i__) {
503 v[step1] = d__[i__] * v[step1];
504 ++step1;
505 /* L140: */
506 }
507 step1 = iv[40];
508 t = v[16] * dv2nrm_(p, &v[step1]);
509 if (v[16] < 1. || t > v[8]) {
510 v[8] = t;
511 }
512
513 /* *** INITIALIZE FOR START OF NEXT ITERATION *** */
514
515 L150:
516 x01 = iv[43];
517 v[13] = v[10];
518 iv[29] = 4;
519 iv[56] = -abs(iv[56]);
520 iv[64] = iv[5];
521
522 /* *** COPY X TO X0 *** */
523
524 dv7cpy_(p, &v[x01], &x[1]);
525
526 /* *** CHECK STOPX AND FUNCTION EVALUATION LIMIT *** */
527
528 L160:
529 if (! stopx_(&dummy)) {
530 goto L180;
531 }
532 iv[1] = 11;
533 goto L190;
534
535 /* *** COME HERE WHEN RESTARTING AFTER FUNC. EVAL. LIMIT OR STOPX. */
536
537 L170:
538 if (v[10] >= v[13]) {
539 goto L180;
540 }
541 v[16] = 1.;
542 k = iv[31];
543 goto L130;
544
545 L180:
546 if (iv[6] < iv[17] + iv[52]) {
547 goto L200;
548 }
549 iv[1] = 9;
550 L190:
551 if (v[10] >= v[13]) {
552 goto L999;
553 }
554
555 /* *** IN CASE OF STOPX OR FUNCTION EVALUATION LIMIT WITH */
556 /* *** IMPROVED V(F), EVALUATE THE GRADIENT AT X. */
557
558 iv[55] = iv[1];
559 goto L430;
560
561 /* . . . . . . . . . . . . . COMPUTE CANDIDATE STEP . . . . . . . . . .
562 */
563
564 L200:
565 step1 = iv[40];
566 w1 = iv[65];
567 h1 = iv[56];
568 t1 = 1.;
569 if (iv[5] == 2) {
570 goto L210;
571 }
572 t1 = 0.;
573
574 /* *** COMPUTE LEVENBERG-MARQUARDT STEP IF POSSIBLE... */
575
576 rmat1 = iv[78];
577 if (rmat1 <= 0) {
578 goto L210;
579 }
580 qtr1 = iv[77];
581 if (qtr1 <= 0) {
582 goto L210;
583 }
584 ipiv1 = iv[76];
585 dl7mst_(&d__[1], &g[1], &iv[75], &iv[ipiv1], &iv[34], p, &v[qtr1], &v[
586 rmat1], &v[step1], &v[1], &v[w1]);
587 /* *** H IS STORED IN THE END OF W AND HAS JUST BEEN OVERWRITTEN,
588 */
589 /* *** SO WE MARK IT INVALID... */
590 iv[56] = -abs(h1);
591 /* *** EVEN IF H WERE STORED ELSEWHERE, IT WOULD BE NECESSARY TO */
592 /* *** MARK INVALID THE INFORMATION DG7QTS MAY HAVE STORED IN V...
593 */
594 iv[33] = -1;
595 goto L260;
596
597 L210:
598 if (h1 > 0) {
599 goto L250;
600 }
601
602 /* *** SET H TO D**-1 * (HC + T1*S) * D**-1. *** */
603
604 h1 = -h1;
605 iv[56] = h1;
606 iv[74] = 0;
607 j = iv[71];
608 if (j > 0) {
609 goto L220;
610 }
611 j = h1;
612 rmat1 = iv[78];
613 dl7sqr_(p, &v[h1], &v[rmat1]);
614 L220:
615 s1 = iv[62];
616 i__1 = *p;
617 for (i__ = 1; i__ <= i__1; ++i__) {
618 t = 1. / d__[i__];
619 i__2 = i__;
620 for (k = 1; k <= i__2; ++k) {
621 v[h1] = t * (v[j] + t1 * v[s1]) / d__[k];
622 ++j;
623 ++h1;
624 ++s1;
625 /* L230: */
626 }
627 /* L240: */
628 }
629 h1 = iv[56];
630 iv[33] = -1;
631
632 /* *** COMPUTE ACTUAL GOLDFELD-QUANDT-TROTTER STEP *** */
633
634 L250:
635 dig1 = iv[37];
636 lmat1 = iv[42];
637 dg7qts_(&d__[1], &v[dig1], &v[h1], &iv[33], &v[lmat1], p, &v[step1], &v[1]
638 , &v[w1]);
639 if (iv[34] > 0) {
640 iv[34] = 0;
641 }
642
643 L260:
644 if (iv[29] != 6) {
645 goto L270;
646 }
647 if (iv[9] != 2) {
648 goto L290;
649 }
650 rstrst = 2;
651 goto L300;
652
653 /* *** CHECK WHETHER EVALUATING F(X0 + STEP) LOOKS WORTHWHILE *** */
654
655 L270:
656 iv[2] = 0;
657 if (v[2] <= 0.) {
658 goto L290;
659 }
660 if (iv[29] != 5) {
661 goto L280;
662 }
663 if (v[16] <= 1.) {
664 goto L280;
665 }
666 if (v[7] > v[11] * 1.2) {
667 goto L280;
668 }
669 x01 = iv[43];
670 step1 = iv[40];
671 dv2axy_(p, &v[step1], &c_b23, &v[x01], &x[1]);
672 if (iv[9] != 2) {
673 goto L290;
674 }
675 rstrst = 0;
676 goto L300;
677
678 /* *** COMPUTE F(X0 + STEP) *** */
679
680 L280:
681 x01 = iv[43];
682 step1 = iv[40];
683 dv2axy_(p, &x[1], &c_b43, &v[step1], &v[x01]);
684 ++iv[6];
685 iv[1] = 1;
686 goto L999;
687
688 /* . . . . . . . . . . . . . ASSESS CANDIDATE STEP . . . . . . . . . . .
689 */
690
691 L290:
692 rstrst = 3;
693 L300:
694 x01 = iv[43];
695 v[17] = drldst_(p, &d__[1], &x[1], &v[x01]);
696 da7sst_(&iv[1], liv, lv, &v[1]);
697 step1 = iv[40];
698 lstgst = iv[41];
699 i__ = iv[9] + 1;
700 switch ((int)i__) {
701 case 1: goto L340;
702 case 2: goto L310;
703 case 3: goto L320;
704 case 4: goto L330;
705 }
706 L310:
707 dv7cpy_(p, &x[1], &v[x01]);
708 goto L340;
709 L320:
710 dv7cpy_(p, &v[lstgst], &v[step1]);
711 goto L340;
712 L330:
713 dv7cpy_(p, &v[step1], &v[lstgst]);
714 dv2axy_(p, &x[1], &c_b43, &v[step1], &v[x01]);
715 v[17] = drldst_(p, &d__[1], &x[1], &v[x01]);
716 iv[9] = rstrst;
717
718 /* *** IF NECESSARY, SWITCH MODELS *** */
719
720 L340:
721 if (iv[12] == 0) {
722 goto L350;
723 }
724 iv[56] = -abs(iv[56]);
725 iv[64] += 2;
726 l = iv[60];
727 dv7cpy_(&c__9, &v[1], &v[l]);
728 L350:
729 l = iv[29] - 4;
730 stpmod = iv[5];
731 if (l > 0) {
732 switch ((int)l) {
733 case 1: goto L370;
734 case 2: goto L380;
735 case 3: goto L390;
736 case 4: goto L390;
737 case 5: goto L390;
738 case 6: goto L390;
739 case 7: goto L390;
740 case 8: goto L390;
741 case 9: goto L500;
742 case 10: goto L440;
743 }
744 }
745
746 /* *** DECIDE WHETHER TO CHANGE MODELS *** */
747
748 e = v[7] - v[11];
749 s1 = iv[62];
750 ds7lvm_(ps, &y[1], &v[s1], &v[step1]);
751 sttsst = dd7tpr_(ps, &v[step1], &y[1]) * .5;
752 if (iv[5] == 1) {
753 sttsst = -sttsst;
754 }
755 if ((d__1 = e + sttsst, abs(d__1)) * v[45] >= abs(e)) {
756 goto L360;
757 }
758
759 /* *** SWITCH MODELS *** */
760
761 iv[5] = 3 - iv[5];
762 if (-2 < l) {
763 goto L400;
764 }
765 iv[56] = -abs(iv[56]);
766 iv[64] += 2;
767 l = iv[60];
768 dv7cpy_(&c__9, &v[l], &v[1]);
769 goto L160;
770
771 L360:
772 if (-3 < l) {
773 goto L400;
774 }
775
776 /* *** RECOMPUTE STEP WITH NEW RADIUS *** */
777
778 L370:
779 v[8] = v[16] * v[2];
780 goto L160;
781
782 /* *** COMPUTE STEP OF LENGTH V(LMAXS) FOR SINGULAR CONVERGENCE TEST */
783
784 L380:
785 v[8] = v[36];
786 goto L200;
787
788 /* *** CONVERGENCE OR FALSE CONVERGENCE *** */
789
790 L390:
791 iv[55] = l;
792 if (v[10] >= v[13]) {
793 goto L510;
794 }
795 if (iv[13] == 14) {
796 goto L510;
797 }
798 iv[13] = 14;
799
800 /* . . . . . . . . . . . . PROCESS ACCEPTABLE STEP . . . . . . . . . . .
801 */
802
803 L400:
804 iv[26] = 0;
805 iv[67] = 0;
806
807 /* *** SEE WHETHER TO SET V(RADFAC) BY GRADIENT TESTS *** */
808
809 if (iv[29] != 3) {
810 goto L430;
811 }
812 step1 = iv[40];
813 temp1 = iv[41];
814 temp2 = iv[65];
815
816 /* *** SET TEMP1 = HESSIAN * STEP FOR USE IN GRADIENT TESTS *** */
817
818 hc1 = iv[71];
819 if (hc1 <= 0) {
820 goto L410;
821 }
822 ds7lvm_(p, &v[temp1], &v[hc1], &v[step1]);
823 goto L420;
824 L410:
825 rmat1 = iv[78];
826 dl7tvm_(p, &v[temp1], &v[rmat1], &v[step1]);
827 dl7vml_(p, &v[temp1], &v[rmat1], &v[temp1]);
828
829 L420:
830 if (stpmod == 1) {
831 goto L430;
832 }
833 s1 = iv[62];
834 ds7lvm_(ps, &v[temp2], &v[s1], &v[step1]);
835 dv2axy_(ps, &v[temp1], &c_b43, &v[temp2], &v[temp1]);
836
837 /* *** SAVE OLD GRADIENT AND COMPUTE NEW ONE *** */
838
839 L430:
840 ++iv[30];
841 g01 = iv[65];
842 dv7cpy_(p, &v[g01], &g[1]);
843 iv[1] = 2;
844 iv[2] = 0;
845 goto L999;
846
847 /* *** INITIALIZATIONS -- G0 = G - G0, ETC. *** */
848
849 L440:
850 g01 = iv[65];
851 dv2axy_(p, &v[g01], &c_b23, &v[g01], &g[1]);
852 step1 = iv[40];
853 temp1 = iv[41];
854 temp2 = iv[65];
855 if (iv[29] != 3) {
856 goto L470;
857 }
858
859 /* *** SET V(RADFAC) BY GRADIENT TESTS *** */
860
861 /* *** SET TEMP1 = D**-1 * (HESSIAN * STEP + (G(X0) - G(X))) ***
862 */
863
864 k = temp1;
865 l = g01;
866 i__1 = *p;
867 for (i__ = 1; i__ <= i__1; ++i__) {
868 v[k] = (v[k] - v[l]) / d__[i__];
869 ++k;
870 ++l;
871 /* L450: */
872 }
873
874 /* *** DO GRADIENT TESTS *** */
875
876 if (dv2nrm_(p, &v[temp1]) <= v[1] * v[29]) {
877 goto L460;
878 }
879 if (dd7tpr_(p, &g[1], &v[step1]) >= v[4] * v[30]) {
880 goto L470;
881 }
882 L460:
883 v[16] = v[23];
884
885 /* *** COMPUTE Y VECTOR NEEDED FOR UPDATING S *** */
886
887 L470:
888 dv2axy_(ps, &y[1], &c_b23, &y[1], &g[1]);
889
890 /* *** DETERMINE SIZING FACTOR V(SIZE) *** */
891
892 /* *** SET TEMP1 = S * STEP *** */
893 s1 = iv[62];
894 ds7lvm_(ps, &v[temp1], &v[s1], &v[step1]);
895
896 t1 = (d__1 = dd7tpr_(ps, &v[step1], &v[temp1]), abs(d__1));
897 t = (d__1 = dd7tpr_(ps, &v[step1], &y[1]), abs(d__1));
898 v[55] = 1.;
899 if (t < t1) {
900 v[55] = t / t1;
901 }
902
903 /* *** SET G0 TO WCHMTD CHOICE OF FLETCHER AND AL-BAALI *** */
904
905 hc1 = iv[71];
906 if (hc1 <= 0) {
907 goto L480;
908 }
909 ds7lvm_(ps, &v[g01], &v[hc1], &v[step1]);
910 goto L490;
911
912 L480:
913 rmat1 = iv[78];
914 dl7tvm_(ps, &v[g01], &v[rmat1], &v[step1]);
915 dl7vml_(ps, &v[g01], &v[rmat1], &v[g01]);
916
917 L490:
918 dv2axy_(ps, &v[g01], &c_b43, &y[1], &v[g01]);
919
920 /* *** UPDATE S *** */
921
922 ds7lup_(&v[s1], &v[47], ps, &v[55], &v[step1], &v[temp1], &v[temp2], &v[
923 g01], &v[56], &y[1]);
924 iv[1] = 2;
925 goto L110;
926
927 /* . . . . . . . . . . . . . . MISC. DETAILS . . . . . . . . . . . . . .
928 */
929
930 /* *** BAD PARAMETERS TO ASSESS *** */
931
932 L500:
933 iv[1] = 64;
934 goto L999;
935
936
937 /* *** CONVERGENCE OBTAINED -- SEE WHETHER TO COMPUTE COVARIANCE *** */
938
939 L510:
940 if (iv[57] == 0) {
941 goto L600;
942 }
943 if (iv[74] != 0) {
944 goto L600;
945 }
946 if (iv[55] >= 7) {
947 goto L600;
948 }
949 if (iv[67] > 0) {
950 goto L600;
951 }
952 if (iv[26] > 0) {
953 goto L600;
954 }
955 if (abs(iv[15]) >= 3) {
956 goto L560;
957 }
958 if (iv[9] == 0) {
959 iv[9] = 2;
960 }
961 goto L530;
962
963 /* *** COMPUTE FINITE-DIFFERENCE HESSIAN FOR COMPUTING COVARIANCE ***
964 */
965
966 L520:
967 iv[9] = 0;
968 L530:
969 df7hes_(&d__[1], &g[1], &i__, &iv[1], liv, lv, p, &v[1], &x[1]);
970 switch ((int)i__) {
971 case 1: goto L540;
972 case 2: goto L550;
973 case 3: goto L580;
974 }
975 L540:
976 ++iv[52];
977 ++iv[6];
978 iv[1] = 1;
979 goto L999;
980
981 L550:
982 ++iv[53];
983 ++iv[30];
984 iv[7] = iv[6] + iv[53];
985 iv[1] = 2;
986 goto L999;
987
988 L560:
989 h1 = abs(iv[56]);
990 iv[56] = -h1;
991 pp1o2 = *p * (*p + 1) / 2;
992 rmat1 = iv[78];
993 if (rmat1 <= 0) {
994 goto L570;
995 }
996 lmat1 = iv[42];
997 dv7cpy_(&pp1o2, &v[lmat1], &v[rmat1]);
998 v[53] = 0.;
999 goto L590;
1000 L570:
1001 hc1 = iv[71];
1002 iv[74] = h1;
1003 i__1 = *p * (*p + 1) / 2;
1004 dv7cpy_(&i__1, &v[h1], &v[hc1]);
1005
1006 /* *** COMPUTE CHOLESKY FACTOR OF FINITE-DIFFERENCE HESSIAN */
1007 /* *** FOR USE IN CALLER*S COVARIANCE CALCULATION... */
1008
1009 L580:
1010 lmat1 = iv[42];
1011 h1 = iv[74];
1012 if (h1 <= 0) {
1013 goto L600;
1014 }
1015 if (iv[55] == 70) {
1016 goto L80;
1017 }
1018 dl7srt_(&c__1, p, &v[lmat1], &v[h1], &i__);
1019 iv[74] = -1;
1020 v[53] = 0.;
1021 if (i__ != 0) {
1022 goto L600;
1023 }
1024
1025 L590:
1026 iv[74] = -1;
1027 step1 = iv[40];
1028 t = dl7svn_(p, &v[lmat1], &v[step1], &v[step1]);
1029 if (t <= 0.) {
1030 goto L600;
1031 }
1032 t /= dl7svx_(p, &v[lmat1], &v[step1], &v[step1]);
1033 if (t > dr7mdc_(&c__4)) {
1034 iv[74] = h1;
1035 }
1036 v[53] = t;
1037
1038 L600:
1039 iv[35] = 0;
1040 iv[1] = iv[55];
1041 iv[55] = 0;
1042 goto L999;
1043
1044 /* *** SPECIAL RETURN FOR MISSING HESSIAN INFORMATION -- BOTH */
1045 /* *** IV(HC) .LE. 0 AND IV(RMAT) .LE. 0 */
1046
1047 L610:
1048 iv[1] = 1400;
1049
1050 L999:
1051 return 0;
1052
1053 /* *** LAST LINE OF DG7LIT FOLLOWS *** */
1054 } /* dg7lit_ */
1055
1056