mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-24 14:45:14 +08:00
Code Issues Packages Projects Releases Wiki Activity

Improve numerical accuracy in LLT and triangular solve by using true scalar divisions (instead of x * (1/y))

Browse Source
This commit is contained in:
Gael Guennebaud 2015-10-18 22:15:01 +02:00
parent e99279f444
commit fe630c9873
2 changed files with 7 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Eigen/src/Cholesky/LLT.h
View File

@ -285,7 +285,7 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower>
return k; return k;
mat.coeffRef(k,k) = x = sqrt(x); mat.coeffRef(k,k) = x = sqrt(x);
if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint();
if (rs>0) A21 *= RealScalar(1)/x; if (rs>0) A21 /= x;
} }
return -1; return -1;
} }

9
Eigen/src/Core/products/TriangularSolverMatrix.h
View File

@ -304,9 +304,12 @@ EIGEN_DONT_INLINE void triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conj
for (Index i=0; i<actual_mc; ++i) for (Index i=0; i<actual_mc; ++i)
r[i] -= a[i] * b; r[i] -= a[i] * b;
} }
Scalar b = (Mode & UnitDiag) ? Scalar(1) : Scalar(1)/conj(rhs(j,j)); if((Mode & UnitDiag)==0)
for (Index i=0; i<actual_mc; ++i) {
r[i] *= b; Scalar b = conj(rhs(j,j));
for (Index i=0; i<actual_mc; ++i)
r[i] /= b;
}
} }
// pack the just computed part of lhs to A // pack the just computed part of lhs to A