diff options
author | Gael Guennebaud <g.gael@free.fr> | 2015-10-18 22:15:01 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2015-10-18 22:15:01 +0200 |
commit | fe630c98739186c1dd3faf7a06b52a64baa5ea23 (patch) | |
tree | dee931dcbc4ddbc803b5ed25e7316cbd01267a2f /Eigen | |
parent | e99279f444232c81901a84b4d15dae5fcf6d7d37 (diff) |
Improve numerical accuracy in LLT and triangular solve by using true scalar divisions (instead of x * (1/y))
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/src/Cholesky/LLT.h | 2 | ||||
-rw-r--r-- | Eigen/src/Core/products/TriangularSolverMatrix.h | 9 |
2 files changed, 7 insertions, 4 deletions
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index dc73304e8..1f0091f3c 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -285,7 +285,7 @@ template<typename Scalar> struct llt_inplace<Scalar, Lower> return k; mat.coeffRef(k,k) = x = sqrt(x); if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); - if (rs>0) A21 *= RealScalar(1)/x; + if (rs>0) A21 /= x; } return -1; } diff --git a/Eigen/src/Core/products/TriangularSolverMatrix.h b/Eigen/src/Core/products/TriangularSolverMatrix.h index a9a198d64..208593718 100644 --- a/Eigen/src/Core/products/TriangularSolverMatrix.h +++ b/Eigen/src/Core/products/TriangularSolverMatrix.h @@ -304,9 +304,12 @@ EIGEN_DONT_INLINE void triangular_solve_matrix<Scalar,Index,OnTheRight,Mode,Conj for (Index i=0; i<actual_mc; ++i) r[i] -= a[i] * b; } - Scalar b = (Mode & UnitDiag) ? Scalar(1) : Scalar(1)/conj(rhs(j,j)); - for (Index i=0; i<actual_mc; ++i) - r[i] *= b; + if((Mode & UnitDiag)==0) + { + 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 |