$treeview $search $mathjax
Eigen-unsupported
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 typedef DenseIndex Index; 00027 /* Local variables */ 00028 Index i, j, k, l, ii, jj; 00029 bool sing; 00030 Scalar temp; 00031 00032 /* Function Body */ 00033 const Index n = r.cols(); 00034 const Scalar tolr = tol * abs(r(0,0)); 00035 Matrix< Scalar, Dynamic, 1 > wa(n); 00036 eigen_assert(ipvt.size()==n); 00037 00038 /* form the inverse of r in the full upper triangle of r. */ 00039 l = -1; 00040 for (k = 0; k < n; ++k) 00041 if (abs(r(k,k)) > tolr) { 00042 r(k,k) = 1. / r(k,k); 00043 for (j = 0; j <= k-1; ++j) { 00044 temp = r(k,k) * r(j,k); 00045 r(j,k) = 0.; 00046 r.col(k).head(j+1) -= r.col(j).head(j+1) * temp; 00047 } 00048 l = k; 00049 } 00050 00051 /* form the full upper triangle of the inverse of (r transpose)*r */ 00052 /* in the full upper triangle of r. */ 00053 for (k = 0; k <= l; ++k) { 00054 for (j = 0; j <= k-1; ++j) 00055 r.col(j).head(j+1) += r.col(k).head(j+1) * r(j,k); 00056 r.col(k).head(k+1) *= r(k,k); 00057 } 00058 00059 /* form the full lower triangle of the covariance matrix */ 00060 /* in the strict lower triangle of r and in wa. */ 00061 for (j = 0; j < n; ++j) { 00062 jj = ipvt[j]; 00063 sing = j > l; 00064 for (i = 0; i <= j; ++i) { 00065 if (sing) 00066 r(i,j) = 0.; 00067 ii = ipvt[i]; 00068 if (ii > jj) 00069 r(ii,jj) = r(i,j); 00070 if (ii < jj) 00071 r(jj,ii) = r(i,j); 00072 } 00073 wa[jj] = r(j,j); 00074 } 00075 00076 /* symmetrize the covariance matrix in r. */ 00077 r.topLeftCorner(n,n).template triangularView<StrictlyUpper>() = r.topLeftCorner(n,n).transpose(); 00078 r.diagonal() = wa; 00079 } 00080 00081 } // end namespace internal 00082 00083 } // end namespace Eigen 00084 00085 #endif // EIGEN_LMCOVAR_H