File indexing completed on 2025-01-18 09:57:05
0001 namespace Eigen {
0002
0003 namespace internal {
0004
0005 template <typename Scalar>
0006 void covar(
0007 Matrix< Scalar, Dynamic, Dynamic > &r,
0008 const VectorXi &ipvt,
0009 Scalar tol = std::sqrt(NumTraits<Scalar>::epsilon()) )
0010 {
0011 using std::abs;
0012 typedef DenseIndex Index;
0013
0014
0015 Index i, j, k, l, ii, jj;
0016 bool sing;
0017 Scalar temp;
0018
0019
0020 const Index n = r.cols();
0021 const Scalar tolr = tol * abs(r(0,0));
0022 Matrix< Scalar, Dynamic, 1 > wa(n);
0023 eigen_assert(ipvt.size()==n);
0024
0025
0026 l = -1;
0027 for (k = 0; k < n; ++k)
0028 if (abs(r(k,k)) > tolr) {
0029 r(k,k) = 1. / r(k,k);
0030 for (j = 0; j <= k-1; ++j) {
0031 temp = r(k,k) * r(j,k);
0032 r(j,k) = 0.;
0033 r.col(k).head(j+1) -= r.col(j).head(j+1) * temp;
0034 }
0035 l = k;
0036 }
0037
0038
0039
0040 for (k = 0; k <= l; ++k) {
0041 for (j = 0; j <= k-1; ++j)
0042 r.col(j).head(j+1) += r.col(k).head(j+1) * r(j,k);
0043 r.col(k).head(k+1) *= r(k,k);
0044 }
0045
0046
0047
0048 for (j = 0; j < n; ++j) {
0049 jj = ipvt[j];
0050 sing = j > l;
0051 for (i = 0; i <= j; ++i) {
0052 if (sing)
0053 r(i,j) = 0.;
0054 ii = ipvt[i];
0055 if (ii > jj)
0056 r(ii,jj) = r(i,j);
0057 if (ii < jj)
0058 r(jj,ii) = r(i,j);
0059 }
0060 wa[jj] = r(j,j);
0061 }
0062
0063
0064 r.topLeftCorner(n,n).template triangularView<StrictlyUpper>() = r.topLeftCorner(n,n).transpose();
0065 r.diagonal() = wa;
0066 }
0067
0068 }
0069
0070 }