GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-10-15 01:21:29 +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 
(grafted from fe630c98739186c1dd3faf7a06b52a64baa5ea23
)
This commit is contained in:
Gael Guennebaud 2015-10-18 22:15:01 +02:00
parent dc0ef2cbed
commit 9ff967199a
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

@ -289,7 +289,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;
}

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

@ -302,9 +302,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