|
Lines 571-587
void LDLT<_MatrixType,_UpLo>::_solve_imp
Link Here
|
| 571 |
using std::abs; |
571 |
using std::abs; |
| 572 |
const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD()); |
572 |
const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD()); |
| 573 |
// In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon |
573 |
// In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon |
| 574 |
// as motivated by LAPACK's xGELSS: |
574 |
// as motivated by LAPACK's xGELSS: |
| 575 |
// RealScalar tolerance = numext::maxi(vecD.array().abs().maxCoeff() * NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest()); |
575 |
// RealScalar tolerance = numext::maxi(vecD.array().abs().maxCoeff() * NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest()); |
| 576 |
// However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest |
576 |
// However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest |
| 577 |
// diagonal element is not well justified and leads to numerical issues in some cases. |
577 |
// diagonal element is not well justified and leads to numerical issues in some cases. |
| 578 |
// Moreover, Lapack's xSYTRS routines use 0 for the tolerance. |
578 |
// Moreover, Lapack's xSYTRS routines use 0 for the tolerance. |
| 579 |
RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest(); |
579 |
RealScalar tolerance = (std::numeric_limits<RealScalar>::min)(); |
| 580 |
|
580 |
|
| 581 |
for (Index i = 0; i < vecD.size(); ++i) |
581 |
for (Index i = 0; i < vecD.size(); ++i) |
| 582 |
{ |
582 |
{ |
| 583 |
if(abs(vecD(i)) > tolerance) |
583 |
if(abs(vecD(i)) > tolerance) |
| 584 |
dst.row(i) /= vecD(i); |
584 |
dst.row(i) /= vecD(i); |
| 585 |
else |
585 |
else |
| 586 |
dst.row(i).setZero(); |
586 |
dst.row(i).setZero(); |
| 587 |
} |
587 |
} |