cleaning : removing #define, use std:min() and such

This commit is contained in:
Thomas Capricelli 2009-08-22 05:29:33 +02:00
parent 20480a5438
commit 11c3762068
6 changed files with 104 additions and 116 deletions

View File

@ -36,18 +36,6 @@ namespace Eigen {
*/ */
//@{ //@{
#define min(a,b) ((a) <= (b) ? (a) : (b))
#define max(a,b) ((a) >= (b) ? (a) : (b))
#define abs(x) ((x) >= 0 ? (x) : -(x))
#define TRUE_ (1)
#define FALSE_ (0)
#define p1 .1
#define p5 .5
#define p25 .25
#define p75 .75
#define p001 .001
#define p0001 1e-4
#if 1 #if 1
#include <cminpack.h> #include <cminpack.h>
#else #else

View File

@ -86,7 +86,7 @@ L20:
/* the jacobian matrix. */ /* the jacobian matrix. */
/* Computing MIN */ /* Computing MIN */
msum = min(nb_of_subdiagonals + nb_of_superdiagonals + 1, n); msum = std::min(nb_of_subdiagonals + nb_of_superdiagonals + 1, n);
/* initialize iteration counter and monitors. */ /* initialize iteration counter and monitors. */
@ -99,7 +99,7 @@ L20:
/* beginning of the outer loop. */ /* beginning of the outer loop. */
L30: L30:
jeval = TRUE_; jeval = true;
/* calculate the jacobian matrix. */ /* calculate the jacobian matrix. */
@ -112,7 +112,7 @@ L30:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
qrfac(n, n, fjac.data(), ldfjac, FALSE_, iwa, 1, wa1.data(), wa2.data(), wa3.data()); qrfac(n, n, fjac.data(), ldfjac, false, iwa, 1, wa1.data(), wa2.data(), wa3.data());
/* on the first iteration and if mode is 1, scale according */ /* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */ /* to the norms of the columns of the initial jacobian. */
@ -173,7 +173,7 @@ 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;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
l = j; l = j;
if (j) { if (j) {
@ -185,7 +185,7 @@ L110:
} }
R[l] = wa1[j]; R[l] = wa1[j];
if (wa1[j] == 0.) { if (wa1[j] == 0.) {
sing = TRUE_; sing = true;
} }
/* L150: */ /* L150: */
} }
@ -201,7 +201,7 @@ L110:
} }
/* Computing MAX */ /* Computing MAX */
for (j = 0; j < n; ++j) for (j = 0; j < n; ++j)
diag[j] = max(diag[j], wa2[j]); diag[j] = std::max(diag[j], wa2[j]);
L170: L170:
/* beginning of the inner loop. */ /* beginning of the inner loop. */
@ -224,7 +224,7 @@ L190:
/* determine the direction p. */ /* determine the direction p. */
dogleg(n, R.data(), lr, diag.data(), qtf.data(), delta, wa1.data(), wa2.data(), wa3.data()); ei_dogleg<Scalar>(n, R.data(), lr, diag.data(), qtf.data(), delta, wa1.data(), wa2.data(), wa3.data());
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
@ -239,7 +239,7 @@ L190:
/* on the first iteration, adjust the initial step bound. */ /* on the first iteration, adjust the initial step bound. */
if (iter == 1) { if (iter == 1) {
delta = min(delta,pnorm); delta = std::min(delta,pnorm);
} }
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
@ -285,26 +285,26 @@ L190:
/* update the step bound. */ /* update the step bound. */
if (ratio >= p1) { if (ratio >= Scalar(.1)) {
goto L230; goto L230;
} }
ncsuc = 0; ncsuc = 0;
++ncfail; ++ncfail;
delta = p5 * delta; delta = Scalar(.5) * delta;
goto L240; goto L240;
L230: L230:
ncfail = 0; ncfail = 0;
++ncsuc; ++ncsuc;
if (ratio >= p5 || ncsuc > 1) /* Computing MAX */ if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
delta = max(delta, pnorm / p5); delta = std::max(delta, pnorm / Scalar(.5));
if (ei_abs(ratio - 1.) <= p1) { if (ei_abs(ratio - 1.) <= Scalar(.1)) {
delta = pnorm / p5; delta = pnorm / Scalar(.5);
} }
L240: L240:
/* test for successful iteration. */ /* test for successful iteration. */
if (ratio < p0001) { if (ratio < Scalar(1e-4)) {
goto L260; goto L260;
} }
@ -324,13 +324,13 @@ L260:
/* determine the progress of the iteration. */ /* determine the progress of the iteration. */
++nslow1; ++nslow1;
if (actred >= p001) { if (actred >= Scalar(.001)) {
nslow1 = 0; nslow1 = 0;
} }
if (jeval) { if (jeval) {
++nslow2; ++nslow2;
} }
if (actred >= p1) { if (actred >= Scalar(.1)) {
nslow2 = 0; nslow2 = 0;
} }
@ -349,7 +349,7 @@ L260:
info = 2; info = 2;
} }
/* Computing MAX */ /* Computing MAX */
if (p1 * max(p1 * delta, pnorm) <= epsilon<Scalar>() * xnorm) if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm)
info = 3; info = 3;
if (nslow2 == 5) if (nslow2 == 5)
info = 4; info = 4;
@ -375,7 +375,7 @@ L260:
} }
wa2[j] = (sum - wa3[j]) / pnorm; wa2[j] = (sum - wa3[j]) / pnorm;
wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm); wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
if (ratio >= p0001) { if (ratio >= Scalar(1e-4)) {
qtf[j] = sum; qtf[j] = sum;
} }
/* L280: */ /* L280: */
@ -389,7 +389,7 @@ L260:
/* end of the inner loop. */ /* end of the inner loop. */
jeval = FALSE_; jeval = false;
goto L180; goto L180;
L290: L290:

View File

@ -88,7 +88,7 @@ L20:
/* beginning of the outer loop. */ /* beginning of the outer loop. */
L30: L30:
jeval = TRUE_; jeval = true;
/* calculate the jacobian matrix. */ /* calculate the jacobian matrix. */
@ -100,7 +100,7 @@ L30:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
qrfac(n, n,fjac.data(), ldfjac, FALSE_, iwa, 1, wa1.data(), wa2.data(), wa3.data()); qrfac(n, n,fjac.data(), ldfjac, false, iwa, 1, wa1.data(), wa2.data(), wa3.data());
/* on the first iteration and if mode is 1, scale according */ /* on the first iteration and if mode is 1, scale according */
/* to the norms of the columns of the initial jacobian. */ /* to the norms of the columns of the initial jacobian. */
@ -161,7 +161,7 @@ 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;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
l = j; l = j;
if (j) { if (j) {
@ -173,7 +173,7 @@ L110:
} }
R[l] = wa1[j]; R[l] = wa1[j];
if (wa1[j] == 0.) { if (wa1[j] == 0.) {
sing = TRUE_; sing = true;
} }
/* L150: */ /* L150: */
} }
@ -189,7 +189,7 @@ L110:
} }
/* Computing MAX */ /* Computing MAX */
for (j = 0; j < n; ++j) for (j = 0; j < n; ++j)
diag[j] = max(diag[j], wa2[j]); diag[j] = std::max(diag[j], wa2[j]);
L170: L170:
/* beginning of the inner loop. */ /* beginning of the inner loop. */
@ -212,7 +212,7 @@ L190:
/* determine the direction p. */ /* determine the direction p. */
dogleg(n, R.data(), lr, diag.data(), qtf.data(), delta, wa1.data(), wa2.data(), wa3.data()); ei_dogleg<Scalar>(n, R.data(), lr, diag.data(), qtf.data(), delta, wa1.data(), wa2.data(), wa3.data());
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
@ -227,7 +227,7 @@ L190:
/* on the first iteration, adjust the initial step bound. */ /* on the first iteration, adjust the initial step bound. */
if (iter == 1) { if (iter == 1) {
delta = min(delta,pnorm); delta = std::min(delta,pnorm);
} }
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
@ -273,26 +273,26 @@ L190:
/* update the step bound. */ /* update the step bound. */
if (ratio >= p1) { if (ratio >= Scalar(.1)) {
goto L230; goto L230;
} }
ncsuc = 0; ncsuc = 0;
++ncfail; ++ncfail;
delta = p5 * delta; delta = Scalar(.5) * delta;
goto L240; goto L240;
L230: L230:
ncfail = 0; ncfail = 0;
++ncsuc; ++ncsuc;
if (ratio >= p5 || ncsuc > 1) /* Computing MAX */ if (ratio >= Scalar(.5) || ncsuc > 1) /* Computing MAX */
delta = max(delta, pnorm / p5); delta = std::max(delta, pnorm / Scalar(.5));
if (ei_abs(ratio - 1.) <= p1) { if (ei_abs(ratio - 1.) <= Scalar(.1)) {
delta = pnorm / p5; delta = pnorm / Scalar(.5);
} }
L240: L240:
/* test for successful iteration. */ /* test for successful iteration. */
if (ratio < p0001) { if (ratio < Scalar(1e-4)) {
goto L260; goto L260;
} }
@ -312,13 +312,13 @@ L260:
/* determine the progress of the iteration. */ /* determine the progress of the iteration. */
++nslow1; ++nslow1;
if (actred >= p001) { if (actred >= Scalar(.001)) {
nslow1 = 0; nslow1 = 0;
} }
if (jeval) { if (jeval) {
++nslow2; ++nslow2;
} }
if (actred >= p1) { if (actred >= Scalar(.1)) {
nslow2 = 0; nslow2 = 0;
} }
@ -337,7 +337,7 @@ L260:
info = 2; info = 2;
} }
/* Computing MAX */ /* Computing MAX */
if (p1 * max(p1 * delta, pnorm) <= epsilon<Scalar>() * xnorm) { if (Scalar(.1) * std::max(Scalar(.1) * delta, pnorm) <= epsilon<Scalar>() * xnorm) {
info = 3; info = 3;
} }
if (nslow2 == 5) { if (nslow2 == 5) {
@ -367,7 +367,7 @@ L260:
} }
wa2[j] = (sum - wa3[j]) / pnorm; wa2[j] = (sum - wa3[j]) / pnorm;
wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm); wa1[j] = diag[j] * (diag[j] * wa1[j] / pnorm);
if (ratio >= p0001) { if (ratio >= Scalar(1e-4)) {
qtf[j] = sum; qtf[j] = sum;
} }
/* L280: */ /* L280: */
@ -381,7 +381,7 @@ L260:
/* end of the inner loop. */ /* end of the inner loop. */
jeval = FALSE_; jeval = false;
goto L180; goto L180;
L290: L290:

View File

@ -104,7 +104,7 @@ L40:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
qrfac(m, n, fjac.data(), ldfjac, TRUE_, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data()); qrfac(m, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1) ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
/* on the first iteration and if mode is 1, scale according */ /* on the first iteration and if mode is 1, scale according */
@ -183,7 +183,7 @@ L120:
/* L140: */ /* L140: */
} }
/* Computing MAX */ /* Computing MAX */
gnorm = max(gnorm, ei_abs(sum / wa2[l])); gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
L150: L150:
/* L160: */ /* L160: */
; ;
@ -205,7 +205,7 @@ L170:
goto L190; goto L190;
} }
for (j = 0; j < n; ++j) for (j = 0; j < n; ++j)
diag[j] = max( diag[j], wa2[j]); diag[j] = std::max( diag[j], wa2[j]);
L190: L190:
/* beginning of the inner loop. */ /* beginning of the inner loop. */
@ -215,7 +215,7 @@ L200:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides) ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
&par, wa1.data(), wa2.data(), wa3.data(), wa4.data()); &par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
ipvt.cwise()-=1; ipvt.cwise()-=1;
@ -232,7 +232,7 @@ L200:
/* on the first iteration, adjust the initial step bound. */ /* on the first iteration, adjust the initial step bound. */
if (iter == 1) { if (iter == 1) {
delta = min(delta,pnorm); delta = std::min(delta,pnorm);
} }
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
@ -247,7 +247,7 @@ L200:
/* compute the scaled actual reduction. */ /* compute the scaled actual reduction. */
actred = -1.; actred = -1.;
if (p1 * fnorm1 < fnorm) /* Computing 2nd power */ if (Scalar(.1) * fnorm1 < fnorm) /* Computing 2nd power */
actred = 1. - ei_abs2(fnorm1 / fnorm); actred = 1. - ei_abs2(fnorm1 / fnorm);
/* compute the scaled predicted reduction and */ /* compute the scaled predicted reduction and */
@ -266,7 +266,7 @@ L200:
temp1 = ei_abs2(wa3.stableNorm() / fnorm); temp1 = ei_abs2(wa3.stableNorm() / fnorm);
temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm); temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm);
/* Computing 2nd power */ /* Computing 2nd power */
prered = temp1 + temp2 / p5; prered = temp1 + temp2 / Scalar(.5);
dirder = -(temp1 + temp2); dirder = -(temp1 + temp2);
/* compute the ratio of the actual to the predicted */ /* compute the ratio of the actual to the predicted */
@ -279,34 +279,34 @@ L200:
/* update the step bound. */ /* update the step bound. */
if (ratio > p25) { if (ratio > Scalar(.25)) {
goto L240; goto L240;
} }
if (actred >= 0.) { if (actred >= 0.) {
temp = p5; temp = Scalar(.5);
} }
if (actred < 0.) { if (actred < 0.) {
temp = p5 * dirder / (dirder + p5 * actred); temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
} }
if (p1 * fnorm1 >= fnorm || temp < p1) { if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
temp = p1; temp = Scalar(.1);
} }
/* Computing MIN */ /* Computing MIN */
delta = temp * min(delta, pnorm/p1); delta = temp * std::min(delta, pnorm/Scalar(.1));
par /= temp; par /= temp;
goto L260; goto L260;
L240: L240:
if (par != 0. && ratio < p75) { if (par != 0. && ratio < Scalar(.75)) {
goto L250; goto L250;
} }
delta = pnorm / p5; delta = pnorm / Scalar(.5);
par = p5 * par; par = Scalar(.5) * par;
L250: L250:
L260: L260:
/* test for successful iteration. */ /* test for successful iteration. */
if (ratio < p0001) { if (ratio < Scalar(1e-4)) {
goto L290; goto L290;
} }
@ -328,13 +328,13 @@ L290:
/* tests for convergence. */ /* tests for convergence. */
if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) { if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1.) {
info = 1; info = 1;
} }
if (delta <= xtol * xnorm) { if (delta <= xtol * xnorm) {
info = 2; info = 2;
} }
if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1. && info
== 2) { == 2) {
info = 3; info = 3;
} }
@ -347,7 +347,7 @@ L290:
if (nfev >= maxfev) { if (nfev >= maxfev) {
info = 5; info = 5;
} }
if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) { if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.) {
info = 6; info = 6;
} }
if (delta <= epsilon<Scalar>() * xnorm) { if (delta <= epsilon<Scalar>() * xnorm) {
@ -362,7 +362,7 @@ L290:
/* end of the inner loop. repeat if iteration unsuccessful. */ /* end of the inner loop. repeat if iteration unsuccessful. */
if (ratio < p0001) { if (ratio < Scalar(1e-4)) {
goto L200; goto L200;
} }

View File

@ -107,7 +107,7 @@ L40:
/* compute the qr factorization of the jacobian. */ /* compute the qr factorization of the jacobian. */
qrfac(m, n, fjac.data(), ldfjac, TRUE_, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data()); qrfac(m, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1) ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
/* on the first iteration and if mode is 1, scale according */ /* on the first iteration and if mode is 1, scale according */
@ -186,7 +186,7 @@ L120:
/* L140: */ /* L140: */
} }
/* Computing MAX */ /* Computing MAX */
gnorm = max(gnorm, ei_abs(sum / wa2[l])); gnorm = std::max(gnorm, ei_abs(sum / wa2[l]));
L150: L150:
/* L160: */ /* L160: */
; ;
@ -208,7 +208,7 @@ L170:
goto L190; goto L190;
} }
for (j = 0; j < n; ++j) /* Computing MAX */ for (j = 0; j < n; ++j) /* Computing MAX */
diag[j] = max(diag[j], wa2[j]); diag[j] = std::max(diag[j], wa2[j]);
L190: L190:
/* beginning of the inner loop. */ /* beginning of the inner loop. */
@ -218,7 +218,7 @@ L200:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides) ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta,
&par, wa1.data(), wa2.data(), wa3.data(), wa4.data()); &par, wa1.data(), wa2.data(), wa3.data(), wa4.data());
ipvt.cwise()-=1; ipvt.cwise()-=1;
@ -235,7 +235,7 @@ L200:
/* on the first iteration, adjust the initial step bound. */ /* on the first iteration, adjust the initial step bound. */
if (iter == 1) { if (iter == 1) {
delta = min(delta,pnorm); delta = std::min(delta,pnorm);
} }
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
@ -250,7 +250,7 @@ L200:
/* compute the scaled actual reduction. */ /* compute the scaled actual reduction. */
actred = -1.; actred = -1.;
if (p1 * fnorm1 < fnorm) /* Computing 2nd power */ if (Scalar(.1) * fnorm1 < fnorm) /* Computing 2nd power */
actred = 1. - ei_abs2(fnorm1 / fnorm); actred = 1. - ei_abs2(fnorm1 / fnorm);
/* compute the scaled predicted reduction and */ /* compute the scaled predicted reduction and */
@ -269,7 +269,7 @@ L200:
temp1 = ei_abs2(wa3.stableNorm() / fnorm); temp1 = ei_abs2(wa3.stableNorm() / fnorm);
temp2 = ei_abs2(ei_sqrt(par) * pnorm / fnorm); temp2 = ei_abs2(ei_sqrt(par) * pnorm / fnorm);
/* Computing 2nd power */ /* Computing 2nd power */
prered = temp1 + temp2 / p5; prered = temp1 + temp2 / Scalar(.5);
dirder = -(temp1 + temp2); dirder = -(temp1 + temp2);
/* compute the ratio of the actual to the predicted */ /* compute the ratio of the actual to the predicted */
@ -282,34 +282,34 @@ L200:
/* update the step bound. */ /* update the step bound. */
if (ratio > p25) { if (ratio > Scalar(.25)) {
goto L240; goto L240;
} }
if (actred >= 0.) { if (actred >= 0.) {
temp = p5; temp = Scalar(.5);
} }
if (actred < 0.) { if (actred < 0.) {
temp = p5 * dirder / (dirder + p5 * actred); temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
} }
if (p1 * fnorm1 >= fnorm || temp < p1) { if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
temp = p1; temp = Scalar(.1);
} }
/* Computing MIN */ /* Computing MIN */
delta = temp * min(delta, pnorm / p1); delta = temp * std::min(delta, pnorm / Scalar(.1));
par /= temp; par /= temp;
goto L260; goto L260;
L240: L240:
if (par != 0. && ratio < p75) { if (par != 0. && ratio < Scalar(.75)) {
goto L250; goto L250;
} }
delta = pnorm / p5; delta = pnorm / Scalar(.5);
par = p5 * par; par = Scalar(.5) * par;
L250: L250:
L260: L260:
/* test for successful iteration. */ /* test for successful iteration. */
if (ratio < p0001) { if (ratio < Scalar(1e-4)) {
goto L290; goto L290;
} }
@ -331,13 +331,13 @@ L290:
/* tests for convergence. */ /* tests for convergence. */
if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) { if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1.) {
info = 1; info = 1;
} }
if (delta <= xtol * xnorm) { if (delta <= xtol * xnorm) {
info = 2; info = 2;
} }
if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1. && info
== 2) { == 2) {
info = 3; info = 3;
} }
@ -350,7 +350,7 @@ L290:
if (nfev >= maxfev) { if (nfev >= maxfev) {
info = 5; info = 5;
} }
if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) { if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.) {
info = 6; info = 6;
} }
if (delta <= epsilon<Scalar>() * xnorm) { if (delta <= epsilon<Scalar>() * xnorm) {
@ -365,7 +365,7 @@ L290:
/* end of the inner loop. repeat if iteration unsuccessful. */ /* end of the inner loop. repeat if iteration unsuccessful. */
if (ratio < p0001) { if (ratio < Scalar(1e-4)) {
goto L200; goto L200;
} }

View File

@ -124,10 +124,10 @@ L40:
/* if the jacobian is rank deficient, call qrfac to */ /* if the jacobian is rank deficient, call qrfac to */
/* reorder its columns and update the components of qtf. */ /* reorder its columns and update the components of qtf. */
sing = FALSE_; sing = false;
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
if (fjac(j,j) == 0.) { if (fjac(j,j) == 0.) {
sing = TRUE_; sing = true;
} }
ipvt[j] = j; ipvt[j] = j;
wa2[j] = fjac.col(j).start(j).stableNorm(); wa2[j] = fjac.col(j).start(j).stableNorm();
@ -139,7 +139,7 @@ L40:
goto L130; goto L130;
} }
ipvt.cwise()+=1; ipvt.cwise()+=1;
qrfac(n, n, fjac.data(), ldfjac, TRUE_, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data()); qrfac(n, n, fjac.data(), ldfjac, true, ipvt.data(), n, wa1.data(), wa2.data(), wa3.data());
ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1) ipvt.cwise()-=1; // qrfac() creates ipvt with fortran convetion (1->n), convert it to c (0->n-1)
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
if (fjac(j,j) == 0.) { if (fjac(j,j) == 0.) {
@ -210,7 +210,7 @@ L170:
/* L180: */ /* L180: */
} }
/* Computing MAX */ /* Computing MAX */
gnorm = max(gnorm, ei_abs(sum/wa2[l])); gnorm = std::max(gnorm, ei_abs(sum/wa2[l]));
L190: L190:
/* L200: */ /* L200: */
; ;
@ -232,7 +232,7 @@ L210:
goto L230; goto L230;
} }
for (j = 0; j < n; ++j) /* Computing MAX */ for (j = 0; j < n; ++j) /* Computing MAX */
diag[j] = max(diag[j], wa2[j]); diag[j] = std::max(diag[j], wa2[j]);
L230: L230:
/* beginning of the inner loop. */ /* beginning of the inner loop. */
@ -242,7 +242,7 @@ L240:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides) ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, &par, ei_lmpar<Scalar>(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, &par,
wa1.data(), wa2.data(), wa3.data(), wa4.data()); wa1.data(), wa2.data(), wa3.data(), wa4.data());
ipvt.cwise()-=1; ipvt.cwise()-=1;
@ -259,7 +259,7 @@ L240:
/* on the first iteration, adjust the initial step bound. */ /* on the first iteration, adjust the initial step bound. */
if (iter == 1) { if (iter == 1) {
delta = min(delta,pnorm); delta = std::min(delta,pnorm);
} }
/* evaluate the function at x + p and calculate its norm. */ /* evaluate the function at x + p and calculate its norm. */
@ -274,7 +274,7 @@ L240:
/* compute the scaled actual reduction. */ /* compute the scaled actual reduction. */
actred = -1.; actred = -1.;
if (p1 * fnorm1 < fnorm) /* Computing 2nd power */ if (Scalar(.1) * fnorm1 < fnorm) /* Computing 2nd power */
actred = 1. - ei_abs2(fnorm1 / fnorm); actred = 1. - ei_abs2(fnorm1 / fnorm);
/* compute the scaled predicted reduction and */ /* compute the scaled predicted reduction and */
@ -293,7 +293,7 @@ L240:
temp1 = ei_abs2(wa3.stableNorm() / fnorm); temp1 = ei_abs2(wa3.stableNorm() / fnorm);
temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm); temp2 = ei_abs2( ei_sqrt(par) * pnorm / fnorm);
/* Computing 2nd power */ /* Computing 2nd power */
prered = temp1 + temp2 / p5; prered = temp1 + temp2 / Scalar(.5);
dirder = -(temp1 + temp2); dirder = -(temp1 + temp2);
/* compute the ratio of the actual to the predicted */ /* compute the ratio of the actual to the predicted */
@ -306,35 +306,35 @@ L240:
/* update the step bound. */ /* update the step bound. */
if (ratio > p25) { if (ratio > Scalar(.25)) {
goto L280; goto L280;
} }
if (actred >= 0.) { if (actred >= 0.) {
temp = p5; temp = Scalar(.5);
} }
if (actred < 0.) { if (actred < 0.) {
temp = p5 * dirder / (dirder + p5 * actred); temp = Scalar(.5) * dirder / (dirder + Scalar(.5) * actred);
} }
if (p1 * fnorm1 >= fnorm || temp < p1) { if (Scalar(.1) * fnorm1 >= fnorm || temp < Scalar(.1)) {
temp = p1; temp = Scalar(.1);
} }
/* Computing MIN */ /* Computing MIN */
delta = temp * min(delta, pnorm / p1); delta = temp * std::min(delta, pnorm / Scalar(.1));
par /= temp; par /= temp;
goto L300; goto L300;
L280: L280:
if (par != 0. && ratio < p75) { if (par != 0. && ratio < Scalar(.75)) {
goto L290; goto L290;
} }
delta = pnorm / p5; delta = pnorm / Scalar(.5);
par = p5 * par; par = Scalar(.5) * par;
L290: L290:
L300: L300:
/* test for successful iteration. */ /* test for successful iteration. */
if (ratio < p0001) { if (ratio < Scalar(1e-4)) {
goto L330; goto L330;
} }
@ -356,13 +356,13 @@ L330:
/* tests for convergence. */ /* tests for convergence. */
if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1.) { if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1.) {
info = 1; info = 1;
} }
if (delta <= xtol * xnorm) { if (delta <= xtol * xnorm) {
info = 2; info = 2;
} }
if (ei_abs(actred) <= ftol && prered <= ftol && p5 * ratio <= 1. && info if (ei_abs(actred) <= ftol && prered <= ftol && Scalar(.5) * ratio <= 1. && info
== 2) { == 2) {
info = 3; info = 3;
} }
@ -375,7 +375,7 @@ L330:
if (nfev >= maxfev) { if (nfev >= maxfev) {
info = 5; info = 5;
} }
if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && p5 * ratio <= 1.) { if (ei_abs(actred) <= epsilon<Scalar>() && prered <= epsilon<Scalar>() && Scalar(.5) * ratio <= 1.) {
info = 6; info = 6;
} }
if (delta <= epsilon<Scalar>() * xnorm) { if (delta <= epsilon<Scalar>() * xnorm) {
@ -390,7 +390,7 @@ L330:
/* end of the inner loop. repeat if iteration unsuccessful. */ /* end of the inner loop. repeat if iteration unsuccessful. */
if (ratio < p0001) { if (ratio < Scalar(1e-4)) {
goto L240; goto L240;
} }