mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

eigenization of ei_covar()

Browse Source
This commit is contained in:
Thomas Capricelli 2009-08-24 17:47:35 +02:00
parent f69869c42a
commit 17905c7399
2 changed files with 45 additions and 60 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

93
unsupported/Eigen/src/NonLinear/covar.h
View File

@ -1,86 +1,75 @@
template <typename Scalar>
void ei_covar(int n, Scalar *r__, int ldr,
const int *ipvt, Scalar tol, Scalar *wa)
void ei_covar(
Matrix< Scalar, Dynamic, Dynamic > &r,
const VectorXi &ipvt,
Scalar tol = ei_sqrt(epsilon<Scalar>()) )
{
/* System generated locals */
int r_dim1, r_offset;
/* Local variables */
int i, j, k, l, ii, jj, km1;
int i, j, k, l, ii, jj;
int sing;
Scalar temp, tolr;
/* Parameter adjustments */
--wa;
--ipvt;
tolr = tol * ei_abs(r__[0]);
r_dim1 = ldr;
r_offset = 1 + r_dim1;
r__ -= r_offset;
Scalar temp;
/* Function Body */
const int n = r.cols();
const Scalar tolr = tol * ei_abs(r[0]);
Matrix< Scalar, Dynamic, 1 > wa(n);
assert(ipvt.size()==n);
/* form the inverse of r in the full upper triangle of r. */
l = 0;
for (k = 1; k <= n; ++k) {
if (ei_abs(r__[k + k * r_dim1]) > tolr) {
r__[k + k * r_dim1] = 1. / r__[k + k * r_dim1];
km1 = k - 1;
if (km1 >= 1)
for (j = 1; j <= km1; ++j) {
temp = r__[k + k * r_dim1] * r__[j + k * r_dim1];
r__[j + k * r_dim1] = 0.;
for (i = 1; i <= j; ++i) {
r__[i + k * r_dim1] -= temp * r__[i + j * r_dim1];
}
l = -1;
for (k = 0; k < n; ++k)
if (ei_abs(r(k,k)) > tolr) {
r(k,k) = 1. / r(k,k);
for (j = 0; j <= k-1; ++j) {
temp = r(k,k) * r(j,k);
r(j,k) = 0.;
for (i = 0; i <= j; ++i) {
r(i,k) -= temp * r(i,j);
}
}
l = k;
}
}
/* form the full upper triangle of the inverse of (r transpose)*r */
/* in the full upper triangle of r. */
if (l >= 1)
for (k = 1; k <= l; ++k) {
km1 = k - 1;
if (km1 >= 1)
for (j = 1; j <= km1; ++j) {
temp = r__[j + k * r_dim1];
for (i = 1; i <= j; ++i)
r__[i + j * r_dim1] += temp * r__[i + k * r_dim1];
}
temp = r__[k + k * r_dim1];
for (i = 1; i <= k; ++i)
r__[i + k * r_dim1] = temp * r__[i + k * r_dim1];
for (k = 0; k <= l; ++k) {
for (j = 0; j <= k-1; ++j) {
temp = r(j,k);
for (i = 0; i <= j; ++i)
r(i,j) += temp * r(i,k);
}
temp = r(k,k);
for (i = 0; i <= k; ++i)
r(i,k) = temp * r(i,k);
}
/* form the full lower triangle of the covariance matrix */
/* in the strict lower triangle of r and in wa. */
for (j = 1; j <= n; ++j) {
jj = ipvt[j];
for (j = 0; j < n; ++j) {
jj = ipvt[j]-1;
sing = j > l;
for (i = 1; i <= j; ++i) {
for (i = 0; i <= j; ++i) {
if (sing)
r__[i + j * r_dim1] = 0.;
ii = ipvt[i];
r(i,j) = 0.;
ii = ipvt[i]-1;
if (ii > jj)
r__[ii + jj * r_dim1] = r__[i + j * r_dim1];
r(ii,jj) = r(i,j);
if (ii < jj)
r__[jj + ii * r_dim1] = r__[i + j * r_dim1];
r(jj,ii) = r(i,j);
}
wa[jj] = r__[j + j * r_dim1];
wa[jj] = r(j,j);
}
/* symmetrize the covariance matrix in r. */
for (j = 1; j <= n; ++j) {
for (i = 1; i <= j; ++i)
r__[i + j * r_dim1] = r__[j + i * r_dim1];
r__[j + j * r_dim1] = wa[j];
for (j = 0; j < n; ++j) {
for (i = 0; i <= j; ++i)
r(i,j) = r(j,i);
r(j,j) = wa[j];
}
}

12
unsupported/test/NonLinear.cpp
View File

@ -166,7 +166,7 @@ void testLmder()
{
const int m=15, n=3;
int info, nfev=0, njev=0;
double fnorm, covfac, covar_ftol;
double fnorm, covfac;
VectorXd x(n), fvec(m), diag(n), qtf;
MatrixXd fjac;
VectorXi ipvt;
@ -192,11 +192,9 @@ void testLmder()
VERIFY_IS_APPROX(x, x_ref);
// check covariance
covar_ftol = epsilon<double>();
covfac = fnorm*fnorm/(m-n);
VectorXd wa(n);
ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
ei_covar<double>(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
ei_covar<double>(fjac, ipvt);
MatrixXd cov_ref(n,n);
cov_ref <<
@ -543,7 +541,7 @@ void testLmdif()
{
const int m=15, n=3;
int info, nfev=0;
double fnorm, covfac, covar_ftol;
double fnorm, covfac;
VectorXd x(n), fvec(m), diag(n), qtf;
MatrixXd fjac;
VectorXi ipvt;
@ -568,11 +566,9 @@ void testLmdif()
VERIFY_IS_APPROX(x, x_ref);
// check covariance
covar_ftol = epsilon<double>();
covfac = fnorm*fnorm/(m-n);
VectorXd wa(n);
ipvt.cwise()+=1; // covar() expects the fortran convention (as qrfac provides)
ei_covar<double>(n, fjac.data(), m, ipvt.data(), covar_ftol, wa.data());
ei_covar<double>(fjac, ipvt);
MatrixXd cov_ref(n,n);
cov_ref <<