![]() |
Eigen-unsupported
3.3.3
|
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