This bugzilla service is closed. All entries have been migrated to https://gitlab.com/libeigen/eigen
View | Details | Raw Unified | Return to bug 1528 | Differences between
and this patch

Collapse All | Expand All

(-)a/Eigen/src/Cholesky/LDLT.h (-7 / +2 lines)
Lines 565-587 void LDLT<_MatrixType,_UpLo>::_solve_imp Link Here
565
565
566
  // dst = L^-1 (P b)
566
  // dst = L^-1 (P b)
567
  matrixL().solveInPlace(dst);
567
  matrixL().solveInPlace(dst);
568
568
569
  // dst = D^-1 (L^-1 P b)
569
  // dst = D^-1 (L^-1 P b)
570
  // more precisely, use pseudo-inverse of D (see bug 241)
570
  // more precisely, use pseudo-inverse of D (see bug 241)
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
  // The paranthesis are there to avoid min from being macro expanded to the cmath macro min(,)
574
  // as motivated by LAPACK's xGELSS:
574
  RealScalar tolerance = (std::numeric_limits<RealScalar>::min)(); 
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
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.
579
  RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest();
580
575
581
  for (Index i = 0; i < vecD.size(); ++i)
576
  for (Index i = 0; i < vecD.size(); ++i)
582
  {
577
  {
583
    if(abs(vecD(i)) > tolerance)
578
    if(abs(vecD(i)) > tolerance)
584
      dst.row(i) /= vecD(i);
579
      dst.row(i) /= vecD(i);
585
    else
580
    else
586
      dst.row(i).setZero();
581
      dst.row(i).setZero();
587
  }
582
  }

Return to bug 1528