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