1 /* minpack/lmder.f -- translated by f2c (version 20050501).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #ifdef __cplusplus
14 extern "C" {
15 #endif
16 #include "v3p_netlib.h"
17
18 /* Table of constant values */
19
20 static integer c__1 = 1;
21 static logical c_true = TRUE_;
22
23 /*< >*/
lmder_(void (* fcn)(v3p_netlib_integer *,v3p_netlib_integer *,v3p_netlib_doublereal *,v3p_netlib_doublereal *,v3p_netlib_doublereal *,v3p_netlib_integer *,v3p_netlib_integer *,void *),integer * m,integer * n,doublereal * x,doublereal * fvec,doublereal * fjac,integer * ldfjac,doublereal * ftol,doublereal * xtol,doublereal * gtol,integer * maxfev,doublereal * diag,integer * mode,doublereal * factor,integer * nprint,integer * info,integer * nfev,integer * njev,integer * ipvt,doublereal * qtf,doublereal * wa1,doublereal * wa2,doublereal * wa3,doublereal * wa4,void * userdata)24 /* Subroutine */ int lmder_(
25 void (*fcn)(v3p_netlib_integer*,
26 v3p_netlib_integer*,
27 v3p_netlib_doublereal*,
28 v3p_netlib_doublereal*,
29 v3p_netlib_doublereal*,
30 v3p_netlib_integer*,
31 v3p_netlib_integer*,
32 void*),
33 integer *m, integer *n, doublereal *x,
34 doublereal *fvec, doublereal *fjac, integer *ldfjac, doublereal *ftol,
35 doublereal *xtol, doublereal *gtol, integer *maxfev, doublereal *
36 diag, integer *mode, doublereal *factor, integer *nprint, integer *
37 info, integer *nfev, integer *njev, integer *ipvt, doublereal *qtf,
38 doublereal *wa1, doublereal *wa2, doublereal *wa3, doublereal *wa4,
39 void* userdata)
40 {
41 /* Initialized data */
42
43 static doublereal one = 1.; /* constant */
44 static doublereal p1 = .1; /* constant */
45 static doublereal p5 = .5; /* constant */
46 static doublereal p25 = .25; /* constant */
47 static doublereal p75 = .75; /* constant */
48 static doublereal p0001 = 1e-4; /* constant */
49 static doublereal zero = 0.; /* constant */
50
51 /* System generated locals */
52 integer fjac_dim1, fjac_offset, i__1, i__2;
53 doublereal d__1, d__2, d__3;
54
55 /* Builtin functions */
56 double sqrt(doublereal);
57
58 /* Local variables */
59 integer i__, j, l;
60 doublereal par, sum;
61 integer iter;
62 doublereal temp=0, temp1, temp2;
63 integer iflag;
64 doublereal delta;
65 extern /* Subroutine */ int qrfac_(integer *, integer *, doublereal *,
66 integer *, logical *, integer *, integer *, doublereal *,
67 doublereal *, doublereal *), lmpar_(integer *, doublereal *,
68 integer *, integer *, doublereal *, doublereal *, doublereal *,
69 doublereal *, doublereal *, doublereal *, doublereal *,
70 doublereal *);
71 doublereal ratio;
72 extern doublereal enorm_(integer *, doublereal *);
73 doublereal fnorm, gnorm, pnorm, xnorm=0, fnorm1, actred, dirder, epsmch,
74 prered;
75 extern doublereal dpmpar_(integer *);
76
77 /*< integer m,n,ldfjac,maxfev,mode,nprint,info,nfev,njev >*/
78 /*< integer ipvt(n) >*/
79 /*< double precision ftol,xtol,gtol,factor >*/
80 /*< >*/
81 /* ********** */
82
83 /* subroutine lmder */
84
85 /* the purpose of lmder is to minimize the sum of the squares of */
86 /* m nonlinear functions in n variables by a modification of */
87 /* the levenberg-marquardt algorithm. the user must provide a */
88 /* subroutine which calculates the functions and the jacobian. */
89
90 /* the subroutine statement is */
91
92 /* subroutine lmder(fcn,m,n,x,fvec,fjac,ldfjac,ftol,xtol,gtol, */
93 /* maxfev,diag,mode,factor,nprint,info,nfev, */
94 /* njev,ipvt,qtf,wa1,wa2,wa3,wa4) */
95
96 /* where */
97
98 /* fcn is the name of the user-supplied subroutine which */
99 /* calculates the functions and the jacobian. fcn must */
100 /* be declared in an external statement in the user */
101 /* calling program, and should be written as follows. */
102
103 /* subroutine fcn(m,n,x,fvec,fjac,ldfjac,iflag) */
104 /* integer m,n,ldfjac,iflag */
105 /* double precision x(n),fvec(m),fjac(ldfjac,n) */
106 /* ---------- */
107 /* if iflag = 1 calculate the functions at x and */
108 /* return this vector in fvec. do not alter fjac. */
109 /* if iflag = 2 calculate the jacobian at x and */
110 /* return this matrix in fjac. do not alter fvec. */
111 /* ---------- */
112 /* return */
113 /* end */
114
115 /* the value of iflag should not be changed by fcn unless */
116 /* the user wants to terminate execution of lmder. */
117 /* in this case set iflag to a negative integer. */
118
119 /* m is a positive integer input variable set to the number */
120 /* of functions. */
121
122 /* n is a positive integer input variable set to the number */
123 /* of variables. n must not exceed m. */
124
125 /* x is an array of length n. on input x must contain */
126 /* an initial estimate of the solution vector. on output x */
127 /* contains the final estimate of the solution vector. */
128
129 /* fvec is an output array of length m which contains */
130 /* the functions evaluated at the output x. */
131
132 /* fjac is an output m by n array. the upper n by n submatrix */
133 /* of fjac contains an upper triangular matrix r with */
134 /* diagonal elements of nonincreasing magnitude such that */
135
136 /* t t t */
137 /* p *(jac *jac)*p = r *r, */
138
139 /* where p is a permutation matrix and jac is the final */
140 /* calculated jacobian. column j of p is column ipvt(j) */
141 /* (see below) of the identity matrix. the lower trapezoidal */
142 /* part of fjac contains information generated during */
143 /* the computation of r. */
144
145 /* ldfjac is a positive integer input variable not less than m */
146 /* which specifies the leading dimension of the array fjac. */
147
148 /* ftol is a nonnegative input variable. termination */
149 /* occurs when both the actual and predicted relative */
150 /* reductions in the sum of squares are at most ftol. */
151 /* therefore, ftol measures the relative error desired */
152 /* in the sum of squares. */
153
154 /* xtol is a nonnegative input variable. termination */
155 /* occurs when the relative error between two consecutive */
156 /* iterates is at most xtol. therefore, xtol measures the */
157 /* relative error desired in the approximate solution. */
158
159 /* gtol is a nonnegative input variable. termination */
160 /* occurs when the cosine of the angle between fvec and */
161 /* any column of the jacobian is at most gtol in absolute */
162 /* value. therefore, gtol measures the orthogonality */
163 /* desired between the function vector and the columns */
164 /* of the jacobian. */
165
166 /* maxfev is a positive integer input variable. termination */
167 /* occurs when the number of calls to fcn with iflag = 1 */
168 /* has reached maxfev. */
169
170 /* diag is an array of length n. if mode = 1 (see */
171 /* below), diag is internally set. if mode = 2, diag */
172 /* must contain positive entries that serve as */
173 /* multiplicative scale factors for the variables. */
174
175 /* mode is an integer input variable. if mode = 1, the */
176 /* variables will be scaled internally. if mode = 2, */
177 /* the scaling is specified by the input diag. other */
178 /* values of mode are equivalent to mode = 1. */
179
180 /* factor is a positive input variable used in determining the */
181 /* initial step bound. this bound is set to the product of */
182 /* factor and the euclidean norm of diag*x if nonzero, or else */
183 /* to factor itself. in most cases factor should lie in the */
184 /* interval (.1,100.).100. is a generally recommended value. */
185
186 /* nprint is an integer input variable that enables controlled */
187 /* printing of iterates if it is positive. in this case, */
188 /* fcn is called with iflag = 0 at the beginning of the first */
189 /* iteration and every nprint iterations thereafter and */
190 /* immediately prior to return, with x, fvec, and fjac */
191 /* available for printing. fvec and fjac should not be */
192 /* altered. if nprint is not positive, no special calls */
193 /* of fcn with iflag = 0 are made. */
194
195 /* info is an integer output variable. if the user has */
196 /* terminated execution, info is set to the (negative) */
197 /* value of iflag. see description of fcn. otherwise, */
198 /* info is set as follows. */
199
200 /* info = 0 improper input parameters. */
201
202 /* info = 1 both actual and predicted relative reductions */
203 /* in the sum of squares are at most ftol. */
204
205 /* info = 2 relative error between two consecutive iterates */
206 /* is at most xtol. */
207
208 /* info = 3 conditions for info = 1 and info = 2 both hold. */
209
210 /* info = 4 the cosine of the angle between fvec and any */
211 /* column of the jacobian is at most gtol in */
212 /* absolute value. */
213
214 /* info = 5 number of calls to fcn with iflag = 1 has */
215 /* reached maxfev. */
216
217 /* info = 6 ftol is too small. no further reduction in */
218 /* the sum of squares is possible. */
219
220 /* info = 7 xtol is too small. no further improvement in */
221 /* the approximate solution x is possible. */
222
223 /* info = 8 gtol is too small. fvec is orthogonal to the */
224 /* columns of the jacobian to machine precision. */
225
226 /* nfev is an integer output variable set to the number of */
227 /* calls to fcn with iflag = 1. */
228
229 /* njev is an integer output variable set to the number of */
230 /* calls to fcn with iflag = 2. */
231
232 /* ipvt is an integer output array of length n. ipvt */
233 /* defines a permutation matrix p such that jac*p = q*r, */
234 /* where jac is the final calculated jacobian, q is */
235 /* orthogonal (not stored), and r is upper triangular */
236 /* with diagonal elements of nonincreasing magnitude. */
237 /* column j of p is column ipvt(j) of the identity matrix. */
238
239 /* qtf is an output array of length n which contains */
240 /* the first n elements of the vector (q transpose)*fvec. */
241
242 /* wa1, wa2, and wa3 are work arrays of length n. */
243
244 /* wa4 is a work array of length m. */
245
246 /* subprograms called */
247
248 /* user-supplied ...... fcn */
249
250 /* minpack-supplied ... dpmpar,enorm,lmpar,qrfac */
251
252 /* fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod */
253
254 /* argonne national laboratory. minpack project. march 1980. */
255 /* burton s. garbow, kenneth e. hillstrom, jorge j. more */
256
257 /* ********** */
258 /*< integer i,iflag,iter,j,l >*/
259 /*< >*/
260 /*< double precision dpmpar,enorm >*/
261 /*< >*/
262 /* Parameter adjustments */
263 --wa4;
264 --fvec;
265 --wa3;
266 --wa2;
267 --wa1;
268 --qtf;
269 --ipvt;
270 --diag;
271 --x;
272 fjac_dim1 = *ldfjac;
273 fjac_offset = 1 + fjac_dim1;
274 fjac -= fjac_offset;
275
276 /* Function Body */
277
278 /* epsmch is the machine precision. */
279
280 /*< epsmch = dpmpar(1) >*/
281 epsmch = dpmpar_(&c__1);
282
283 /*< info = 0 >*/
284 *info = 0;
285 /*< iflag = 0 >*/
286 iflag = 0;
287 /*< nfev = 0 >*/
288 *nfev = 0;
289 /*< njev = 0 >*/
290 *njev = 0;
291
292 /* check the input parameters for errors. */
293
294 /*< >*/
295 if (*n <= 0 || *m < *n || *ldfjac < *m || *ftol < zero || *xtol < zero ||
296 *gtol < zero || *maxfev <= 0 || *factor <= zero) {
297 goto L300;
298 }
299 /*< if (mode .ne. 2) go to 20 >*/
300 if (*mode != 2) {
301 goto L20;
302 }
303 /*< do 10 j = 1, n >*/
304 i__1 = *n;
305 for (j = 1; j <= i__1; ++j) {
306 /*< if (diag(j) .le. zero) go to 300 >*/
307 if (diag[j] <= zero) {
308 goto L300;
309 }
310 /*< 10 continue >*/
311 /* L10: */
312 }
313 /*< 20 continue >*/
314 L20:
315
316 /* evaluate the function at the starting point */
317 /* and calculate its norm. */
318
319 /*< iflag = 1 >*/
320 iflag = 1;
321 /*< call fcn(m,n,x,fvec,fjac,ldfjac,iflag) >*/
322 (*fcn)(m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, &iflag,
323 userdata);
324 /*< nfev = 1 >*/
325 *nfev = 1;
326 /*< if (iflag .lt. 0) go to 300 >*/
327 if (iflag < 0) {
328 goto L300;
329 }
330 /*< fnorm = enorm(m,fvec) >*/
331 fnorm = enorm_(m, &fvec[1]);
332
333 /* initialize levenberg-marquardt parameter and iteration counter. */
334
335 /*< par = zero >*/
336 par = zero;
337 /*< iter = 1 >*/
338 iter = 1;
339
340 /* beginning of the outer loop. */
341
342 /*< 30 continue >*/
343 L30:
344
345 /* calculate the jacobian matrix. */
346
347 /*< iflag = 2 >*/
348 iflag = 2;
349 /*< call fcn(m,n,x,fvec,fjac,ldfjac,iflag) >*/
350 (*fcn)(m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, &iflag,
351 userdata);
352 /*< njev = njev + 1 >*/
353 ++(*njev);
354 /*< if (iflag .lt. 0) go to 300 >*/
355 if (iflag < 0) {
356 goto L300;
357 }
358
359 /* if requested, call fcn to enable printing of iterates. */
360
361 /*< if (nprint .le. 0) go to 40 >*/
362 if (*nprint <= 0) {
363 goto L40;
364 }
365 /*< iflag = 0 >*/
366 iflag = 0;
367 /*< >*/
368 if ((iter - 1) % *nprint == 0) {
369 (*fcn)(m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, &iflag,
370 userdata);
371 }
372 /*< if (iflag .lt. 0) go to 300 >*/
373 if (iflag < 0) {
374 goto L300;
375 }
376 /*< 40 continue >*/
377 L40:
378
379 /* compute the qr factorization of the jacobian. */
380
381 /*< call qrfac(m,n,fjac,ldfjac,.true.,ipvt,n,wa1,wa2,wa3) >*/
382 qrfac_(m, n, &fjac[fjac_offset], ldfjac, &c_true, &ipvt[1], n, &wa1[1], &
383 wa2[1], &wa3[1]);
384
385 /* on the first iteration and if mode is 1, scale according */
386 /* to the norms of the columns of the initial jacobian. */
387
388 /*< if (iter .ne. 1) go to 80 >*/
389 if (iter != 1) {
390 goto L80;
391 }
392 /*< if (mode .eq. 2) go to 60 >*/
393 if (*mode == 2) {
394 goto L60;
395 }
396 /*< do 50 j = 1, n >*/
397 i__1 = *n;
398 for (j = 1; j <= i__1; ++j) {
399 /*< diag(j) = wa2(j) >*/
400 diag[j] = wa2[j];
401 /*< if (wa2(j) .eq. zero) diag(j) = one >*/
402 if (wa2[j] == zero) {
403 diag[j] = one;
404 }
405 /*< 50 continue >*/
406 /* L50: */
407 }
408 /*< 60 continue >*/
409 L60:
410
411 /* on the first iteration, calculate the norm of the scaled x */
412 /* and initialize the step bound delta. */
413
414 /*< do 70 j = 1, n >*/
415 i__1 = *n;
416 for (j = 1; j <= i__1; ++j) {
417 /*< wa3(j) = diag(j)*x(j) >*/
418 wa3[j] = diag[j] * x[j];
419 /*< 70 continue >*/
420 /* L70: */
421 }
422 /*< xnorm = enorm(n,wa3) >*/
423 xnorm = enorm_(n, &wa3[1]);
424 /*< delta = factor*xnorm >*/
425 delta = *factor * xnorm;
426 /*< if (delta .eq. zero) delta = factor >*/
427 if (delta == zero) {
428 delta = *factor;
429 }
430 /*< 80 continue >*/
431 L80:
432
433 /* form (q transpose)*fvec and store the first n components in */
434 /* qtf. */
435
436 /*< do 90 i = 1, m >*/
437 i__1 = *m;
438 for (i__ = 1; i__ <= i__1; ++i__) {
439 /*< wa4(i) = fvec(i) >*/
440 wa4[i__] = fvec[i__];
441 /*< 90 continue >*/
442 /* L90: */
443 }
444 /*< do 130 j = 1, n >*/
445 i__1 = *n;
446 for (j = 1; j <= i__1; ++j) {
447 /*< if (fjac(j,j) .eq. zero) go to 120 >*/
448 if (fjac[j + j * fjac_dim1] == zero) {
449 goto L120;
450 }
451 /*< sum = zero >*/
452 sum = zero;
453 /*< do 100 i = j, m >*/
454 i__2 = *m;
455 for (i__ = j; i__ <= i__2; ++i__) {
456 /*< sum = sum + fjac(i,j)*wa4(i) >*/
457 sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
458 /*< 100 continue >*/
459 /* L100: */
460 }
461 /*< temp = -sum/fjac(j,j) >*/
462 temp = -sum / fjac[j + j * fjac_dim1];
463 /*< do 110 i = j, m >*/
464 i__2 = *m;
465 for (i__ = j; i__ <= i__2; ++i__) {
466 /*< wa4(i) = wa4(i) + fjac(i,j)*temp >*/
467 wa4[i__] += fjac[i__ + j * fjac_dim1] * temp;
468 /*< 110 continue >*/
469 /* L110: */
470 }
471 /*< 120 continue >*/
472 L120:
473 /*< fjac(j,j) = wa1(j) >*/
474 fjac[j + j * fjac_dim1] = wa1[j];
475 /*< qtf(j) = wa4(j) >*/
476 qtf[j] = wa4[j];
477 /*< 130 continue >*/
478 /* L130: */
479 }
480
481 /* compute the norm of the scaled gradient. */
482
483 /*< gnorm = zero >*/
484 gnorm = zero;
485 /*< if (fnorm .eq. zero) go to 170 >*/
486 if (fnorm == zero) {
487 goto L170;
488 }
489 /*< do 160 j = 1, n >*/
490 i__1 = *n;
491 for (j = 1; j <= i__1; ++j) {
492 /*< l = ipvt(j) >*/
493 l = ipvt[j];
494 /*< if (wa2(l) .eq. zero) go to 150 >*/
495 if (wa2[l] == zero) {
496 goto L150;
497 }
498 /*< sum = zero >*/
499 sum = zero;
500 /*< do 140 i = 1, j >*/
501 i__2 = j;
502 for (i__ = 1; i__ <= i__2; ++i__) {
503 /*< sum = sum + fjac(i,j)*(qtf(i)/fnorm) >*/
504 sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
505 /*< 140 continue >*/
506 /* L140: */
507 }
508 /*< gnorm = dmax1(gnorm,dabs(sum/wa2(l))) >*/
509 /* Computing MAX */
510 d__2 = gnorm, d__3 = (d__1 = sum / wa2[l], abs(d__1));
511 gnorm = max(d__2,d__3);
512 /*< 150 continue >*/
513 L150:
514 /*< 160 continue >*/
515 /* L160: */
516 ;
517 }
518 /*< 170 continue >*/
519 L170:
520
521 /* test for convergence of the gradient norm. */
522
523 /*< if (gnorm .le. gtol) info = 4 >*/
524 if (gnorm <= *gtol) {
525 *info = 4;
526 }
527 /*< if (info .ne. 0) go to 300 >*/
528 if (*info != 0) {
529 goto L300;
530 }
531
532 /* rescale if necessary. */
533
534 /*< if (mode .eq. 2) go to 190 >*/
535 if (*mode == 2) {
536 goto L190;
537 }
538 /*< do 180 j = 1, n >*/
539 i__1 = *n;
540 for (j = 1; j <= i__1; ++j) {
541 /*< diag(j) = dmax1(diag(j),wa2(j)) >*/
542 /* Computing MAX */
543 d__1 = diag[j], d__2 = wa2[j];
544 diag[j] = max(d__1,d__2);
545 /*< 180 continue >*/
546 /* L180: */
547 }
548 /*< 190 continue >*/
549 L190:
550
551 /* beginning of the inner loop. */
552
553 /*< 200 continue >*/
554 L200:
555
556 /* determine the levenberg-marquardt parameter. */
557
558 /*< >*/
559 lmpar_(n, &fjac[fjac_offset], ldfjac, &ipvt[1], &diag[1], &qtf[1], &delta,
560 &par, &wa1[1], &wa2[1], &wa3[1], &wa4[1]);
561
562 /* store the direction p and x + p. calculate the norm of p. */
563
564 /*< do 210 j = 1, n >*/
565 i__1 = *n;
566 for (j = 1; j <= i__1; ++j) {
567 /*< wa1(j) = -wa1(j) >*/
568 wa1[j] = -wa1[j];
569 /*< wa2(j) = x(j) + wa1(j) >*/
570 wa2[j] = x[j] + wa1[j];
571 /*< wa3(j) = diag(j)*wa1(j) >*/
572 wa3[j] = diag[j] * wa1[j];
573 /*< 210 continue >*/
574 /* L210: */
575 }
576 /*< pnorm = enorm(n,wa3) >*/
577 pnorm = enorm_(n, &wa3[1]);
578
579 /* on the first iteration, adjust the initial step bound. */
580
581 /*< if (iter .eq. 1) delta = dmin1(delta,pnorm) >*/
582 if (iter == 1) {
583 delta = min(delta,pnorm);
584 }
585
586 /* evaluate the function at x + p and calculate its norm. */
587
588 /*< iflag = 1 >*/
589 iflag = 1;
590 /*< call fcn(m,n,wa2,wa4,fjac,ldfjac,iflag) >*/
591 (*fcn)(m, n, &wa2[1], &wa4[1], &fjac[fjac_offset], ldfjac, &iflag,
592 userdata);
593 /*< nfev = nfev + 1 >*/
594 ++(*nfev);
595 /*< if (iflag .lt. 0) go to 300 >*/
596 if (iflag < 0) {
597 goto L300;
598 }
599 /*< fnorm1 = enorm(m,wa4) >*/
600 fnorm1 = enorm_(m, &wa4[1]);
601
602 /* compute the scaled actual reduction. */
603
604 /*< actred = -one >*/
605 actred = -one;
606 /*< if (p1*fnorm1 .lt. fnorm) actred = one - (fnorm1/fnorm)**2 >*/
607 if (p1 * fnorm1 < fnorm) {
608 /* Computing 2nd power */
609 d__1 = fnorm1 / fnorm;
610 actred = one - d__1 * d__1;
611 }
612
613 /* compute the scaled predicted reduction and */
614 /* the scaled directional derivative. */
615
616 /*< do 230 j = 1, n >*/
617 i__1 = *n;
618 for (j = 1; j <= i__1; ++j) {
619 /*< wa3(j) = zero >*/
620 wa3[j] = zero;
621 /*< l = ipvt(j) >*/
622 l = ipvt[j];
623 /*< temp = wa1(l) >*/
624 temp = wa1[l];
625 /*< do 220 i = 1, j >*/
626 i__2 = j;
627 for (i__ = 1; i__ <= i__2; ++i__) {
628 /*< wa3(i) = wa3(i) + fjac(i,j)*temp >*/
629 wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
630 /*< 220 continue >*/
631 /* L220: */
632 }
633 /*< 230 continue >*/
634 /* L230: */
635 }
636 /*< temp1 = enorm(n,wa3)/fnorm >*/
637 temp1 = enorm_(n, &wa3[1]) / fnorm;
638 /*< temp2 = (dsqrt(par)*pnorm)/fnorm >*/
639 temp2 = sqrt(par) * pnorm / fnorm;
640 /*< prered = temp1**2 + temp2**2/p5 >*/
641 /* Computing 2nd power */
642 d__1 = temp1;
643 /* Computing 2nd power */
644 d__2 = temp2;
645 prered = d__1 * d__1 + d__2 * d__2 / p5;
646 /*< dirder = -(temp1**2 + temp2**2) >*/
647 /* Computing 2nd power */
648 d__1 = temp1;
649 /* Computing 2nd power */
650 d__2 = temp2;
651 dirder = -(d__1 * d__1 + d__2 * d__2);
652
653 /* compute the ratio of the actual to the predicted */
654 /* reduction. */
655
656 /*< ratio = zero >*/
657 ratio = zero;
658 /*< if (prered .ne. zero) ratio = actred/prered >*/
659 if (prered != zero) {
660 ratio = actred / prered;
661 }
662
663 /* update the step bound. */
664
665 /*< if (ratio .gt. p25) go to 240 >*/
666 if (ratio > p25) {
667 goto L240;
668 }
669 /*< if (actred .ge. zero) temp = p5 >*/
670 if (actred >= zero) {
671 temp = p5;
672 }
673 /*< >*/
674 if (actred < zero) {
675 temp = p5 * dirder / (dirder + p5 * actred);
676 }
677 /*< if (p1*fnorm1 .ge. fnorm .or. temp .lt. p1) temp = p1 >*/
678 if (p1 * fnorm1 >= fnorm || temp < p1) {
679 temp = p1;
680 }
681 /*< delta = temp*dmin1(delta,pnorm/p1) >*/
682 /* Computing MIN */
683 d__1 = delta, d__2 = pnorm / p1;
684 delta = temp * min(d__1,d__2);
685 /*< par = par/temp >*/
686 par /= temp;
687 /*< go to 260 >*/
688 goto L260;
689 /*< 240 continue >*/
690 L240:
691 /*< if (par .ne. zero .and. ratio .lt. p75) go to 250 >*/
692 if (par != zero && ratio < p75) {
693 goto L250;
694 }
695 /*< delta = pnorm/p5 >*/
696 delta = pnorm / p5;
697 /*< par = p5*par >*/
698 par = p5 * par;
699 /*< 250 continue >*/
700 L250:
701 /*< 260 continue >*/
702 L260:
703
704 /* test for successful iteration. */
705
706 /*< if (ratio .lt. p0001) go to 290 >*/
707 if (ratio < p0001) {
708 goto L290;
709 }
710
711 /* successful iteration. update x, fvec, and their norms. */
712
713 /*< do 270 j = 1, n >*/
714 i__1 = *n;
715 for (j = 1; j <= i__1; ++j) {
716 /*< x(j) = wa2(j) >*/
717 x[j] = wa2[j];
718 /*< wa2(j) = diag(j)*x(j) >*/
719 wa2[j] = diag[j] * x[j];
720 /*< 270 continue >*/
721 /* L270: */
722 }
723 /*< do 280 i = 1, m >*/
724 i__1 = *m;
725 for (i__ = 1; i__ <= i__1; ++i__) {
726 /*< fvec(i) = wa4(i) >*/
727 fvec[i__] = wa4[i__];
728 /*< 280 continue >*/
729 /* L280: */
730 }
731 /*< xnorm = enorm(n,wa2) >*/
732 xnorm = enorm_(n, &wa2[1]);
733 /*< fnorm = fnorm1 >*/
734 fnorm = fnorm1;
735 /*< iter = iter + 1 >*/
736 ++iter;
737 /*< 290 continue >*/
738 L290:
739
740 /* tests for convergence. */
741
742 /*< >*/
743 if (abs(actred) <= *ftol && prered <= *ftol && p5 * ratio <= one) {
744 *info = 1;
745 }
746 /*< if (delta .le. xtol*xnorm) info = 2 >*/
747 if (delta <= *xtol * xnorm) {
748 *info = 2;
749 }
750 /*< >*/
751 if (abs(actred) <= *ftol && prered <= *ftol && p5 * ratio <= one && *info
752 == 2) {
753 *info = 3;
754 }
755 /*< if (info .ne. 0) go to 300 >*/
756 if (*info != 0) {
757 goto L300;
758 }
759
760 /* tests for termination and stringent tolerances. */
761
762 /*< if (nfev .ge. maxfev) info = 5 >*/
763 if (*nfev >= *maxfev) {
764 *info = 5;
765 }
766 /*< >*/
767 if (abs(actred) <= epsmch && prered <= epsmch && p5 * ratio <= one) {
768 *info = 6;
769 }
770 /*< if (delta .le. epsmch*xnorm) info = 7 >*/
771 if (delta <= epsmch * xnorm) {
772 *info = 7;
773 }
774 /*< if (gnorm .le. epsmch) info = 8 >*/
775 if (gnorm <= epsmch) {
776 *info = 8;
777 }
778 /*< if (info .ne. 0) go to 300 >*/
779 if (*info != 0) {
780 goto L300;
781 }
782
783 /* end of the inner loop. repeat if iteration unsuccessful. */
784
785 /*< if (ratio .lt. p0001) go to 200 >*/
786 if (ratio < p0001) {
787 goto L200;
788 }
789
790 /* end of the outer loop. */
791
792 /*< go to 30 >*/
793 goto L30;
794 /*< 300 continue >*/
795 L300:
796
797 /* termination, either normal or user imposed. */
798
799 /*< if (iflag .lt. 0) info = iflag >*/
800 if (iflag < 0) {
801 *info = iflag;
802 }
803 /*< iflag = 0 >*/
804 iflag = 0;
805 /*< if (nprint .gt. 0) call fcn(m,n,x,fvec,fjac,ldfjac,iflag) >*/
806 if (*nprint > 0) {
807 (*fcn)(m, n, &x[1], &fvec[1], &fjac[fjac_offset], ldfjac, &iflag,
808 userdata);
809 }
810 /*< return >*/
811 return 0;
812
813 /* last card of subroutine lmder. */
814
815 /*< end >*/
816 } /* lmder_ */
817
818 #ifdef __cplusplus
819 }
820 #endif
821