bug #620: fix robustness issue in JacobiSVD::solve (also fix a perf. issue)

This commit is contained in:
Gael Guennebaud 2013-06-24 13:08:09 +02:00
parent d1d7a1ade9
commit 8bbde351e7

View File

@ -850,17 +850,13 @@ struct solve_retval<JacobiSVD<_MatrixType, QRPreconditioner>, Rhs>
// A = U S V^* // A = U S V^*
// So A^{-1} = V S^{-1} U^* // So A^{-1} = V S^{-1} U^*
Matrix<Scalar, Dynamic, Rhs::ColsAtCompileTime, 0, _MatrixType::MaxRowsAtCompileTime, Rhs::MaxColsAtCompileTime> tmp;
Index diagSize = (std::min)(dec().rows(), dec().cols()); Index diagSize = (std::min)(dec().rows(), dec().cols());
typename JacobiSVDType::SingularValuesType invertedSingVals(diagSize);
Index nonzeroSingVals = dec().nonzeroSingularValues(); Index nonzeroSingVals = dec().nonzeroSingularValues();
invertedSingVals.head(nonzeroSingVals) = dec().singularValues().head(nonzeroSingVals).array().inverse();
invertedSingVals.tail(diagSize - nonzeroSingVals).setZero(); tmp.noalias() = dec().matrixU().leftCols(nonzeroSingVals).adjoint() * rhs();
tmp = dec().singularValues().head(nonzeroSingVals).asDiagonal().inverse() * tmp;
dst = dec().matrixV().leftCols(diagSize) dst = dec().matrixV().leftCols(nonzeroSingVals) * tmp;
* invertedSingVals.asDiagonal()
* dec().matrixU().leftCols(diagSize).adjoint()
* rhs();
} }
}; };
} // end namespace internal } // end namespace internal