aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen
diff options
context:
space:
mode:
Diffstat (limited to 'Eigen')
-rw-r--r--Eigen/src/SVD/BDCSVD.h32
-rw-r--r--Eigen/src/SVD/JacobiSVD.h9
-rw-r--r--Eigen/src/SVD/SVDBase.h41
3 files changed, 64 insertions, 18 deletions
diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h
index e0c4456c7..17f8e4436 100644
--- a/Eigen/src/SVD/BDCSVD.h
+++ b/Eigen/src/SVD/BDCSVD.h
@@ -208,6 +208,7 @@ protected:
using Base::m_computeThinV;
using Base::m_matrixU;
using Base::m_matrixV;
+ using Base::m_info;
using Base::m_isInitialized;
using Base::m_nonzeroSingularValues;
@@ -256,16 +257,25 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
{
// FIXME this line involves temporaries
JacobiSVD<MatrixType> jsvd(matrix,computationOptions);
- if(computeU()) m_matrixU = jsvd.matrixU();
- if(computeV()) m_matrixV = jsvd.matrixV();
- m_singularValues = jsvd.singularValues();
- m_nonzeroSingularValues = jsvd.nonzeroSingularValues();
m_isInitialized = true;
+ m_info = jsvd.info();
+ if (m_info == Success || m_info == NoConvergence) {
+ if(computeU()) m_matrixU = jsvd.matrixU();
+ if(computeV()) m_matrixV = jsvd.matrixV();
+ m_singularValues = jsvd.singularValues();
+ m_nonzeroSingularValues = jsvd.nonzeroSingularValues();
+ }
return *this;
}
//**** step 0 - Copy the input matrix and apply scaling to reduce over/under-flows
- RealScalar scale = matrix.cwiseAbs().maxCoeff();
+ RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
+ if (!(numext::isfinite)(scale)) {
+ m_isInitialized = true;
+ m_info = InvalidInput;
+ return *this;
+ }
+
if(scale==Literal(0)) scale = Literal(1);
MatrixX copy;
if (m_isTranspose) copy = matrix.adjoint()/scale;
@@ -282,7 +292,11 @@ BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
m_computed.template bottomRows<1>().setZero();
divide(0, m_diagSize - 1, 0, 0, 0);
-
+ if (m_info != Success && m_info != NoConvergence) {
+ m_isInitialized = true;
+ return *this;
+ }
+
//**** step 3 - Copy singular values and vectors
for (int i=0; i<m_diagSize; i++)
{
@@ -394,7 +408,7 @@ void BDCSVD<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> A, co
//@param shift : Each time one takes the left submatrix, one must add 1 to the shift. Why? Because! We actually want the last column of the U submatrix
// to become the first column (*coeff) and to shift all the other columns to the right. There are more details on the reference paper.
template<typename MatrixType>
-void BDCSVD<MatrixType>::divide (Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift)
+void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift)
{
// requires rows = cols + 1;
using std::pow;
@@ -414,6 +428,8 @@ void BDCSVD<MatrixType>::divide (Eigen::Index firstCol, Eigen::Index lastCol, Ei
{
// FIXME this line involves temporaries
JacobiSVD<MatrixXr> b(m_computed.block(firstCol, firstCol, n + 1, n), ComputeFullU | (m_compV ? ComputeFullV : 0));
+ m_info = b.info();
+ if (m_info != Success && m_info != NoConvergence) return;
if (m_compU)
m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = b.matrixU();
else
@@ -433,7 +449,9 @@ void BDCSVD<MatrixType>::divide (Eigen::Index firstCol, Eigen::Index lastCol, Ei
// and the divide of the right submatrice reads one column of the left submatrice. That's why we need to treat the
// right submatrix before the left one.
divide(k + 1 + firstCol, lastCol, k + 1 + firstRowW, k + 1 + firstColW, shift);
+ if (m_info != Success && m_info != NoConvergence) return;
divide(firstCol, k - 1 + firstCol, firstRowW, firstColW + 1, shift + 1);
+ if (m_info != Success && m_info != NoConvergence) return;
if (m_compU)
{
diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h
index 2b6891105..a22a2e5c3 100644
--- a/Eigen/src/SVD/JacobiSVD.h
+++ b/Eigen/src/SVD/JacobiSVD.h
@@ -585,6 +585,7 @@ template<typename _MatrixType, int QRPreconditioner> class JacobiSVD
using Base::m_matrixU;
using Base::m_matrixV;
using Base::m_singularValues;
+ using Base::m_info;
using Base::m_isInitialized;
using Base::m_isAllocated;
using Base::m_usePrescribedThreshold;
@@ -625,6 +626,7 @@ void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Eigen::Index rows, Eigen:
m_rows = rows;
m_cols = cols;
+ m_info = Success;
m_isInitialized = false;
m_isAllocated = true;
m_computationOptions = computationOptions;
@@ -674,7 +676,12 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
// Scaling factor to reduce over/under-flows
- RealScalar scale = matrix.cwiseAbs().maxCoeff();
+ RealScalar scale = matrix.cwiseAbs().template maxCoeff<PropagateNaN>();
+ if (!(numext::isfinite)(scale)) {
+ m_isInitialized = true;
+ m_info = InvalidInput;
+ return *this;
+ }
if(scale==RealScalar(0)) scale = RealScalar(1);
/*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
diff --git a/Eigen/src/SVD/SVDBase.h b/Eigen/src/SVD/SVDBase.h
index 34d5c9dd3..13c8394ae 100644
--- a/Eigen/src/SVD/SVDBase.h
+++ b/Eigen/src/SVD/SVDBase.h
@@ -52,7 +52,7 @@ template<typename Derived> struct traits<SVDBase<Derived> >
* singular vectors. Asking for \em thin \a U or \a V means asking for only their \a m first columns to be formed. So \a U is then a n-by-m matrix,
* and \a V is then a p-by-m matrix. Notice that thin \a U and \a V are all you need for (least squares) solving.
*
- * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to
+ * If the input matrix has inf or nan coefficients, the result of the computation is undefined, and \a info() will return \a InvalidInput, but the computation is guaranteed to
* terminate in finite (and reasonable) time.
* \sa class BDCSVD, class JacobiSVD
*/
@@ -97,7 +97,7 @@ public:
*/
const MatrixUType& matrixU() const
{
- eigen_assert(m_isInitialized && "SVD is not initialized.");
+ _check_compute_assertions();
eigen_assert(computeU() && "This SVD decomposition didn't compute U. Did you ask for it?");
return m_matrixU;
}
@@ -113,7 +113,7 @@ public:
*/
const MatrixVType& matrixV() const
{
- eigen_assert(m_isInitialized && "SVD is not initialized.");
+ _check_compute_assertions();
eigen_assert(computeV() && "This SVD decomposition didn't compute V. Did you ask for it?");
return m_matrixV;
}
@@ -125,14 +125,14 @@ public:
*/
const SingularValuesType& singularValues() const
{
- eigen_assert(m_isInitialized && "SVD is not initialized.");
+ _check_compute_assertions();
return m_singularValues;
}
/** \returns the number of singular values that are not exactly 0 */
Index nonzeroSingularValues() const
{
- eigen_assert(m_isInitialized && "SVD is not initialized.");
+ _check_compute_assertions();
return m_nonzeroSingularValues;
}
@@ -145,7 +145,7 @@ public:
inline Index rank() const
{
using std::abs;
- eigen_assert(m_isInitialized && "JacobiSVD is not initialized.");
+ _check_compute_assertions();
if(m_singularValues.size()==0) return 0;
RealScalar premultiplied_threshold = numext::maxi<RealScalar>(m_singularValues.coeff(0) * threshold(), (std::numeric_limits<RealScalar>::min)());
Index i = m_nonzeroSingularValues-1;
@@ -224,6 +224,18 @@ public:
solve(const MatrixBase<Rhs>& b) const;
#endif
+
+ /** \brief Reports whether previous computation was successful.
+ *
+ * \returns \c Success if computation was successful.
+ */
+ EIGEN_DEVICE_FUNC
+ ComputationInfo info() const
+ {
+ eigen_assert(m_isInitialized && "SVD is not initialized.");
+ return m_info;
+ }
+
#ifndef EIGEN_PARSED_BY_DOXYGEN
template<typename RhsType, typename DstType>
void _solve_impl(const RhsType &rhs, DstType &dst) const;
@@ -233,26 +245,33 @@ public:
#endif
protected:
-
+
static void check_template_parameters()
{
EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
}
+ void _check_compute_assertions() const {
+ eigen_assert(m_isInitialized && "SVD is not initialized.");
+ eigen_assert(m_info != InvalidInput && "SVD failed due to invalid input.");
+ eigen_assert(m_info != NumericalIssue && "SVD failed due to invalid input.");
+ }
+
template<bool Transpose_, typename Rhs>
void _check_solve_assertion(const Rhs& b) const {
EIGEN_ONLY_USED_FOR_DEBUG(b);
- eigen_assert(m_isInitialized && "SVD is not initialized.");
+ _check_compute_assertions();
eigen_assert(computeU() && computeV() && "SVDBase::solve(): Both unitaries U and V are required to be computed (thin unitaries suffice).");
eigen_assert((Transpose_?cols():rows())==b.rows() && "SVDBase::solve(): invalid number of rows of the right hand side matrix b");
}
-
+
// return true if already allocated
bool allocate(Index rows, Index cols, unsigned int computationOptions) ;
MatrixUType m_matrixU;
MatrixVType m_matrixV;
SingularValuesType m_singularValues;
+ ComputationInfo m_info;
bool m_isInitialized, m_isAllocated, m_usePrescribedThreshold;
bool m_computeFullU, m_computeThinU;
bool m_computeFullV, m_computeThinV;
@@ -265,7 +284,8 @@ protected:
* Default constructor of SVDBase
*/
SVDBase()
- : m_isInitialized(false),
+ : m_info(Success),
+ m_isInitialized(false),
m_isAllocated(false),
m_usePrescribedThreshold(false),
m_computeFullU(false),
@@ -327,6 +347,7 @@ bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computat
m_rows = rows;
m_cols = cols;
+ m_info = Success;
m_isInitialized = false;
m_isAllocated = true;
m_computationOptions = computationOptions;
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-17 03:02:34 +0000