LMcovar.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // This code initially comes from MINPACK whose original authors are:
00005 // Copyright Jorge More - Argonne National Laboratory
00006 // Copyright Burt Garbow - Argonne National Laboratory
00007 // Copyright Ken Hillstrom - Argonne National Laboratory
00008 //
00009 // This Source Code Form is subject to the terms of the Minpack license
00010 // (a BSD-like license) described in the campaigned CopyrightMINPACK.txt file.
00011 
00012 #ifndef EIGEN_LMCOVAR_H
00013 #define EIGEN_LMCOVAR_H
00014 
00015 namespace Eigen { 
00016 
00017 namespace internal {
00018 
00019 template <typename Scalar>
00020 void covar(
00021         Matrix< Scalar, Dynamic, Dynamic > &r,
00022         const VectorXi& ipvt,
00023         Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) )
00024 {
00025     using std::abs;
00026     /* Local variables */
00027     Index i, j, k, l, ii, jj;
00028     bool sing;
00029     Scalar temp;
00030 
00031     /* Function Body */
00032     const Index n = r.cols();
00033     const Scalar tolr = tol * abs(r(0,0));
00034     Matrix< Scalar, Dynamic, 1 > wa(n);
00035     eigen_assert(ipvt.size()==n);
00036 
00037     /* form the inverse of r in the full upper triangle of r. */
00038     l = -1;
00039     for (k = 0; k < n; ++k)
00040         if (abs(r(k,k)) > tolr) {
00041             r(k,k) = 1. / r(k,k);
00042             for (j = 0; j <= k-1; ++j) {
00043                 temp = r(k,k) * r(j,k);
00044                 r(j,k) = 0.;
00045                 r.col(k).head(j+1) -= r.col(j).head(j+1) * temp;
00046             }
00047             l = k;
00048         }
00049 
00050     /* form the full upper triangle of the inverse of (r transpose)*r */
00051     /* in the full upper triangle of r. */
00052     for (k = 0; k <= l; ++k) {
00053         for (j = 0; j <= k-1; ++j)
00054             r.col(j).head(j+1) += r.col(k).head(j+1) * r(j,k);
00055         r.col(k).head(k+1) *= r(k,k);
00056     }
00057 
00058     /* form the full lower triangle of the covariance matrix */
00059     /* in the strict lower triangle of r and in wa. */
00060     for (j = 0; j < n; ++j) {
00061         jj = ipvt[j];
00062         sing = j > l;
00063         for (i = 0; i <= j; ++i) {
00064             if (sing)
00065                 r(i,j) = 0.;
00066             ii = ipvt[i];
00067             if (ii > jj)
00068                 r(ii,jj) = r(i,j);
00069             if (ii < jj)
00070                 r(jj,ii) = r(i,j);
00071         }
00072         wa[jj] = r(j,j);
00073     }
00074 
00075     /* symmetrize the covariance matrix in r. */
00076     r.topLeftCorner(n,n).template triangularView<StrictlyUpper>() = r.topLeftCorner(n,n).transpose();
00077     r.diagonal() = wa;
00078 }
00079 
00080 } // end namespace internal
00081 
00082 } // end namespace Eigen
00083 
00084 #endif // EIGEN_LMCOVAR_H
 All Classes Functions Variables Typedefs Enumerator