mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2024-12-15 07:10:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

be (hopefully) smarter with indices convention : we keep the c convention

Browse Source 
(0->n-1) as much as possible, and only convert at borders with
fortran-expecting methods, that will eventually dissapear.
This commit is contained in:
Thomas Capricelli 2009-08-21 04:16:37 +02:00
parent 524e112ee5
commit 2e3d17c3ce
3 changed files with 14 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
unsupported/Eigen/src/NonLinear/lmder.h
View File

@ -105,6 +105,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)
/* 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. */
@ -172,7 +173,7 @@ L120:
goto L170; goto L170;
} }
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
l = ipvt[j]-1; l = ipvt[j];
if (wa2[l] == 0.) { if (wa2[l] == 0.) {
goto L150; goto L150;
} }
@ -213,8 +214,10 @@ L200:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, lmpar(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;
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
@ -252,7 +255,7 @@ L200:
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
wa3[j] = 0.; wa3[j] = 0.;
l = ipvt[j]-1; l = ipvt[j];
temp = wa1[l]; temp = wa1[l];
for (i = 0; i <= j; ++i) { for (i = 0; i <= j; ++i) {
wa3[i] += fjac(i,j) * temp; wa3[i] += fjac(i,j) * temp;

10
unsupported/Eigen/src/NonLinear/lmstr.h
View File

@ -129,7 +129,7 @@ L40:
if (fjac(j,j) == 0.) { if (fjac(j,j) == 0.) {
sing = TRUE_; sing = TRUE_;
} }
ipvt[j] = j+1; ipvt[j] = j;
wa2[j] = fjac.col(j).start(j).stableNorm(); wa2[j] = fjac.col(j).start(j).stableNorm();
// wa2[j] = ei_enorm<Scalar>(j, &fjac[j * ldfjac + 1]); // wa2[j] = ei_enorm<Scalar>(j, &fjac[j * ldfjac + 1]);
// sum += fjac[i + j * ldfjac] * (qtf[i] / fnorm); // sum += fjac[i + j * ldfjac] * (qtf[i] / fnorm);
@ -138,7 +138,9 @@ L40:
if (! sing) { if (! sing) {
goto L130; goto L130;
} }
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)
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
if (fjac(j,j) == 0.) { if (fjac(j,j) == 0.) {
goto L110; goto L110;
@ -198,7 +200,7 @@ L170:
goto L210; goto L210;
} }
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
l = ipvt[j]-1; l = ipvt[j];
if (wa2[l] == 0.) { if (wa2[l] == 0.) {
goto L190; goto L190;
} }
@ -239,8 +241,10 @@ L240:
/* determine the levenberg-marquardt parameter. */ /* determine the levenberg-marquardt parameter. */
ipvt.cwise()+=1; // lmpar() expects the fortran convention (as qrfac provides)
lmpar(n, fjac.data(), ldfjac, ipvt.data(), diag.data(), qtf.data(), delta, &par, lmpar(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;
/* store the direction p and x + p. calculate the norm of p. */ /* store the direction p and x + p. calculate the norm of p. */
@ -278,7 +282,7 @@ L240:
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
wa3[j] = 0.; wa3[j] = 0.;
l = ipvt[j]-1; l = ipvt[j];
temp = wa1[l]; temp = wa1[l];
for (i = 0; i <= j; ++i) { for (i = 0; i <= j; ++i) {
wa3[i] += fjac(i,j) * temp; wa3[i] += fjac(i,j) * temp;

2
unsupported/test/NonLinear.cpp
View File

@ -246,6 +246,7 @@ void testLmder()
covar_ftol = epsilon<double>(); covar_ftol = epsilon<double>();
covfac = fnorm*fnorm/(m-n); covfac = fnorm*fnorm/(m-n);
VectorXd wa(n); VectorXd wa(n);
ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data()); covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
MatrixXd cov_ref(n,n); MatrixXd cov_ref(n,n);
@ -761,6 +762,7 @@ void testLmdif()
covar_ftol = epsilon<double>(); covar_ftol = epsilon<double>();
covfac = fnorm*fnorm/(m-n); covfac = fnorm*fnorm/(m-n);
VectorXd wa(n); VectorXd wa(n);
// ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data()); covar(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
MatrixXd cov_ref(n,n); MatrixXd cov_ref(n,n);