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 (-1 / +1 lines)
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
  }

Return to bug 1528