1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // This code initially comes from MINPACK whose original authors are: 5 // Copyright Jorge More - Argonne National Laboratory 6 // Copyright Burt Garbow - Argonne National Laboratory 7 // Copyright Ken Hillstrom - Argonne National Laboratory 8 // 9 // This Source Code Form is subject to the terms of the Minpack license 10 // (a BSD-like license) described in the campaigned CopyrightMINPACK.txt file. 11 12 #ifndef EIGEN_LMCOVAR_H 13 #define EIGEN_LMCOVAR_H 14 15 namespace Eigen { 16 17 namespace internal { 18 19 template <typename Scalar> 20 void covar( 21 Matrix< Scalar, Dynamic, Dynamic > &r, 22 const VectorXi& ipvt, 23 Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) ) 24 { 25 using std::abs; 26 typedef DenseIndex Index; 27 /* Local variables */ 28 Index i, j, k, l, ii, jj; 29 bool sing; 30 Scalar temp; 31 32 /* Function Body */ 33 const Index n = r.cols(); 34 const Scalar tolr = tol * abs(r(0,0)); 35 Matrix< Scalar, Dynamic, 1 > wa(n); 36 eigen_assert(ipvt.size()==n); 37 38 /* form the inverse of r in the full upper triangle of r. */ 39 l = -1; 40 for (k = 0; k < n; ++k) 41 if (abs(r(k,k)) > tolr) { 42 r(k,k) = 1. / r(k,k); 43 for (j = 0; j <= k-1; ++j) { 44 temp = r(k,k) * r(j,k); 45 r(j,k) = 0.; 46 r.col(k).head(j+1) -= r.col(j).head(j+1) * temp; 47 } 48 l = k; 49 } 50 51 /* form the full upper triangle of the inverse of (r transpose)*r */ 52 /* in the full upper triangle of r. */ 53 for (k = 0; k <= l; ++k) { 54 for (j = 0; j <= k-1; ++j) 55 r.col(j).head(j+1) += r.col(k).head(j+1) * r(j,k); 56 r.col(k).head(k+1) *= r(k,k); 57 } 58 59 /* form the full lower triangle of the covariance matrix */ 60 /* in the strict lower triangle of r and in wa. */ 61 for (j = 0; j < n; ++j) { 62 jj = ipvt[j]; 63 sing = j > l; 64 for (i = 0; i <= j; ++i) { 65 if (sing) 66 r(i,j) = 0.; 67 ii = ipvt[i]; 68 if (ii > jj) 69 r(ii,jj) = r(i,j); 70 if (ii < jj) 71 r(jj,ii) = r(i,j); 72 } 73 wa[jj] = r(j,j); 74 } 75 76 /* symmetrize the covariance matrix in r. */ 77 r.topLeftCorner(n,n).template triangularView<StrictlyUpper>() = r.topLeftCorner(n,n).transpose(); 78 r.diagonal() = wa; 79 } 80 81 } // end namespace internal 82 83 } // end namespace Eigen 84 85 #endif // EIGEN_LMCOVAR_H 86
