Back to home page

EIC code displayed by LXR

 
 

    


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     /* Local variables */
0015     Index i, j, k, l, ii, jj;
0016     bool sing;
0017     Scalar temp;
0018 
0019     /* Function Body */
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     /* form the inverse of r in the full upper triangle of r. */
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     /* form the full upper triangle of the inverse of (r transpose)*r */
0039     /* in the full upper triangle of r. */
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     /* form the full lower triangle of the covariance matrix */
0047     /* in the strict lower triangle of r and in wa. */
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     /* symmetrize the covariance matrix in r. */
0064     r.topLeftCorner(n,n).template triangularView<StrictlyUpper>() = r.topLeftCorner(n,n).transpose();
0065     r.diagonal() = wa;
0066 }
0067 
0068 } // end namespace internal
0069 
0070 } // end namespace Eigen