Fix -Wmaybe-uninitialized in SVD

This commit is contained in:
Charles Schlosser 2023-07-25 22:22:17 +00:00
parent 4e9e493b4a
commit 9995c3da6f
8 changed files with 200 additions and 154 deletions

View File

@ -107,6 +107,7 @@ class BDCSVD : public SVDBase<BDCSVD<MatrixType_, Options_> > {
public:
using Base::rows;
using Base::cols;
using Base::diagSize;
using Base::computeU;
using Base::computeV;
@ -270,7 +271,6 @@ private:
using Base::m_computationOptions;
using Base::m_computeThinU;
using Base::m_computeThinV;
using Base::m_diagSize;
using Base::m_info;
using Base::m_isInitialized;
using Base::m_matrixU;
@ -291,7 +291,7 @@ void BDCSVD<MatrixType, Options>::allocate(Index rows, Index cols, unsigned int
if (cols < m_algoswap)
internal::allocate_small_svd<MatrixType, ComputationOptions>::run(smallSvd, rows, cols, computationOptions);
m_computed = MatrixXr::Zero(m_diagSize + 1, m_diagSize );
m_computed = MatrixXr::Zero(diagSize() + 1, diagSize() );
m_compU = computeV();
m_compV = computeU();
m_isTranspose = (cols > rows);
@ -307,20 +307,20 @@ void BDCSVD<MatrixType, Options>::allocate(Index rows, Index cols, unsigned int
m_useQrDecomp = !disableQrDecomp && ((rows / kMinAspectRatio > cols) || (cols / kMinAspectRatio > rows));
if (m_useQrDecomp) {
qrDecomp = HouseholderQR<MatrixX>((std::max)(rows, cols), (std::min)(rows, cols));
reducedTriangle = MatrixX(m_diagSize, m_diagSize);
reducedTriangle = MatrixX(diagSize(), diagSize());
}
copyWorkspace = MatrixX(m_isTranspose ? cols : rows, m_isTranspose ? rows : cols);
bid = internal::UpperBidiagonalization<MatrixX>(m_useQrDecomp ? m_diagSize : copyWorkspace.rows(),
m_useQrDecomp ? m_diagSize : copyWorkspace.cols());
bid = internal::UpperBidiagonalization<MatrixX>(m_useQrDecomp ? diagSize() : copyWorkspace.rows(),
m_useQrDecomp ? diagSize() : copyWorkspace.cols());
if (m_compU) m_naiveU = MatrixXr::Zero(m_diagSize + 1, m_diagSize + 1 );
else m_naiveU = MatrixXr::Zero(2, m_diagSize + 1 );
if (m_compU) m_naiveU = MatrixXr::Zero(diagSize() + 1, diagSize() + 1 );
else m_naiveU = MatrixXr::Zero(2, diagSize() + 1 );
if (m_compV) m_naiveV = MatrixXr::Zero(m_diagSize, m_diagSize);
if (m_compV) m_naiveV = MatrixXr::Zero(diagSize(), diagSize());
m_workspace.resize((m_diagSize+1)*(m_diagSize+1)*3);
m_workspaceI.resize(3*m_diagSize);
m_workspace.resize((diagSize()+1)*(diagSize()+1)*3);
m_workspaceI.resize(3*diagSize());
} // end allocate
template <typename MatrixType, int Options>
@ -369,7 +369,7 @@ BDCSVD<MatrixType, Options>& BDCSVD<MatrixType, Options>::compute_impl(const Mat
// bidiagonalize the input matrix directly.
if (m_useQrDecomp) {
qrDecomp.compute(copyWorkspace);
reducedTriangle = qrDecomp.matrixQR().topRows(m_diagSize);
reducedTriangle = qrDecomp.matrixQR().topRows(diagSize());
reducedTriangle.template triangularView<StrictlyLower>().setZero();
bid.compute(reducedTriangle);
} else {
@ -380,26 +380,26 @@ BDCSVD<MatrixType, Options>& BDCSVD<MatrixType, Options>::compute_impl(const Mat
m_naiveU.setZero();
m_naiveV.setZero();
// FIXME this line involves a temporary matrix
m_computed.topRows(m_diagSize) = bid.bidiagonal().toDenseMatrix().transpose();
m_computed.topRows(diagSize()) = bid.bidiagonal().toDenseMatrix().transpose();
m_computed.template bottomRows<1>().setZero();
divide(0, m_diagSize - 1, 0, 0, 0);
divide(0, 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++)
for (int i=0; i<diagSize(); i++)
{
RealScalar a = abs(m_computed.coeff(i, i));
m_singularValues.coeffRef(i) = a * scale;
if (a<considerZero)
{
m_nonzeroSingularValues = i;
m_singularValues.tail(m_diagSize - i - 1).setZero();
m_singularValues.tail(diagSize() - i - 1).setZero();
break;
}
else if (i == m_diagSize - 1)
else if (i == diagSize() - 1)
{
m_nonzeroSingularValues = i + 1;
break;
@ -426,20 +426,20 @@ void BDCSVD<MatrixType, Options>::copyUV(const HouseholderU& householderU, const
// Note exchange of U and V: m_matrixU is set from m_naiveV and vice versa
if (computeU())
{
Index Ucols = m_computeThinU ? m_diagSize : rows();
Index Ucols = m_computeThinU ? diagSize() : rows();
m_matrixU = MatrixX::Identity(rows(), Ucols);
m_matrixU.topLeftCorner(m_diagSize, m_diagSize) = naiveV.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
m_matrixU.topLeftCorner(diagSize(), diagSize()) = naiveV.template cast<Scalar>().topLeftCorner(diagSize(), diagSize());
// FIXME the following conditionals involve temporary buffers
if (m_useQrDecomp) m_matrixU.topLeftCorner(householderU.cols(), m_diagSize).applyOnTheLeft(householderU);
if (m_useQrDecomp) m_matrixU.topLeftCorner(householderU.cols(), diagSize()).applyOnTheLeft(householderU);
else m_matrixU.applyOnTheLeft(householderU);
}
if (computeV())
{
Index Vcols = m_computeThinV ? m_diagSize : cols();
Index Vcols = m_computeThinV ? diagSize() : cols();
m_matrixV = MatrixX::Identity(cols(), Vcols);
m_matrixV.topLeftCorner(m_diagSize, m_diagSize) = naiveU.template cast<Scalar>().topLeftCorner(m_diagSize, m_diagSize);
m_matrixV.topLeftCorner(diagSize(), diagSize()) = naiveU.template cast<Scalar>().topLeftCorner(diagSize(), diagSize());
// FIXME the following conditionals involve temporary buffers
if (m_useQrDecomp) m_matrixV.topLeftCorner(householderV.cols(), m_diagSize).applyOnTheLeft(householderV);
if (m_useQrDecomp) m_matrixV.topLeftCorner(householderV.cols(), diagSize()).applyOnTheLeft(householderV);
else m_matrixV.applyOnTheLeft(householderV);
}
}

View File

@ -67,8 +67,8 @@ public:
// prepare arguments to ?gesdd
const lapack_int matrix_order = lapack_storage_of(matrix);
const char jobz = (SVD::m_computeFullU || SVD::m_computeFullV) ? 'A' : (SVD::m_computeThinU || SVD::m_computeThinV) ? 'S' : 'N';
const lapack_int u_cols = (jobz == 'A') ? to_lapack(SVD::m_rows) : (jobz == 'S') ? to_lapack(SVD::m_diagSize) : 1;
const lapack_int vt_rows = (jobz == 'A') ? to_lapack(SVD::m_cols) : (jobz == 'S') ? to_lapack(SVD::m_diagSize) : 1;
const lapack_int u_cols = (jobz == 'A') ? to_lapack(SVD::rows()) : (jobz == 'S') ? to_lapack(SVD::diagSize()) : 1;
const lapack_int vt_rows = (jobz == 'A') ? to_lapack(SVD::cols()) : (jobz == 'S') ? to_lapack(SVD::diagSize()) : 1;
lapack_int ldu, ldvt;
Scalar *u, *vt, dummy;
MatrixType localU;
@ -76,20 +76,20 @@ public:
ldu = to_lapack(SVD::m_matrixU.outerStride());
u = SVD::m_matrixU.data();
} else if (SVD::computeV()) {
localU.resize(SVD::m_rows, u_cols);
localU.resize(SVD::rows(), u_cols);
ldu = to_lapack(localU.outerStride());
u = localU.data();
} else { ldu=1; u=&dummy; }
MatrixType localV;
if (SVD::computeU() || SVD::computeV()) {
localV.resize(vt_rows, SVD::m_cols);
localV.resize(vt_rows, SVD::cols());
ldvt = to_lapack(localV.outerStride());
vt = localV.data();
} else { ldvt=1; vt=&dummy; }
MatrixType temp; temp = matrix;
// actual call to ?gesdd
lapack_int info = gesdd( matrix_order, jobz, to_lapack(SVD::m_rows), to_lapack(SVD::m_cols),
lapack_int info = gesdd( matrix_order, jobz, to_lapack(SVD::rows()), to_lapack(SVD::cols()),
to_lapack(temp.data()), to_lapack(temp.outerStride()), (RealScalar*)SVD::m_singularValues.data(),
to_lapack(u), ldu, to_lapack(vt), ldvt);

View File

@ -519,7 +519,7 @@ class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_> > {
typedef typename Base::Scalar Scalar;
typedef typename Base::RealScalar RealScalar;
typedef typename Base::Index Index;
enum {
enum : int {
Options = Options_,
QRPreconditioner = internal::get_qr_preconditioner(Options),
RowsAtCompileTime = Base::RowsAtCompileTime,
@ -625,6 +625,7 @@ class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_> > {
using Base::computeV;
using Base::rows;
using Base::cols;
using Base::diagSize;
using Base::rank;
private:
@ -632,13 +633,11 @@ class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_> > {
JacobiSVD& compute_impl(const MatrixType& matrix, unsigned int computationOptions);
protected:
using Base::m_cols;
using Base::m_computationOptions;
using Base::m_computeFullU;
using Base::m_computeFullV;
using Base::m_computeThinU;
using Base::m_computeThinV;
using Base::m_diagSize;
using Base::m_info;
using Base::m_isAllocated;
using Base::m_isInitialized;
@ -646,7 +645,6 @@ class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_> > {
using Base::m_matrixV;
using Base::m_nonzeroSingularValues;
using Base::m_prescribedThreshold;
using Base::m_rows;
using Base::m_singularValues;
using Base::m_usePrescribedThreshold;
using Base::ShouldComputeThinU;
@ -671,18 +669,18 @@ class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options_> > {
};
template <typename MatrixType, int Options>
void JacobiSVD<MatrixType, Options>::allocate(Index rows, Index cols, unsigned int computationOptions) {
if (Base::allocate(rows, cols, computationOptions)) return;
void JacobiSVD<MatrixType, Options>::allocate(Index rows_, Index cols_, unsigned int computationOptions_) {
if (Base::allocate(rows_, cols_, computationOptions_)) return;
eigen_assert(!(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)) &&
!(ShouldComputeThinU && int(QRPreconditioner) == int(FullPivHouseholderQRPreconditioner)) &&
"JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. "
"Use the ColPivHouseholderQR preconditioner instead.");
m_workMatrix.resize(m_diagSize, m_diagSize);
if(m_cols>m_rows) m_qr_precond_morecols.allocate(*this);
if(m_rows>m_cols) m_qr_precond_morerows.allocate(*this);
if(m_rows!=m_cols) m_scaledMatrix.resize(rows,cols);
m_workMatrix.resize(diagSize(), diagSize());
if(cols()>rows()) m_qr_precond_morecols.allocate(*this);
if(rows()>cols()) m_qr_precond_morerows.allocate(*this);
if(rows()!=cols()) m_scaledMatrix.resize(rows(),cols());
}
template <typename MatrixType, int Options>
@ -711,7 +709,7 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(con
/*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */
if(m_rows!=m_cols)
if(rows() != cols())
{
m_scaledMatrix = matrix / scale;
m_qr_precond_morecols.run(*this, m_scaledMatrix);
@ -719,11 +717,11 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(con
}
else
{
m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize) / scale;
if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows);
if(m_computeThinU) m_matrixU.setIdentity(m_rows,m_diagSize);
if(m_computeFullV) m_matrixV.setIdentity(m_cols,m_cols);
if(m_computeThinV) m_matrixV.setIdentity(m_cols, m_diagSize);
m_workMatrix = matrix.template topLeftCorner<DiagSizeAtCompileTime,DiagSizeAtCompileTime>(diagSize(),diagSize()) / scale;
if(m_computeFullU) m_matrixU.setIdentity(rows(),rows());
if(m_computeThinU) m_matrixU.setIdentity(rows(),diagSize());
if(m_computeFullV) m_matrixV.setIdentity(cols(),cols());
if(m_computeThinV) m_matrixV.setIdentity(cols(),diagSize());
}
/*** step 2. The main Jacobi SVD iteration. ***/
@ -736,7 +734,7 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(con
// do a sweep: for all index pairs (p,q), perform SVD of the corresponding 2x2 sub-matrix
for(Index p = 1; p < m_diagSize; ++p)
for(Index p = 1; p < diagSize(); ++p)
{
for(Index q = 0; q < p; ++q)
{
@ -771,7 +769,7 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(con
/*** step 3. The work matrix is now diagonal, so ensure it's positive so its diagonal entries are the singular values ***/
for(Index i = 0; i < m_diagSize; ++i)
for(Index i = 0; i < diagSize(); ++i)
{
// For a complex matrix, some diagonal coefficients might note have been
// treated by svd_precondition_2x2_block_to_be_real, and the imaginary part
@ -795,11 +793,11 @@ JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::compute_impl(con
/*** step 4. Sort singular values in descending order and compute the number of nonzero singular values ***/
m_nonzeroSingularValues = m_diagSize;
for(Index i = 0; i < m_diagSize; i++)
m_nonzeroSingularValues = diagSize();
for(Index i = 0; i < diagSize(); i++)
{
Index pos;
RealScalar maxRemainingSingularValue = m_singularValues.tail(m_diagSize-i).maxCoeff(&pos);
RealScalar maxRemainingSingularValue = m_singularValues.tail(diagSize()-i).maxCoeff(&pos);
if(numext::is_exactly_zero(maxRemainingSingularValue))
{
m_nonzeroSingularValues = i;

View File

@ -51,7 +51,7 @@ JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTION
allocate(matrix.rows(), matrix.cols(), computationOptions); \
\
/*const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon();*/ \
m_nonzeroSingularValues = m_diagSize; \
m_nonzeroSingularValues = diagSize(); \
\
lapack_int lda = internal::convert_index<lapack_int>(matrix.outerStride()), ldu, ldvt; \
lapack_int matrix_order = LAPACKE_COLROW; \
@ -64,15 +64,15 @@ JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTION
u = (LAPACKE_TYPE*)m_matrixU.data(); \
} else { ldu=1; u=&dummy; }\
MatrixType localV; \
lapack_int vt_rows = (m_computeFullV) ? internal::convert_index<lapack_int>(m_cols) : (m_computeThinV) ? internal::convert_index<lapack_int>(m_diagSize) : 1; \
lapack_int vt_rows = (m_computeFullV) ? internal::convert_index<lapack_int>(cols()) : (m_computeThinV) ? internal::convert_index<lapack_int>(diagSize()) : 1; \
if (computeV()) { \
localV.resize(vt_rows, m_cols); \
localV.resize(vt_rows, cols()); \
ldvt = internal::convert_index<lapack_int>(localV.outerStride()); \
vt = (LAPACKE_TYPE*)localV.data(); \
} else { ldvt=1; vt=&dummy; }\
Matrix<LAPACKE_RTYPE, Dynamic, Dynamic> superb; superb.resize(m_diagSize, 1); \
Matrix<LAPACKE_RTYPE, Dynamic, Dynamic> superb; superb.resize(diagSize(), 1); \
MatrixType m_temp; m_temp = matrix; \
lapack_int info = LAPACKE_##LAPACKE_PREFIX##gesvd( matrix_order, jobu, jobvt, internal::convert_index<lapack_int>(m_rows), internal::convert_index<lapack_int>(m_cols), (LAPACKE_TYPE*)m_temp.data(), lda, (LAPACKE_RTYPE*)m_singularValues.data(), u, ldu, vt, ldvt, superb.data()); \
lapack_int info = LAPACKE_##LAPACKE_PREFIX##gesvd( matrix_order, jobu, jobvt, internal::convert_index<lapack_int>(rows()), internal::convert_index<lapack_int>(cols()), (LAPACKE_TYPE*)m_temp.data(), lda, (LAPACKE_RTYPE*)m_singularValues.data(), u, ldu, vt, ldvt, superb.data()); \
/* Check the result of the LAPACK call */ \
if (info < 0 || !m_singularValues.allFinite()) { \
m_info = InvalidInput; \
@ -82,7 +82,7 @@ JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTION
m_info = Success; \
if (computeV()) m_matrixV = localV.adjoint(); \
} \
/* for(int i=0;i<m_diagSize;i++) if (m_singularValues.coeffRef(i) < precision) { m_nonzeroSingularValues--; m_singularValues.coeffRef(i)=RealScalar(0);}*/ \
/* for(int i=0;i<diagSize();i++) if (m_singularValues.coeffRef(i) < precision) { m_nonzeroSingularValues--; m_singularValues.coeffRef(i)=RealScalar(0);}*/ \
m_isInitialized = true; \
return *this; \
}

View File

@ -283,8 +283,9 @@ public:
/** \returns true if \a V (full or thin) is asked for in this SVD decomposition */
inline bool computeV() const { return m_computeFullV || m_computeThinV; }
inline Index rows() const { return m_rows; }
inline Index cols() const { return m_cols; }
inline Index rows() const { return m_rows.value(); }
inline Index cols() const { return m_cols.value(); }
inline Index diagSize() const { return m_diagSize.value(); }
#ifdef EIGEN_PARSED_BY_DOXYGEN
/** \returns a (least squares) solution of \f$ A x = b \f$ using the current SVD decomposition of A.
@ -348,7 +349,10 @@ protected:
bool m_computeFullU, m_computeThinU;
bool m_computeFullV, m_computeThinV;
unsigned int m_computationOptions;
Index m_nonzeroSingularValues, m_rows, m_cols, m_diagSize;
Index m_nonzeroSingularValues;
internal::variable_if_dynamic<Index, RowsAtCompileTime> m_rows;
internal::variable_if_dynamic<Index, ColsAtCompileTime> m_cols;
internal::variable_if_dynamic<Index, DiagSizeAtCompileTime> m_diagSize;
RealScalar m_prescribedThreshold;
/** \brief Default Constructor.
@ -369,9 +373,9 @@ protected:
m_computeThinV(false),
m_computationOptions(0),
m_nonzeroSingularValues(0),
m_rows(-1),
m_cols(-1),
m_diagSize(0),
m_rows(RowsAtCompileTime),
m_cols(ColsAtCompileTime),
m_diagSize(DiagSizeAtCompileTime),
m_prescribedThreshold(0) {}
};
@ -411,15 +415,15 @@ bool SVDBase<Derived>::allocate(Index rows, Index cols, unsigned int computation
eigen_assert(rows >= 0 && cols >= 0);
if (m_isAllocated &&
rows == m_rows &&
cols == m_cols &&
rows == m_rows.value() &&
cols == m_cols.value() &&
computationOptions == m_computationOptions)
{
return true;
}
m_rows = rows;
m_cols = cols;
m_rows.setValue(rows);
m_cols.setValue(cols);
m_info = Success;
m_isInitialized = false;
m_isAllocated = true;
@ -432,12 +436,12 @@ bool SVDBase<Derived>::allocate(Index rows, Index cols, unsigned int computation
eigen_assert(!(m_computeFullU && m_computeThinU) && "SVDBase: you can't ask for both full and thin U");
eigen_assert(!(m_computeFullV && m_computeThinV) && "SVDBase: you can't ask for both full and thin V");
m_diagSize = (std::min)(m_rows, m_cols);
m_singularValues.resize(m_diagSize);
m_diagSize.setValue(numext::mini(m_rows.value(), m_cols.value()));
m_singularValues.resize(m_diagSize.value());
if(RowsAtCompileTime==Dynamic)
m_matrixU.resize(m_rows, m_computeFullU ? m_rows : m_computeThinU ? m_diagSize : 0);
m_matrixU.resize(m_rows.value(), m_computeFullU ? m_rows.value() : m_computeThinU ? m_diagSize.value() : 0);
if(ColsAtCompileTime==Dynamic)
m_matrixV.resize(m_cols, m_computeFullV ? m_cols : m_computeThinV ? m_diagSize : 0);
m_matrixV.resize(m_cols.value(), m_computeFullV ? m_cols.value() : m_computeThinV ? m_diagSize.value() : 0);
return false;
}

View File

@ -69,8 +69,13 @@ void compare_bdc_jacobi_instance(bool structure_as_m, int algoswap = 16)
}
template <typename MatrixType>
void bdcsvd_all_options(const MatrixType& input = MatrixType()) {
svd_option_checks<MatrixType, 0>(input);
void bdcsvd_thin_options(const MatrixType& input = MatrixType()) {
svd_thin_option_checks<MatrixType, 0>(input);
}
template <typename MatrixType>
void bdcsvd_full_options(const MatrixType& input = MatrixType()) {
svd_option_checks_full_only<MatrixType, 0>(input);
}
template <typename MatrixType>
@ -82,13 +87,15 @@ void bdcsvd_verify_assert(const MatrixType& input = MatrixType()) {
EIGEN_DECLARE_TEST(bdcsvd)
{
CALL_SUBTEST_1((bdcsvd_verify_assert<Matrix3f>()));
CALL_SUBTEST_1((bdcsvd_verify_assert<Matrix4d>()));
CALL_SUBTEST_2((bdcsvd_verify_assert<Matrix<float, 10, 7>>()));
CALL_SUBTEST_2((bdcsvd_verify_assert<Matrix<float, 7, 10>>()));
CALL_SUBTEST_3((bdcsvd_verify_assert<Matrix<std::complex<double>, 6, 9>>()));
CALL_SUBTEST_2((bdcsvd_verify_assert<Matrix4d>()));
CALL_SUBTEST_3((bdcsvd_verify_assert<Matrix<float, 10, 7>>()));
CALL_SUBTEST_4((bdcsvd_verify_assert<Matrix<float, 7, 10>>()));
CALL_SUBTEST_5((bdcsvd_verify_assert<Matrix<std::complex<double>, 6, 9>>()));
CALL_SUBTEST_4((svd_all_trivial_2x2(bdcsvd_all_options<Matrix2cd>)));
CALL_SUBTEST_5((svd_all_trivial_2x2(bdcsvd_all_options<Matrix2d>)));
CALL_SUBTEST_6((svd_all_trivial_2x2(bdcsvd_thin_options<Matrix2cd>)));
CALL_SUBTEST_7((svd_all_trivial_2x2(bdcsvd_full_options<Matrix2cd>)));
CALL_SUBTEST_8((svd_all_trivial_2x2(bdcsvd_thin_options<Matrix2d>)));
CALL_SUBTEST_9((svd_all_trivial_2x2(bdcsvd_full_options<Matrix2d>)));
for (int i = 0; i < g_repeat; i++) {
int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2),
@ -97,58 +104,68 @@ EIGEN_DECLARE_TEST(bdcsvd)
TEST_SET_BUT_UNUSED_VARIABLE(r)
TEST_SET_BUT_UNUSED_VARIABLE(c)
CALL_SUBTEST_6((compare_bdc_jacobi<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_7((compare_bdc_jacobi<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_8((compare_bdc_jacobi<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_10((compare_bdc_jacobi<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_11((compare_bdc_jacobi<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_12((compare_bdc_jacobi<MatrixXcd>(MatrixXcd(r, c))));
// Test on inf/nan matrix
CALL_SUBTEST_9((svd_inf_nan<MatrixXf>()));
CALL_SUBTEST_10((svd_inf_nan<MatrixXd>()));
CALL_SUBTEST_13((svd_inf_nan<MatrixXf>()));
CALL_SUBTEST_14((svd_inf_nan<MatrixXd>()));
// Verify some computations using all combinations of the Options template parameter.
CALL_SUBTEST_11((bdcsvd_all_options<Matrix3f>()));
CALL_SUBTEST_12((bdcsvd_all_options<Matrix<float, 2, 3>>()));
CALL_SUBTEST_13((bdcsvd_all_options<MatrixXd>(MatrixXd(20, 17))));
CALL_SUBTEST_14((bdcsvd_all_options<MatrixXd>(MatrixXd(17, 20))));
CALL_SUBTEST_15((bdcsvd_all_options<Matrix<double, Dynamic, 15>>(Matrix<double, Dynamic, 15>(r, 15))));
CALL_SUBTEST_16((bdcsvd_all_options<Matrix<double, 13, Dynamic>>(Matrix<double, 13, Dynamic>(13, c))));
CALL_SUBTEST_17((bdcsvd_all_options<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_18((bdcsvd_all_options<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_19((bdcsvd_all_options<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_20((bdcsvd_all_options<Matrix<double, Dynamic, Dynamic, RowMajor>>(Matrix<double, Dynamic, Dynamic, RowMajor>(20, 27))));
CALL_SUBTEST_21((bdcsvd_all_options<Matrix<double, Dynamic, Dynamic, RowMajor>>(Matrix<double, Dynamic, Dynamic, RowMajor>(27, 20))));
CALL_SUBTEST_22((
CALL_SUBTEST_15((bdcsvd_thin_options<Matrix3f>()));
CALL_SUBTEST_16((bdcsvd_full_options<Matrix3f>()));
CALL_SUBTEST_17((bdcsvd_thin_options<Matrix<float, 2, 3>>()));
CALL_SUBTEST_18((bdcsvd_full_options<Matrix<float, 2, 3>>()));
CALL_SUBTEST_19((bdcsvd_thin_options<MatrixXd>(MatrixXd(20, 17))));
CALL_SUBTEST_20((bdcsvd_full_options<MatrixXd>(MatrixXd(20, 17))));
CALL_SUBTEST_21((bdcsvd_thin_options<MatrixXd>(MatrixXd(17, 20))));
CALL_SUBTEST_22((bdcsvd_full_options<MatrixXd>(MatrixXd(17, 20))));
CALL_SUBTEST_23((bdcsvd_thin_options<Matrix<double, Dynamic, 15>>(Matrix<double, Dynamic, 15>(r, 15))));
CALL_SUBTEST_24((bdcsvd_full_options<Matrix<double, Dynamic, 15>>(Matrix<double, Dynamic, 15>(r, 15))));
CALL_SUBTEST_25((bdcsvd_thin_options<Matrix<double, 13, Dynamic>>(Matrix<double, 13, Dynamic>(13, c))));
CALL_SUBTEST_26((bdcsvd_full_options<Matrix<double, 13, Dynamic>>(Matrix<double, 13, Dynamic>(13, c))));
CALL_SUBTEST_27((bdcsvd_thin_options<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_28((bdcsvd_full_options<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_29((bdcsvd_thin_options<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_30((bdcsvd_full_options<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_31((bdcsvd_thin_options<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_32((bdcsvd_full_options<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_33((bdcsvd_thin_options<Matrix<double, Dynamic, Dynamic, RowMajor>>(Matrix<double, Dynamic, Dynamic, RowMajor>(20, 27))));
CALL_SUBTEST_34((bdcsvd_full_options<Matrix<double, Dynamic, Dynamic, RowMajor>>(Matrix<double, Dynamic, Dynamic, RowMajor>(20, 27))));
CALL_SUBTEST_35((bdcsvd_thin_options<Matrix<double, Dynamic, Dynamic, RowMajor>>(Matrix<double, Dynamic, Dynamic, RowMajor>(27, 20))));
CALL_SUBTEST_36((bdcsvd_full_options<Matrix<double, Dynamic, Dynamic, RowMajor>>(Matrix<double, Dynamic, Dynamic, RowMajor>(27, 20))));
CALL_SUBTEST_37((
svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, ColMajor, 20, 35>, ColPivHouseholderQRPreconditioner>(
r, c)));
CALL_SUBTEST_22(
CALL_SUBTEST_38(
(svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, ColMajor, 35, 20>, HouseholderQRPreconditioner>(r,
c)));
CALL_SUBTEST_22((
CALL_SUBTEST_39((
svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, RowMajor, 20, 35>, ColPivHouseholderQRPreconditioner>(
r, c)));
CALL_SUBTEST_22(
CALL_SUBTEST_40(
(svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, RowMajor, 35, 20>, HouseholderQRPreconditioner>(r,
c)));
}
// test matrixbase method
CALL_SUBTEST_23(( bdcsvd_method<Matrix2cd>() ));
CALL_SUBTEST_23(( bdcsvd_method<Matrix3f>() ));
CALL_SUBTEST_41(( bdcsvd_method<Matrix2cd>() ));
CALL_SUBTEST_42(( bdcsvd_method<Matrix3f>() ));
// Test problem size constructors
CALL_SUBTEST_24( BDCSVD<MatrixXf>(10,10) );
CALL_SUBTEST_43( BDCSVD<MatrixXf>(10,10) );
// Check that preallocation avoids subsequent mallocs
// Disabled because not supported by BDCSVD
// CALL_SUBTEST_9( svd_preallocate<void>() );
CALL_SUBTEST_25( svd_underoverflow<void>() );
CALL_SUBTEST_44( svd_underoverflow<void>() );
// Without total deflation issues.
CALL_SUBTEST_26(( compare_bdc_jacobi_instance(true) ));
CALL_SUBTEST_26(( compare_bdc_jacobi_instance(false) ));
CALL_SUBTEST_45(( compare_bdc_jacobi_instance(true) ));
CALL_SUBTEST_46(( compare_bdc_jacobi_instance(false) ));
// With total deflation issues before, when it shouldn't be triggered.
CALL_SUBTEST_27(( compare_bdc_jacobi_instance(true, 3) ));
CALL_SUBTEST_27(( compare_bdc_jacobi_instance(false, 3) ));
CALL_SUBTEST_47(( compare_bdc_jacobi_instance(true, 3) ));
CALL_SUBTEST_48(( compare_bdc_jacobi_instance(false, 3) ));
}

View File

@ -35,12 +35,23 @@ void jacobisvd_method()
}
template <typename MatrixType>
void jacobisvd_all_options(const MatrixType& input = MatrixType()) {
void jacobisvd_thin_options(const MatrixType& input = MatrixType()) {
MatrixType m(input.rows(), input.cols());
svd_fill_random(m);
svd_option_checks<MatrixType, 0>(m);
svd_option_checks<MatrixType, ColPivHouseholderQRPreconditioner>(m);
svd_option_checks<MatrixType, HouseholderQRPreconditioner>(m);
svd_thin_option_checks<MatrixType, 0>(m);
svd_thin_option_checks<MatrixType, ColPivHouseholderQRPreconditioner>(m);
svd_thin_option_checks<MatrixType, HouseholderQRPreconditioner>(m);
}
template <typename MatrixType>
void jacobisvd_full_options(const MatrixType& input = MatrixType()) {
MatrixType m(input.rows(), input.cols());
svd_fill_random(m);
svd_option_checks_full_only<MatrixType, 0>(m);
svd_option_checks_full_only<MatrixType, ColPivHouseholderQRPreconditioner>(m);
svd_option_checks_full_only<MatrixType, HouseholderQRPreconditioner>(m);
svd_option_checks_full_only<MatrixType, FullPivHouseholderQRPreconditioner>(
m); // FullPiv only used when computing full unitaries
}
@ -103,18 +114,18 @@ void msvc_workaround()
EIGEN_DECLARE_TEST(jacobisvd)
{
CALL_SUBTEST_1((jacobisvd_verify_inputs<Matrix4d>()));
CALL_SUBTEST_1((jacobisvd_verify_inputs(Matrix<float, 5, Dynamic>(5, 6))));
CALL_SUBTEST_1((jacobisvd_verify_inputs<Matrix<std::complex<double>, 7, 5>>()));
CALL_SUBTEST_2((jacobisvd_verify_inputs(Matrix<float, 5, Dynamic>(5, 6))));
CALL_SUBTEST_3((jacobisvd_verify_inputs<Matrix<std::complex<double>, 7, 5>>()));
CALL_SUBTEST_2((jacobisvd_verify_assert<Matrix3f>()));
CALL_SUBTEST_2((jacobisvd_verify_assert<Matrix4d>()));
CALL_SUBTEST_2((jacobisvd_verify_assert<Matrix<float, 10, 12>>()));
CALL_SUBTEST_2((jacobisvd_verify_assert<Matrix<float, 12, 10>>()));
CALL_SUBTEST_2((jacobisvd_verify_assert<MatrixXf>(MatrixXf(10, 12))));
CALL_SUBTEST_2((jacobisvd_verify_assert<MatrixXcd>(MatrixXcd(7, 5))));
CALL_SUBTEST_4((jacobisvd_verify_assert<Matrix3f>()));
CALL_SUBTEST_5((jacobisvd_verify_assert<Matrix4d>()));
CALL_SUBTEST_6((jacobisvd_verify_assert<Matrix<float, 10, 12>>()));
CALL_SUBTEST_7((jacobisvd_verify_assert<Matrix<float, 12, 10>>()));
CALL_SUBTEST_8((jacobisvd_verify_assert<MatrixXf>(MatrixXf(10, 12))));
CALL_SUBTEST_9((jacobisvd_verify_assert<MatrixXcd>(MatrixXcd(7, 5))));
CALL_SUBTEST_3(svd_all_trivial_2x2(jacobisvd_all_options<Matrix2cd>));
CALL_SUBTEST_4(svd_all_trivial_2x2(jacobisvd_all_options<Matrix2d>));
CALL_SUBTEST_10(svd_all_trivial_2x2(jacobisvd_thin_options<Matrix2cd>));
CALL_SUBTEST_11(svd_all_trivial_2x2(jacobisvd_thin_options<Matrix2d>));
for (int i = 0; i < g_repeat; i++) {
int r = internal::random<int>(1, 30),
@ -123,64 +134,83 @@ EIGEN_DECLARE_TEST(jacobisvd)
TEST_SET_BUT_UNUSED_VARIABLE(r)
TEST_SET_BUT_UNUSED_VARIABLE(c)
CALL_SUBTEST_5((jacobisvd_all_options<Matrix3f>()));
CALL_SUBTEST_6((jacobisvd_all_options<Matrix4d>()));
CALL_SUBTEST_7((jacobisvd_all_options<Matrix<float, 2, 3>>()));
CALL_SUBTEST_8((jacobisvd_all_options<Matrix<double, 4, 7>>()));
CALL_SUBTEST_9((jacobisvd_all_options<Matrix<double, 7, 4>>()));
CALL_SUBTEST_10((jacobisvd_all_options<Matrix<double, Dynamic, 5>>(Matrix<double, Dynamic, 5>(r, 5))));
CALL_SUBTEST_11((jacobisvd_all_options<Matrix<double, 5, Dynamic>>(Matrix<double, 5, Dynamic>(5, c))));
CALL_SUBTEST_12((jacobisvd_all_options<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_13((jacobisvd_all_options<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_14((jacobisvd_all_options<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_15((jacobisvd_all_options<Matrix<double, 5, 7, RowMajor>>()));
CALL_SUBTEST_16((jacobisvd_all_options<Matrix<double, 7, 5, RowMajor>>()));
CALL_SUBTEST_12((jacobisvd_thin_options<Matrix3f>()));
CALL_SUBTEST_13((jacobisvd_full_options<Matrix3f>()));
CALL_SUBTEST_14((jacobisvd_thin_options<Matrix4d>()));
CALL_SUBTEST_15((jacobisvd_full_options<Matrix4d>()));
CALL_SUBTEST_16((jacobisvd_thin_options<Matrix<float, 2, 3>>()));
CALL_SUBTEST_17((jacobisvd_full_options<Matrix<float, 2, 3>>()));
CALL_SUBTEST_18((jacobisvd_thin_options<Matrix<double, 4, 7>>()));
CALL_SUBTEST_19((jacobisvd_full_options<Matrix<double, 4, 7>>()));
CALL_SUBTEST_20((jacobisvd_thin_options<Matrix<double, 7, 4>>()));
CALL_SUBTEST_21((jacobisvd_full_options<Matrix<double, 7, 4>>()));
CALL_SUBTEST_22((jacobisvd_thin_options<Matrix<double, Dynamic, 5>>(Matrix<double, Dynamic, 5>(r, 5))));
CALL_SUBTEST_23((jacobisvd_full_options<Matrix<double, Dynamic, 5>>(Matrix<double, Dynamic, 5>(r, 5))));
CALL_SUBTEST_24((jacobisvd_thin_options<Matrix<double, 5, Dynamic>>(Matrix<double, 5, Dynamic>(5, c))));
CALL_SUBTEST_25((jacobisvd_full_options<Matrix<double, 5, Dynamic>>(Matrix<double, 5, Dynamic>(5, c))));
CALL_SUBTEST_26((jacobisvd_thin_options<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_27((jacobisvd_full_options<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_28((jacobisvd_thin_options<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_29((jacobisvd_full_options<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_30((jacobisvd_thin_options<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_31((jacobisvd_full_options<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_32((jacobisvd_thin_options<Matrix<double, 5, 7, RowMajor>>()));
CALL_SUBTEST_33((jacobisvd_full_options<Matrix<double, 5, 7, RowMajor>>()));
CALL_SUBTEST_34((jacobisvd_thin_options<Matrix<double, 7, 5, RowMajor>>()));
CALL_SUBTEST_35((jacobisvd_full_options<Matrix<double, 7, 5, RowMajor>>()));
MatrixXcd noQRTest = MatrixXcd(r, r);
svd_fill_random(noQRTest);
CALL_SUBTEST_17((svd_option_checks<MatrixXcd, NoQRPreconditioner>(noQRTest)));
CALL_SUBTEST_36((svd_thin_option_checks<MatrixXcd, NoQRPreconditioner>(noQRTest)));
CALL_SUBTEST_36((svd_option_checks_full_only<MatrixXcd, NoQRPreconditioner>(noQRTest)));
CALL_SUBTEST_18((
CALL_SUBTEST_37((
svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, ColMajor, 13, 15>, ColPivHouseholderQRPreconditioner>(
r, c)));
CALL_SUBTEST_18(
CALL_SUBTEST_38(
(svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, ColMajor, 15, 13>, HouseholderQRPreconditioner>(r,
c)));
CALL_SUBTEST_18((
CALL_SUBTEST_39((
svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, RowMajor, 13, 15>, ColPivHouseholderQRPreconditioner>(
r, c)));
CALL_SUBTEST_18(
CALL_SUBTEST_40(
(svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, RowMajor, 15, 13>, HouseholderQRPreconditioner>(r,
c)));
// Test on inf/nan matrix
CALL_SUBTEST_19((svd_inf_nan<MatrixXf>()));
CALL_SUBTEST_19((svd_inf_nan<MatrixXd>()));
CALL_SUBTEST_41((svd_inf_nan<MatrixXf>()));
CALL_SUBTEST_42((svd_inf_nan<MatrixXd>()));
CALL_SUBTEST_20((jacobisvd_verify_assert<Matrix<double, 6, 1>>()));
CALL_SUBTEST_20((jacobisvd_verify_assert<Matrix<double, 1, 6>>()));
CALL_SUBTEST_20((jacobisvd_verify_assert<Matrix<double, Dynamic, 1>>(Matrix<double, Dynamic, 1>(r))));
CALL_SUBTEST_20((jacobisvd_verify_assert<Matrix<double, 1, Dynamic>>(Matrix<double, 1, Dynamic>(c))));
CALL_SUBTEST_43((jacobisvd_verify_assert<Matrix<double, 6, 1>>()));
CALL_SUBTEST_44((jacobisvd_verify_assert<Matrix<double, 1, 6>>()));
CALL_SUBTEST_45((jacobisvd_verify_assert<Matrix<double, Dynamic, 1>>(Matrix<double, Dynamic, 1>(r))));
CALL_SUBTEST_46((jacobisvd_verify_assert<Matrix<double, 1, Dynamic>>(Matrix<double, 1, Dynamic>(c))));
}
CALL_SUBTEST_21((jacobisvd_all_options<MatrixXd>(
CALL_SUBTEST_47((jacobisvd_thin_options<MatrixXd>(
MatrixXd(internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 2),
internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 2)))));
CALL_SUBTEST_22((jacobisvd_all_options<MatrixXcd>(
CALL_SUBTEST_48((jacobisvd_full_options<MatrixXd>(
MatrixXd(internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 2),
internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 2)))));
CALL_SUBTEST_49((jacobisvd_thin_options<MatrixXcd>(
MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 3),
internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 3)))));
CALL_SUBTEST_50((jacobisvd_full_options<MatrixXcd>(
MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 3),
internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 3)))));
// test matrixbase method
CALL_SUBTEST_23(( jacobisvd_method<Matrix2cd>() ));
CALL_SUBTEST_23(( jacobisvd_method<Matrix3f>() ));
CALL_SUBTEST_51(( jacobisvd_method<Matrix2cd>() ));
CALL_SUBTEST_52(( jacobisvd_method<Matrix3f>() ));
// Test problem size constructors
CALL_SUBTEST_24( JacobiSVD<MatrixXf>(10,10) );
CALL_SUBTEST_53( JacobiSVD<MatrixXf>(10,10) );
// Check that preallocation avoids subsequent mallocs
CALL_SUBTEST_25( svd_preallocate<void>() );
CALL_SUBTEST_54( svd_preallocate<void>() );
CALL_SUBTEST_26( svd_underoverflow<void>() );
CALL_SUBTEST_55( svd_underoverflow<void>() );
msvc_workaround();
}

View File

@ -484,18 +484,15 @@ void svd_compute_checks(const MatrixType& m) {
}
template <typename MatrixType, int QRPreconditioner = 0>
void svd_option_checks(const MatrixType& input) {
void svd_thin_option_checks(const MatrixType& input) {
MatrixType m(input.rows(), input.cols());
svd_fill_random(m);
svd_compute_checks<MatrixType, QRPreconditioner>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeThinU>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeThinV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeThinU | ComputeThinV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullU>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullU | ComputeFullV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeThinU | ComputeFullV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullU | ComputeThinV>(m);