cleaning f2c mess

This commit is contained in:
Thomas Capricelli 2009-08-20 23:33:45 +02:00
parent 275a658ec5
commit b1e0662785
5 changed files with 86 additions and 167 deletions

View File

@ -9,7 +9,7 @@ int hybrd_template(minpack_func_nn fcn, void *p, int n, Scalar *x, Scalar *
/* Initialized data */ /* Initialized data */
/* System generated locals */ /* System generated locals */
int fjac_dim1, fjac_offset, i__1, i__2; int fjac_dim1, fjac_offset;
Scalar d__1, d__2; Scalar d__1, d__2;
/* Local variables */ /* Local variables */
@ -59,8 +59,7 @@ int hybrd_template(minpack_func_nn fcn, void *p, int n, Scalar *x, Scalar *
if (mode != 2) { if (mode != 2) {
goto L20; goto L20;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (diag[j] <= 0.) { if (diag[j] <= 0.) {
goto L300; goto L300;
} }
@ -82,8 +81,7 @@ L20:
/* the jacobian matrix. */ /* the jacobian matrix. */
/* Computing MIN */ /* Computing MIN */
i__1 = ml + mu + 1; msum = min(ml + mu + 1, n);
msum = min(i__1,n);
/* initialize iteration counter and monitors. */ /* initialize iteration counter and monitors. */
@ -121,8 +119,7 @@ L30:
if (mode == 2) { if (mode == 2) {
goto L50; goto L50;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
diag[j] = wa2[j]; diag[j] = wa2[j];
if (wa2[j] == 0.) { if (wa2[j] == 0.) {
diag[j] = 1.; diag[j] = 1.;
@ -134,8 +131,7 @@ L50:
/* on the first iteration, calculate the norm of the scaled x */ /* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */ /* and initialize the step bound delta. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa3[j] = diag[j] * x[j]; wa3[j] = diag[j] * x[j];
/* L60: */ /* L60: */
} }
@ -148,25 +144,21 @@ L70:
/* form (q transpose)*fvec and store in qtf. */ /* form (q transpose)*fvec and store in qtf. */
i__1 = n; for (i__ = 1; i__ <= n; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
qtf[i__] = fvec[i__]; qtf[i__] = fvec[i__];
/* L80: */ /* L80: */
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (fjac[j + j * fjac_dim1] == 0.) { if (fjac[j + j * fjac_dim1] == 0.) {
goto L110; goto L110;
} }
sum = 0.; sum = 0.;
i__2 = n; for (i__ = j; i__ <= n; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * qtf[i__]; sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
/* L90: */ /* L90: */
} }
temp = -sum / fjac[j + j * fjac_dim1]; temp = -sum / fjac[j + j * fjac_dim1];
i__2 = n; for (i__ = j; i__ <= n; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
qtf[i__] += fjac[i__ + j * fjac_dim1] * temp; qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
/* L100: */ /* L100: */
} }
@ -178,15 +170,13 @@ L110:
/* copy the triangular factor of the qr factorization into r. */ /* copy the triangular factor of the qr factorization into r. */
sing = FALSE_; sing = FALSE_;
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
l = j; l = j;
jm1 = j - 1; jm1 = j - 1;
if (jm1 < 1) { if (jm1 < 1) {
goto L140; goto L140;
} }
i__2 = jm1; for (i__ = 1; i__ <= jm1; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
r__[l] = fjac[i__ + j * fjac_dim1]; r__[l] = fjac[i__ + j * fjac_dim1];
l = l + n - i__; l = l + n - i__;
/* L130: */ /* L130: */
@ -208,8 +198,7 @@ L140:
if (mode == 2) { if (mode == 2) {
goto L170; goto L170;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
/* Computing MAX */ /* Computing MAX */
d__1 = diag[j], d__2 = wa2[j]; d__1 = diag[j], d__2 = wa2[j];
diag[j] = max(d__1,d__2); diag[j] = max(d__1,d__2);
@ -242,8 +231,7 @@ L190:
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa1[j] = -wa1[j]; wa1[j] = -wa1[j];
wa2[j] = x[j] + wa1[j]; wa2[j] = x[j] + wa1[j];
wa3[j] = diag[j] * wa1[j]; wa3[j] = diag[j] * wa1[j];
@ -278,11 +266,9 @@ L190:
/* compute the scaled predicted reduction. */ /* compute the scaled predicted reduction. */
l = 1; l = 1;
i__1 = n; for (i__ = 1; i__ <= n; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
sum = 0.; sum = 0.;
i__2 = n; for (j = i__; j <= n; ++j) {
for (j = i__; j <= i__2; ++j) {
sum += r__[l] * wa1[j]; sum += r__[l] * wa1[j];
++l; ++l;
/* L210: */ /* L210: */
@ -336,8 +322,7 @@ L240:
/* successful iteration. update x, fvec, and their norms. */ /* successful iteration. update x, fvec, and their norms. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
x[j] = wa2[j]; x[j] = wa2[j];
wa2[j] = diag[j] * x[j]; wa2[j] = diag[j] * x[j];
fvec[j] = wa4[j]; fvec[j] = wa4[j];
@ -400,11 +385,9 @@ L260:
/* calculate the rank one modification to the jacobian */ /* calculate the rank one modification to the jacobian */
/* and update qtf if necessary. */ /* and update qtf if necessary. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
sum = 0.; sum = 0.;
i__2 = n; for (i__ = 1; i__ <= n; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * wa4[i__]; sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
/* L270: */ /* L270: */
} }

View File

@ -10,7 +10,7 @@ int hybrj_template(minpack_funcder_nn fcn, void *p, int n, Scalar *x, Scalar *
/* Initialized data */ /* Initialized data */
/* System generated locals */ /* System generated locals */
int fjac_dim1, fjac_offset, i__1, i__2; int fjac_dim1, fjac_offset;
Scalar d__1, d__2; Scalar d__1, d__2;
/* Local variables */ /* Local variables */
@ -61,8 +61,7 @@ int hybrj_template(minpack_funcder_nn fcn, void *p, int n, Scalar *x, Scalar *
if (mode != 2) { if (mode != 2) {
goto L20; goto L20;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (diag[j] <= 0.) { if (diag[j] <= 0.) {
goto L300; goto L300;
} }
@ -115,8 +114,7 @@ L30:
if (mode == 2) { if (mode == 2) {
goto L50; goto L50;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
diag[j] = wa2[j]; diag[j] = wa2[j];
if (wa2[j] == 0.) { if (wa2[j] == 0.) {
diag[j] = 1.; diag[j] = 1.;
@ -128,8 +126,7 @@ L50:
/* on the first iteration, calculate the norm of the scaled x */ /* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */ /* and initialize the step bound delta. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa3[j] = diag[j] * x[j]; wa3[j] = diag[j] * x[j];
/* L60: */ /* L60: */
} }
@ -142,25 +139,21 @@ L70:
/* form (q transpose)*fvec and store in qtf. */ /* form (q transpose)*fvec and store in qtf. */
i__1 = n; for (i__ = 1; i__ <= n; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
qtf[i__] = fvec[i__]; qtf[i__] = fvec[i__];
/* L80: */ /* L80: */
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (fjac[j + j * fjac_dim1] == 0.) { if (fjac[j + j * fjac_dim1] == 0.) {
goto L110; goto L110;
} }
sum = 0.; sum = 0.;
i__2 = n; for (i__ = j; i__ <= n; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * qtf[i__]; sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
/* L90: */ /* L90: */
} }
temp = -sum / fjac[j + j * fjac_dim1]; temp = -sum / fjac[j + j * fjac_dim1];
i__2 = n; for (i__ = j; i__ <= n; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
qtf[i__] += fjac[i__ + j * fjac_dim1] * temp; qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
/* L100: */ /* L100: */
} }
@ -172,15 +165,13 @@ L110:
/* copy the triangular factor of the qr factorization into r. */ /* copy the triangular factor of the qr factorization into r. */
sing = FALSE_; sing = FALSE_;
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
l = j; l = j;
jm1 = j - 1; jm1 = j - 1;
if (jm1 < 1) { if (jm1 < 1) {
goto L140; goto L140;
} }
i__2 = jm1; for (i__ = 1; i__ <= jm1; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
r__[l] = fjac[i__ + j * fjac_dim1]; r__[l] = fjac[i__ + j * fjac_dim1];
l = l + n - i__; l = l + n - i__;
/* L130: */ /* L130: */
@ -202,8 +193,7 @@ L140:
if (mode == 2) { if (mode == 2) {
goto L170; goto L170;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
/* Computing MAX */ /* Computing MAX */
d__1 = diag[j], d__2 = wa2[j]; d__1 = diag[j], d__2 = wa2[j];
diag[j] = max(d__1,d__2); diag[j] = max(d__1,d__2);
@ -236,8 +226,7 @@ L190:
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa1[j] = -wa1[j]; wa1[j] = -wa1[j];
wa2[j] = x[j] + wa1[j]; wa2[j] = x[j] + wa1[j];
wa3[j] = diag[j] * wa1[j]; wa3[j] = diag[j] * wa1[j];
@ -272,11 +261,9 @@ L190:
/* compute the scaled predicted reduction. */ /* compute the scaled predicted reduction. */
l = 1; l = 1;
i__1 = n; for (i__ = 1; i__ <= n; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
sum = 0.; sum = 0.;
i__2 = n; for (j = i__; j <= n; ++j) {
for (j = i__; j <= i__2; ++j) {
sum += r__[l] * wa1[j]; sum += r__[l] * wa1[j];
++l; ++l;
/* L210: */ /* L210: */
@ -330,8 +317,7 @@ L240:
/* successful iteration. update x, fvec, and their norms. */ /* successful iteration. update x, fvec, and their norms. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
x[j] = wa2[j]; x[j] = wa2[j];
wa2[j] = diag[j] * x[j]; wa2[j] = diag[j] * x[j];
fvec[j] = wa4[j]; fvec[j] = wa4[j];
@ -393,11 +379,9 @@ L260:
/* calculate the rank one modification to the jacobian */ /* calculate the rank one modification to the jacobian */
/* and update qtf if necessary. */ /* and update qtf if necessary. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
sum = 0.; sum = 0.;
i__2 = n; for (i__ = 1; i__ <= n; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * wa4[i__]; sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
/* L270: */ /* L270: */
} }

View File

@ -10,7 +10,7 @@ int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, Scalar *x,
/* Initialized data */ /* Initialized data */
/* System generated locals */ /* System generated locals */
int fjac_dim1, fjac_offset, i__1, i__2; int fjac_dim1, fjac_offset;
Scalar d__1, d__2, d__3; Scalar d__1, d__2, d__3;
/* Local variables */ /* Local variables */
@ -54,8 +54,7 @@ int lmder_template(minpack_funcder_mn fcn, void *p, int m, int n, Scalar *x,
if (mode != 2) { if (mode != 2) {
goto L20; goto L20;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (diag[j] <= 0.) { if (diag[j] <= 0.) {
goto L300; goto L300;
} }
@ -118,8 +117,7 @@ L40:
if (mode == 2) { if (mode == 2) {
goto L60; goto L60;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
diag[j] = wa2[j]; diag[j] = wa2[j];
if (wa2[j] == 0.) { if (wa2[j] == 0.) {
diag[j] = 1.; diag[j] = 1.;
@ -131,8 +129,7 @@ L60:
/* on the first iteration, calculate the norm of the scaled x */ /* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */ /* and initialize the step bound delta. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa3[j] = diag[j] * x[j]; wa3[j] = diag[j] * x[j];
/* L70: */ /* L70: */
} }
@ -146,25 +143,21 @@ L80:
/* form (q transpose)*fvec and store the first n components in */ /* form (q transpose)*fvec and store the first n components in */
/* qtf. */ /* qtf. */
i__1 = m; for (i__ = 1; i__ <= m; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
wa4[i__] = fvec[i__]; wa4[i__] = fvec[i__];
/* L90: */ /* L90: */
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (fjac[j + j * fjac_dim1] == 0.) { if (fjac[j + j * fjac_dim1] == 0.) {
goto L120; goto L120;
} }
sum = 0.; sum = 0.;
i__2 = m; for (i__ = j; i__ <= m; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * wa4[i__]; sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
/* L100: */ /* L100: */
} }
temp = -sum / fjac[j + j * fjac_dim1]; temp = -sum / fjac[j + j * fjac_dim1];
i__2 = m; for (i__ = j; i__ <= m; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
wa4[i__] += fjac[i__ + j * fjac_dim1] * temp; wa4[i__] += fjac[i__ + j * fjac_dim1] * temp;
/* L110: */ /* L110: */
} }
@ -180,15 +173,13 @@ L120:
if (fnorm == 0.) { if (fnorm == 0.) {
goto L170; goto L170;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
l = ipvt[j]; l = ipvt[j];
if (wa2[l] == 0.) { if (wa2[l] == 0.) {
goto L150; goto L150;
} }
sum = 0.; sum = 0.;
i__2 = j; for (i__ = 1; i__ <= j; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm); sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
/* L140: */ /* L140: */
} }
@ -215,8 +206,7 @@ L170:
if (mode == 2) { if (mode == 2) {
goto L190; goto L190;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
/* Computing MAX */ /* Computing MAX */
d__1 = diag[j], d__2 = wa2[j]; d__1 = diag[j], d__2 = wa2[j];
diag[j] = max(d__1,d__2); diag[j] = max(d__1,d__2);
@ -235,8 +225,7 @@ L200:
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa1[j] = -wa1[j]; wa1[j] = -wa1[j];
wa2[j] = x[j] + wa1[j]; wa2[j] = x[j] + wa1[j];
wa3[j] = diag[j] * wa1[j]; wa3[j] = diag[j] * wa1[j];
@ -271,13 +260,11 @@ L200:
/* compute the scaled predicted reduction and */ /* compute the scaled predicted reduction and */
/* the scaled directional derivative. */ /* the scaled directional derivative. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa3[j] = 0.; wa3[j] = 0.;
l = ipvt[j]; l = ipvt[j];
temp = wa1[l]; temp = wa1[l];
i__2 = j; for (i__ = 1; i__ <= j; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
wa3[i__] += fjac[i__ + j * fjac_dim1] * temp; wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
/* L220: */ /* L220: */
} }
@ -340,14 +327,12 @@ L260:
/* successful iteration. update x, fvec, and their norms. */ /* successful iteration. update x, fvec, and their norms. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
x[j] = wa2[j]; x[j] = wa2[j];
wa2[j] = diag[j] * x[j]; wa2[j] = diag[j] * x[j];
/* L270: */ /* L270: */
} }
i__1 = m; for (i__ = 1; i__ <= m; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
fvec[i__] = wa4[i__]; fvec[i__] = wa4[i__];
/* L280: */ /* L280: */
} }

View File

@ -11,7 +11,7 @@ int lmdif_template(minpack_func_mn fcn, void *p, int m, int n, Scalar *x,
/* Initialized data */ /* Initialized data */
/* System generated locals */ /* System generated locals */
int fjac_dim1, fjac_offset, i__1, i__2; int fjac_dim1, fjac_offset;
Scalar d__1, d__2, d__3; Scalar d__1, d__2, d__3;
/* Local variables */ /* Local variables */
@ -55,8 +55,7 @@ int lmdif_template(minpack_func_mn fcn, void *p, int m, int n, Scalar *x,
if (mode != 2) { if (mode != 2) {
goto L20; goto L20;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (diag[j] <= 0.) { if (diag[j] <= 0.) {
goto L300; goto L300;
} }
@ -120,8 +119,7 @@ L40:
if (mode == 2) { if (mode == 2) {
goto L60; goto L60;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
diag[j] = wa2[j]; diag[j] = wa2[j];
if (wa2[j] == 0.) { if (wa2[j] == 0.) {
diag[j] = 1.; diag[j] = 1.;
@ -133,8 +131,7 @@ L60:
/* on the first iteration, calculate the norm of the scaled x */ /* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */ /* and initialize the step bound delta. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa3[j] = diag[j] * x[j]; wa3[j] = diag[j] * x[j];
/* L70: */ /* L70: */
} }
@ -148,25 +145,21 @@ L80:
/* form (q transpose)*fvec and store the first n components in */ /* form (q transpose)*fvec and store the first n components in */
/* qtf. */ /* qtf. */
i__1 = m; for (i__ = 1; i__ <= m; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
wa4[i__] = fvec[i__]; wa4[i__] = fvec[i__];
/* L90: */ /* L90: */
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (fjac[j + j * fjac_dim1] == 0.) { if (fjac[j + j * fjac_dim1] == 0.) {
goto L120; goto L120;
} }
sum = 0.; sum = 0.;
i__2 = m; for (i__ = j; i__ <= m; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * wa4[i__]; sum += fjac[i__ + j * fjac_dim1] * wa4[i__];
/* L100: */ /* L100: */
} }
temp = -sum / fjac[j + j * fjac_dim1]; temp = -sum / fjac[j + j * fjac_dim1];
i__2 = m; for (i__ = j; i__ <= m; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
wa4[i__] += fjac[i__ + j * fjac_dim1] * temp; wa4[i__] += fjac[i__ + j * fjac_dim1] * temp;
/* L110: */ /* L110: */
} }
@ -182,15 +175,13 @@ L120:
if (fnorm == 0.) { if (fnorm == 0.) {
goto L170; goto L170;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
l = ipvt[j]; l = ipvt[j];
if (wa2[l] == 0.) { if (wa2[l] == 0.) {
goto L150; goto L150;
} }
sum = 0.; sum = 0.;
i__2 = j; for (i__ = 1; i__ <= j; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm); sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
/* L140: */ /* L140: */
} }
@ -217,8 +208,7 @@ L170:
if (mode == 2) { if (mode == 2) {
goto L190; goto L190;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
/* Computing MAX */ /* Computing MAX */
d__1 = diag[j], d__2 = wa2[j]; d__1 = diag[j], d__2 = wa2[j];
diag[j] = max(d__1,d__2); diag[j] = max(d__1,d__2);
@ -237,8 +227,7 @@ L200:
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa1[j] = -wa1[j]; wa1[j] = -wa1[j];
wa2[j] = x[j] + wa1[j]; wa2[j] = x[j] + wa1[j];
wa3[j] = diag[j] * wa1[j]; wa3[j] = diag[j] * wa1[j];
@ -273,13 +262,11 @@ L200:
/* compute the scaled predicted reduction and */ /* compute the scaled predicted reduction and */
/* the scaled directional derivative. */ /* the scaled directional derivative. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa3[j] = 0.; wa3[j] = 0.;
l = ipvt[j]; l = ipvt[j];
temp = wa1[l]; temp = wa1[l];
i__2 = j; for (i__ = 1; i__ <= j; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
wa3[i__] += fjac[i__ + j * fjac_dim1] * temp; wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
/* L220: */ /* L220: */
} }
@ -342,14 +329,12 @@ L260:
/* successful iteration. update x, fvec, and their norms. */ /* successful iteration. update x, fvec, and their norms. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
x[j] = wa2[j]; x[j] = wa2[j];
wa2[j] = diag[j] * x[j]; wa2[j] = diag[j] * x[j];
/* L270: */ /* L270: */
} }
i__1 = m; for (i__ = 1; i__ <= m; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
fvec[i__] = wa4[i__]; fvec[i__] = wa4[i__];
/* L280: */ /* L280: */
} }

View File

@ -10,7 +10,7 @@ int lmstr_template(minpack_funcderstr_mn fcn, void *p, int m, int n, Scalar *x,
/* Initialized data */ /* Initialized data */
/* System generated locals */ /* System generated locals */
int fjac_dim1, fjac_offset, i__1, i__2; int fjac_dim1, fjac_offset;
Scalar d__1, d__2, d__3; Scalar d__1, d__2, d__3;
/* Local variables */ /* Local variables */
@ -55,8 +55,7 @@ int lmstr_template(minpack_funcderstr_mn fcn, void *p, int m, int n, Scalar *x,
if (mode != 2) { if (mode != 2) {
goto L20; goto L20;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (diag[j] <= 0.) { if (diag[j] <= 0.) {
goto L340; goto L340;
} }
@ -102,19 +101,16 @@ L40:
/* forming (q transpose)*fvec and storing the first */ /* forming (q transpose)*fvec and storing the first */
/* n components in qtf. */ /* n components in qtf. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
qtf[j] = 0.; qtf[j] = 0.;
i__2 = n; for (i__ = 1; i__ <= n; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
fjac[i__ + j * fjac_dim1] = 0.; fjac[i__ + j * fjac_dim1] = 0.;
/* L50: */ /* L50: */
} }
/* L60: */ /* L60: */
} }
iflag = 2; iflag = 2;
i__1 = m; for (i__ = 1; i__ <= m; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
if ((*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], iflag) < 0) { if ((*fcn)(p, m, n, &x[1], &fvec[1], &wa3[1], iflag) < 0) {
goto L340; goto L340;
} }
@ -130,8 +126,7 @@ L40:
/* reorder its columns and update the components of qtf. */ /* reorder its columns and update the components of qtf. */
sing = FALSE_; sing = FALSE_;
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (fjac[j + j * fjac_dim1] == 0.) { if (fjac[j + j * fjac_dim1] == 0.) {
sing = TRUE_; sing = TRUE_;
} }
@ -144,20 +139,17 @@ L40:
} }
qrfac(n, n, &fjac[fjac_offset], ldfjac, TRUE_, &ipvt[1], n, &wa1[1], & qrfac(n, n, &fjac[fjac_offset], ldfjac, TRUE_, &ipvt[1], n, &wa1[1], &
wa2[1], &wa3[1]); wa2[1], &wa3[1]);
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
if (fjac[j + j * fjac_dim1] == 0.) { if (fjac[j + j * fjac_dim1] == 0.) {
goto L110; goto L110;
} }
sum = 0.; sum = 0.;
i__2 = n; for (i__ = j; i__ <= n; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * qtf[i__]; sum += fjac[i__ + j * fjac_dim1] * qtf[i__];
/* L90: */ /* L90: */
} }
temp = -sum / fjac[j + j * fjac_dim1]; temp = -sum / fjac[j + j * fjac_dim1];
i__2 = n; for (i__ = j; i__ <= n; ++i__) {
for (i__ = j; i__ <= i__2; ++i__) {
qtf[i__] += fjac[i__ + j * fjac_dim1] * temp; qtf[i__] += fjac[i__ + j * fjac_dim1] * temp;
/* L100: */ /* L100: */
} }
@ -176,8 +168,7 @@ L130:
if (mode == 2) { if (mode == 2) {
goto L150; goto L150;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
diag[j] = wa2[j]; diag[j] = wa2[j];
if (wa2[j] == 0.) { if (wa2[j] == 0.) {
diag[j] = 1.; diag[j] = 1.;
@ -189,8 +180,7 @@ L150:
/* on the first iteration, calculate the norm of the scaled x */ /* on the first iteration, calculate the norm of the scaled x */
/* and initialize the step bound delta. */ /* and initialize the step bound delta. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa3[j] = diag[j] * x[j]; wa3[j] = diag[j] * x[j];
/* L160: */ /* L160: */
} }
@ -207,15 +197,13 @@ L170:
if (fnorm == 0.) { if (fnorm == 0.) {
goto L210; goto L210;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
l = ipvt[j]; l = ipvt[j];
if (wa2[l] == 0.) { if (wa2[l] == 0.) {
goto L190; goto L190;
} }
sum = 0.; sum = 0.;
i__2 = j; for (i__ = 1; i__ <= j; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm); sum += fjac[i__ + j * fjac_dim1] * (qtf[i__] / fnorm);
/* L180: */ /* L180: */
} }
@ -242,8 +230,7 @@ L210:
if (mode == 2) { if (mode == 2) {
goto L230; goto L230;
} }
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
/* Computing MAX */ /* Computing MAX */
d__1 = diag[j], d__2 = wa2[j]; d__1 = diag[j], d__2 = wa2[j];
diag[j] = max(d__1,d__2); diag[j] = max(d__1,d__2);
@ -262,8 +249,7 @@ L240:
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa1[j] = -wa1[j]; wa1[j] = -wa1[j];
wa2[j] = x[j] + wa1[j]; wa2[j] = x[j] + wa1[j];
wa3[j] = diag[j] * wa1[j]; wa3[j] = diag[j] * wa1[j];
@ -298,13 +284,11 @@ L240:
/* compute the scaled predicted reduction and */ /* compute the scaled predicted reduction and */
/* the scaled directional derivative. */ /* the scaled directional derivative. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
wa3[j] = 0.; wa3[j] = 0.;
l = ipvt[j]; l = ipvt[j];
temp = wa1[l]; temp = wa1[l];
i__2 = j; for (i__ = 1; i__ <= j; ++i__) {
for (i__ = 1; i__ <= i__2; ++i__) {
wa3[i__] += fjac[i__ + j * fjac_dim1] * temp; wa3[i__] += fjac[i__ + j * fjac_dim1] * temp;
/* L260: */ /* L260: */
} }
@ -367,14 +351,12 @@ L300:
/* successful iteration. update x, fvec, and their norms. */ /* successful iteration. update x, fvec, and their norms. */
i__1 = n; for (j = 1; j <= n; ++j) {
for (j = 1; j <= i__1; ++j) {
x[j] = wa2[j]; x[j] = wa2[j];
wa2[j] = diag[j] * x[j]; wa2[j] = diag[j] * x[j];
/* L310: */ /* L310: */
} }
i__1 = m; for (i__ = 1; i__ <= m; ++i__) {
for (i__ = 1; i__ <= i__1; ++i__) {
fvec[i__] = wa4[i__]; fvec[i__] = wa4[i__];
/* L320: */ /* L320: */
} }