diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h index 89a411010..fa311bf00 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h @@ -375,8 +375,15 @@ template class MatrixBase /////////// SVD module /////////// - inline JacobiSVD jacobiSvd(unsigned int computationOptions = 0) const; - inline BDCSVD bdcSvd(unsigned int computationOptions = 0) const; + template + inline JacobiSVD jacobiSvd() const; + template + inline JacobiSVD jacobiSvd(unsigned int computationOptions) const; + + template + inline BDCSVD bdcSvd() const; + template + inline BDCSVD bdcSvd(unsigned int computationOptions) const; /////////// Geometry module /////////// diff --git a/Eigen/src/Core/util/Constants.h b/Eigen/src/Core/util/Constants.h index 72fe3307c..16ed58560 100644 --- a/Eigen/src/Core/util/Constants.h +++ b/Eigen/src/Core/util/Constants.h @@ -423,14 +423,14 @@ enum DecompositionOptions { /** \ingroup enums * Possible values for the \p QRPreconditioner template parameter of JacobiSVD. */ enum QRPreconditioners { - /** Do not specify what is to be done if the SVD of a non-square matrix is asked for. */ - NoQRPreconditioner, - /** Use a QR decomposition without pivoting as the first step. */ - HouseholderQRPreconditioner, /** Use a QR decomposition with column pivoting as the first step. */ - ColPivHouseholderQRPreconditioner, + ColPivHouseholderQRPreconditioner = 0x0, + /** Do not specify what is to be done if the SVD of a non-square matrix is asked for. */ + NoQRPreconditioner = 0x40, + /** Use a QR decomposition without pivoting as the first step. */ + HouseholderQRPreconditioner = 0x80, /** Use a QR decomposition with full pivoting as the first step. */ - FullPivHouseholderQRPreconditioner + FullPivHouseholderQRPreconditioner = 0xC0 }; #ifdef Success diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index f3df3728b..051cb303e 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -260,8 +260,8 @@ template class ColPivHouseholderQR; template class FullPivHouseholderQR; template class CompleteOrthogonalDecomposition; template class SVDBase; -template class JacobiSVD; -template class BDCSVD; +template class JacobiSVD; +template class BDCSVD; template class LLT; template class LDLT; template class HouseholderSequence; diff --git a/Eigen/src/Geometry/Hyperplane.h b/Eigen/src/Geometry/Hyperplane.h index b39c9cae9..ad6aae9ef 100644 --- a/Eigen/src/Geometry/Hyperplane.h +++ b/Eigen/src/Geometry/Hyperplane.h @@ -108,7 +108,7 @@ public: if(norm <= v0.norm() * v1.norm() * NumTraits::epsilon()) { Matrix m; m << v0.transpose(), v1.transpose(); - JacobiSVD > svd(m, ComputeFullV); + JacobiSVD, ComputeFullV> svd(m); result.normal() = svd.matrixV().col(2); } else diff --git a/Eigen/src/Geometry/Quaternion.h b/Eigen/src/Geometry/Quaternion.h index e291bf30f..52e8a81f0 100644 --- a/Eigen/src/Geometry/Quaternion.h +++ b/Eigen/src/Geometry/Quaternion.h @@ -651,7 +651,7 @@ EIGEN_DEVICE_FUNC inline Derived& QuaternionBase::setFromTwoVectors(con { c = numext::maxi(c,Scalar(-1)); Matrix m; m << v0.transpose(), v1.transpose(); - JacobiSVD > svd(m, ComputeFullV); + JacobiSVD, ComputeFullV> svd(m); Vector3 axis = svd.matrixV().col(2); Scalar w2 = (Scalar(1)+c)*Scalar(0.5); diff --git a/Eigen/src/Geometry/Transform.h b/Eigen/src/Geometry/Transform.h index 78d51e619..e147d8070 100644 --- a/Eigen/src/Geometry/Transform.h +++ b/Eigen/src/Geometry/Transform.h @@ -1105,7 +1105,7 @@ template EIGEN_DEVICE_FUNC void Transform::computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const { // Note that JacobiSVD is faster than BDCSVD for small matrices. - JacobiSVD svd(linear(), ComputeFullU | ComputeFullV); + JacobiSVD svd(linear()); Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant() < Scalar(0) ? Scalar(-1) : Scalar(1); // so x has absolute value 1 VectorType sv(svd.singularValues()); @@ -1135,7 +1135,7 @@ template EIGEN_DEVICE_FUNC void Transform::computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const { // Note that JacobiSVD is faster than BDCSVD for small matrices. - JacobiSVD svd(linear(), ComputeFullU | ComputeFullV); + JacobiSVD svd(linear()); Scalar x = (svd.matrixU() * svd.matrixV().adjoint()).determinant() < Scalar(0) ? Scalar(-1) : Scalar(1); // so x has absolute value 1 VectorType sv(svd.singularValues()); diff --git a/Eigen/src/Geometry/Umeyama.h b/Eigen/src/Geometry/Umeyama.h index 98adc4bc7..804978720 100644 --- a/Eigen/src/Geometry/Umeyama.h +++ b/Eigen/src/Geometry/Umeyama.h @@ -127,7 +127,7 @@ umeyama(const MatrixBase& src, const MatrixBase& dst, boo // Eq. (38) const MatrixType sigma = one_over_n * dst_demean * src_demean.transpose(); - JacobiSVD svd(sigma, ComputeFullU | ComputeFullV); + JacobiSVD svd(sigma); // Initialize the resulting transformation with an identity matrix... TransformationMatrixType Rt = TransformationMatrixType::Identity(m+1,m+1); diff --git a/Eigen/src/SVD/BDCSVD.h b/Eigen/src/SVD/BDCSVD.h index 7b2064980..a9f7b2939 100644 --- a/Eigen/src/SVD/BDCSVD.h +++ b/Eigen/src/SVD/BDCSVD.h @@ -39,19 +39,32 @@ namespace Eigen { IOFormat bdcsvdfmt(8, 0, ", ", "\n", " [", "]"); #endif -template class BDCSVD; +template +class BDCSVD; namespace internal { -template -struct traits > - : traits -{ +template +struct traits > : svd_traits { typedef MatrixType_ MatrixType; }; -} // end namespace internal +template +struct allocate_small_svd { + static void run(JacobiSVD& smallSvd, Index rows, Index cols, unsigned int computationOptions) { + (void)computationOptions; + smallSvd = JacobiSVD(rows, cols); + } +}; +template +struct allocate_small_svd { + static void run(JacobiSVD& smallSvd, Index rows, Index cols, unsigned int computationOptions) { + smallSvd = JacobiSVD(rows, cols, computationOptions); + } +}; + +} // end namespace internal /** \ingroup SVD_Module * @@ -62,6 +75,11 @@ struct traits > * * \tparam MatrixType_ the type of the matrix of which we are computing the SVD decomposition * + * \tparam Options this optional parameter allows one to specify options for computing unitaries \a U and \a V. + * Possible values are #ComputeThinU, #ComputeThinV, #ComputeFullU, #ComputeFullV. + * It is not possible to request both the thin and full version of \a U or \a V. + * By default, unitaries are not computed. + * * This class first reduces the input matrix to bi-diagonal form using class UpperBidiagonalization, * and then performs a divide-and-conquer diagonalization. Small blocks are diagonalized using class JacobiSVD. * You can control the switching size with the setSwitchSize() method, default is 16. @@ -75,9 +93,8 @@ struct traits > * * \sa class JacobiSVD */ -template -class BDCSVD : public SVDBase > -{ +template +class BDCSVD : public SVDBase > { typedef SVDBase Base; public: @@ -120,61 +137,65 @@ public: BDCSVD() : m_algoswap(16), m_isTranspose(false), m_compU(false), m_compV(false), m_numIters(0) {} + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem size and \a Options template parameter. + * \sa BDCSVD() + */ + BDCSVD(Index rows, Index cols) : m_algoswap(16), m_numIters(0) { + allocate(rows, cols, internal::get_computation_options(Options)); + } /** \brief Default Constructor with memory preallocation * * Like the default constructor but with preallocation of the internal data - * according to the specified problem size. + * according to the specified problem size and the \a computationOptions. + * + * Note: This constructor is deprecated. + * One \b cannot request unitiaries using both the \a Options template parameter + * and the constructor. If possible, prefer using the \a Options template parameter. + * + * \param computationOptions specifification for computing Thin/Full unitaries U/V * \sa BDCSVD() */ - BDCSVD(Index rows, Index cols, unsigned int computationOptions = 0) - : m_algoswap(16), m_numIters(0) - { + BDCSVD(Index rows, Index cols, unsigned int computationOptions) : m_algoswap(16), m_numIters(0) { + internal::check_svd_constructor_assertions(computationOptions); allocate(rows, cols, computationOptions); } - /** \brief Constructor performing the decomposition of given matrix. + /** \brief Constructor performing the decomposition of given matrix, using the custom options specified + * with the \a Options template paramter. * * \param matrix the matrix to decompose - * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. - * By default, none is computed. This is a bit - field, the possible bits are #ComputeFullU, #ComputeThinU, - * #ComputeFullV, #ComputeThinV. - * - * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not - * available with the (non - default) FullPivHouseholderQR preconditioner. */ - BDCSVD(const MatrixType& matrix, unsigned int computationOptions = 0) - : m_algoswap(16), m_numIters(0) - { - compute(matrix, computationOptions); + BDCSVD(const MatrixType& matrix) : m_algoswap(16), m_numIters(0) { + compute_impl(matrix, internal::get_computation_options(Options)); } - ~BDCSVD() - { - } - - /** \brief Method performing the decomposition of given matrix using custom options. + /** \brief Constructor performing the decomposition of given matrix using specified options + * for computing unitaries. + * + * Note: This constructor is deprecated. + * One \b cannot request unitiaries using both the \a Options template parameter + * and the constructor. If possible, prefer using the \a Options template parameter. * * \param matrix the matrix to decompose - * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. - * By default, none is computed. This is a bit - field, the possible bits are #ComputeFullU, #ComputeThinU, - * #ComputeFullV, #ComputeThinV. - * - * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not - * available with the (non - default) FullPivHouseholderQR preconditioner. + * \param computationOptions specifification for computing Thin/Full unitaries U/V */ - BDCSVD& compute(const MatrixType& matrix, unsigned int computationOptions); + BDCSVD(const MatrixType& matrix, unsigned int computationOptions) : m_algoswap(16), m_numIters(0) { + internal::check_svd_constructor_assertions(computationOptions); + compute_impl(matrix, computationOptions); + } - /** \brief Method performing the decomposition of given matrix using current options. + ~BDCSVD() {} + + /** \brief Method performing the decomposition of given matrix. Computes Thin/Full unitaries U/V if specified + * using the \a Options template parameter or the class constructor. * * \param matrix the matrix to decompose - * - * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int). */ - BDCSVD& compute(const MatrixType& matrix) - { - return compute(matrix, this->m_computationOptions); - } + BDCSVD& compute(const MatrixType& matrix) { return compute_impl(matrix, m_computationOptions); } void setSwitchSize(int s) { @@ -184,6 +205,7 @@ public: private: void allocate(Index rows, Index cols, unsigned int computationOptions); + BDCSVD& compute_impl(const MatrixType& matrix, unsigned int computationOptions); void divide(Index firstCol, Index lastCol, Index firstRowW, Index firstColW, Index shift); void computeSVDofM(Index firstCol, Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V); void computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm, VectorType& singVals, ArrayRef shifts, ArrayRef mus); @@ -196,8 +218,10 @@ private: void copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naivev); void structured_update(Block A, const MatrixXr &B, Index n1); static RealScalar secularEq(RealScalar x, const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const ArrayRef& diagShifted, RealScalar shift); + template + void computeBaseCase(SVDType& svd, Index n, Index firstCol, Index firstRowW, Index firstColW, Index shift); -protected: + protected: MatrixXr m_naiveU, m_naiveV; MatrixXr m_computed; Index m_nRec; @@ -205,36 +229,36 @@ protected: ArrayXi m_workspaceI; int m_algoswap; bool m_isTranspose, m_compU, m_compV; + JacobiSVD smallSvd; - using Base::m_singularValues; - using Base::m_diagSize; - using Base::m_computeFullU; - using Base::m_computeFullV; + using Base::m_computationOptions; using Base::m_computeThinU; using Base::m_computeThinV; - using Base::m_matrixU; - using Base::m_matrixV; + using Base::m_diagSize; using Base::m_info; using Base::m_isInitialized; + using Base::m_matrixU; + using Base::m_matrixV; using Base::m_nonzeroSingularValues; + using Base::m_singularValues; -public: + public: int m_numIters; -}; //end class BDCSVD - +}; // end class BDCSVD // Method to allocate and initialize matrix and attributes -template -void BDCSVD::allocate(Eigen::Index rows, Eigen::Index cols, unsigned int computationOptions) -{ - m_isTranspose = (cols > rows); - +template +void BDCSVD::allocate(Eigen::Index rows, Eigen::Index cols, unsigned int computationOptions) { if (Base::allocate(rows, cols, computationOptions)) return; + if (cols < m_algoswap) + internal::allocate_small_svd::run(smallSvd, rows, cols, computationOptions); + m_computed = MatrixXr::Zero(m_diagSize + 1, m_diagSize ); m_compU = computeV(); m_compV = computeU(); + m_isTranspose = (cols > rows); if (m_isTranspose) std::swap(m_compU, m_compV); @@ -245,31 +269,31 @@ void BDCSVD::allocate(Eigen::Index rows, Eigen::Index cols, unsigned m_workspace.resize((m_diagSize+1)*(m_diagSize+1)*3); m_workspaceI.resize(3*m_diagSize); -}// end allocate +} // end allocate -template -BDCSVD& BDCSVD::compute(const MatrixType& matrix, unsigned int computationOptions) -{ +template +BDCSVD& BDCSVD::compute_impl(const MatrixType& matrix, + unsigned int computationOptions) { #ifdef EIGEN_BDCSVD_DEBUG_VERBOSE std::cout << "\n\n\n======================================================================================================================\n\n\n"; #endif - allocate(matrix.rows(), matrix.cols(), computationOptions); using std::abs; + allocate(matrix.rows(), matrix.cols(), computationOptions); + const RealScalar considerZero = (std::numeric_limits::min)(); //**** step -1 - If the problem is too small, directly falls back to JacobiSVD and return if(matrix.cols() < m_algoswap) { - // FIXME this line involves temporaries - JacobiSVD jsvd(matrix,computationOptions); + smallSvd.compute(matrix); m_isInitialized = true; - m_info = jsvd.info(); + m_info = smallSvd.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(); + if (computeU()) m_matrixU = smallSvd.matrixU(); + if (computeV()) m_matrixV = smallSvd.matrixV(); + m_singularValues = smallSvd.singularValues(); + m_nonzeroSingularValues = smallSvd.nonzeroSingularValues(); } return *this; } @@ -330,13 +354,12 @@ BDCSVD& BDCSVD::compute(const MatrixType& matrix, unsign m_isInitialized = true; return *this; -}// end compute +} // end compute - -template -template -void BDCSVD::copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naiveV) -{ +template +template +void BDCSVD::copyUV(const HouseholderU& householderU, const HouseholderV& householderV, + const NaiveU& naiveU, const NaiveV& naiveV) { // Note exchange of U and V: m_matrixU is set from m_naiveV and vice versa if (computeU()) { @@ -362,9 +385,8 @@ void BDCSVD::copyUV(const HouseholderU &householderU, const Househol * We can thus pack them prior to the the matrix product. However, this is only worth the effort if the matrix is large * enough. */ -template -void BDCSVD::structured_update(Block A, const MatrixXr &B, Index n1) -{ +template +void BDCSVD::structured_update(Block A, const MatrixXr& B, Index n1) { Index n = A.rows(); if(n>100) { @@ -403,8 +425,26 @@ void BDCSVD::structured_update(Block A, co } } -// The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods takes as argument the -// place of the submatrix we are currently working on. +template +template +void BDCSVD::computeBaseCase(SVDType& svd, Index n, Index firstCol, Index firstRowW, + Index firstColW, Index shift) { + svd.compute(m_computed.block(firstCol, firstCol, n + 1, n)); + m_info = svd.info(); + if (m_info != Success && m_info != NoConvergence) return; + if (m_compU) + m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = svd.matrixU(); + else { + m_naiveU.row(0).segment(firstCol, n + 1).real() = svd.matrixU().row(0); + m_naiveU.row(1).segment(firstCol, n + 1).real() = svd.matrixU().row(n); + } + if (m_compV) m_naiveV.block(firstRowW, firstColW, n, n).real() = svd.matrixV(); + m_computed.block(firstCol + shift, firstCol + shift, n + 1, n).setZero(); + m_computed.diagonal().segment(firstCol + shift, n) = svd.singularValues().head(n); +} + +// The divide algorithm is done "in place", we are always working on subsets of the same matrix. The divide methods +// takes as argument the place of the submatrix we are currently working on. //@param firstCol : The Index of the first column of the submatrix of m_computed and for m_naiveU; //@param lastCol : The Index of the last column of the submatrix of m_computed and for m_naiveU; @@ -413,9 +453,9 @@ void BDCSVD::structured_update(Block A, co //@param firstColW : Same as firstRowW with the column. //@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 -void BDCSVD::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift) -{ +template +void BDCSVD::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, + Eigen::Index firstColW, Eigen::Index shift) { // requires rows = cols + 1; using std::pow; using std::sqrt; @@ -432,20 +472,14 @@ void BDCSVD::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig // matrices. if (n < m_algoswap) { - // FIXME this line involves temporaries - JacobiSVD 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 - { - m_naiveU.row(0).segment(firstCol, n + 1).real() = b.matrixU().row(0); - m_naiveU.row(1).segment(firstCol, n + 1).real() = b.matrixU().row(n); + // FIXME this block involves temporaries + if (m_compV) { + JacobiSVD baseSvd; + computeBaseCase(baseSvd, n, firstCol, firstRowW, firstColW, shift); + } else { + JacobiSVD baseSvd; + computeBaseCase(baseSvd, n, firstCol, firstRowW, firstColW, shift); } - if (m_compV) m_naiveV.block(firstRowW, firstColW, n, n).real() = b.matrixV(); - m_computed.block(firstCol + shift, firstCol + shift, n + 1, n).setZero(); - m_computed.diagonal().segment(firstCol + shift, n) = b.singularValues().head(n); return; } // We use the divide and conquer algorithm @@ -587,7 +621,7 @@ void BDCSVD::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig m_computed.block(firstCol + shift, firstCol + shift, n, n).setZero(); m_computed.block(firstCol + shift, firstCol + shift, n, n).diagonal() = singVals; -}// end divide +} // end divide // Compute SVD of m_computed.block(firstCol, firstCol, n + 1, n); this block only has non-zeros in // the first column and on the diagonal and has undergone deflation, so diagonal is in increasing @@ -597,9 +631,9 @@ void BDCSVD::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig // TODO Opportunities for optimization: better root finding algo, better stopping criterion, better // handling of round-off errors, be consistent in ordering // For instance, to solve the secular equation using FMM, see http://www.stat.uchicago.edu/~lekheng/courses/302/classics/greengard-rokhlin.pdf -template -void BDCSVD::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V) -{ +template +void BDCSVD::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, MatrixXr& U, + VectorType& singVals, MatrixXr& V) { const RealScalar considerZero = (std::numeric_limits::min)(); using std::abs; ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n); @@ -728,9 +762,10 @@ void BDCSVD::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma #endif } -template -typename BDCSVD::RealScalar BDCSVD::secularEq(RealScalar mu, const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const ArrayRef& diagShifted, RealScalar shift) -{ +template +typename BDCSVD::RealScalar BDCSVD::secularEq( + RealScalar mu, const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm, const ArrayRef& diagShifted, + RealScalar shift) { Index m = perm.size(); RealScalar res = Literal(1); for(Index i=0; i::RealScalar BDCSVD::secularEq(RealScalar res += (col0(j) / (diagShifted(j) - mu)) * (col0(j) / (diag(j) + shift + mu)); } return res; - } -template -void BDCSVD::computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, - VectorType& singVals, ArrayRef shifts, ArrayRef mus) -{ +template +void BDCSVD::computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm, + VectorType& singVals, ArrayRef shifts, ArrayRef mus) { using std::abs; using std::swap; using std::sqrt; @@ -993,13 +1026,11 @@ void BDCSVD::computeSingVals(const ArrayRef& col0, const ArrayRef& d } } - // zhat is perturbation of col0 for which singular vectors can be computed stably (see Section 3.1) -template -void BDCSVD::perturbCol0 - (const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const VectorType& singVals, - const ArrayRef& shifts, const ArrayRef& mus, ArrayRef zhat) -{ +template +void BDCSVD::perturbCol0(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef& perm, + const VectorType& singVals, const ArrayRef& shifts, const ArrayRef& mus, + ArrayRef zhat) { using std::sqrt; Index n = col0.size(); Index m = perm.size(); @@ -1075,11 +1106,10 @@ void BDCSVD::perturbCol0 } // compute singular vectors -template -void BDCSVD::computeSingVecs - (const ArrayRef& zhat, const ArrayRef& diag, const IndicesRef &perm, const VectorType& singVals, - const ArrayRef& shifts, const ArrayRef& mus, MatrixXr& U, MatrixXr& V) -{ +template +void BDCSVD::computeSingVecs(const ArrayRef& zhat, const ArrayRef& diag, const IndicesRef& perm, + const VectorType& singVals, const ArrayRef& shifts, + const ArrayRef& mus, MatrixXr& U, MatrixXr& V) { Index n = zhat.size(); Index m = perm.size(); @@ -1117,13 +1147,12 @@ void BDCSVD::computeSingVecs U.col(n) = VectorType::Unit(n+1, n); } - // page 12_13 // i >= 1, di almost null and zi non null. // We use a rotation to zero out zi applied to the left of M -template -void BDCSVD::deflation43(Eigen::Index firstCol, Eigen::Index shift, Eigen::Index i, Eigen::Index size) -{ +template +void BDCSVD::deflation43(Eigen::Index firstCol, Eigen::Index shift, Eigen::Index i, + Eigen::Index size) { using std::abs; using std::sqrt; using std::pow; @@ -1143,16 +1172,16 @@ void BDCSVD::deflation43(Eigen::Index firstCol, Eigen::Index shift, JacobiRotation J(c/r,-s/r); if (m_compU) m_naiveU.middleRows(firstCol, size+1).applyOnTheRight(firstCol, firstCol+i, J); else m_naiveU.applyOnTheRight(firstCol, firstCol+i, J); -}// end deflation 43 - +} // end deflation 43 // page 13 // i,j >= 1, i!=j and |di - dj| < epsilon * norm2(M) // We apply two rotations to have zj = 0; // TODO deflation44 is still broken and not properly tested -template -void BDCSVD::deflation44(Eigen::Index firstColu , Eigen::Index firstColm, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index i, Eigen::Index j, Eigen::Index size) -{ +template +void BDCSVD::deflation44(Eigen::Index firstColu, Eigen::Index firstColm, Eigen::Index firstRowW, + Eigen::Index firstColW, Eigen::Index i, Eigen::Index j, + Eigen::Index size) { using std::abs; using std::sqrt; using std::conj; @@ -1186,13 +1215,12 @@ void BDCSVD::deflation44(Eigen::Index firstColu , Eigen::Index first if (m_compU) m_naiveU.middleRows(firstColu, size+1).applyOnTheRight(firstColu + i, firstColu + j, J); else m_naiveU.applyOnTheRight(firstColu+i, firstColu+j, J); if (m_compV) m_naiveV.middleRows(firstRowW, size).applyOnTheRight(firstColW + i, firstColW + j, J); -}// end deflation 44 - +} // end deflation 44 // acts on block from (firstCol+shift, firstCol+shift) to (lastCol+shift, lastCol+shift) [inclusive] -template -void BDCSVD::deflation(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index k, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift) -{ +template +void BDCSVD::deflation(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index k, + Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift) { using std::sqrt; using std::abs; const Index length = lastCol + 1 - firstCol; @@ -1360,7 +1388,7 @@ void BDCSVD::deflation(Eigen::Index firstCol, Eigen::Index lastCol, assert(m_naiveV.allFinite()); assert(m_computed.allFinite()); #endif -}//end deflation +} // end deflation /** \svd_module * @@ -1368,11 +1396,23 @@ void BDCSVD::deflation(Eigen::Index firstCol, Eigen::Index lastCol, * * \sa class BDCSVD */ -template -BDCSVD::PlainObject> -MatrixBase::bdcSvd(unsigned int computationOptions) const -{ - return BDCSVD(*this, computationOptions); +template +template +BDCSVD::PlainObject, Options> MatrixBase::bdcSvd() const { + return BDCSVD(*this); +} + +/** \svd_module + * + * \return the singular value decomposition of \c *this computed by Divide & Conquer algorithm + * + * \sa class BDCSVD + */ +template +template +BDCSVD::PlainObject, Options> MatrixBase::bdcSvd( + unsigned int computationOptions) const { + return BDCSVD(*this, computationOptions); } } // end namespace Eigen diff --git a/Eigen/src/SVD/JacobiSVD.h b/Eigen/src/SVD/JacobiSVD.h index c57b20199..20a276cdc 100644 --- a/Eigen/src/SVD/JacobiSVD.h +++ b/Eigen/src/SVD/JacobiSVD.h @@ -13,13 +13,13 @@ #include "./InternalHeaderCheck.h" -namespace Eigen { +namespace Eigen { namespace internal { + // forward declaration (needed by ICC) // the empty body is required by MSVC -template::IsComplex> +template ::IsComplex> struct svd_precondition_2x2_block_to_be_real {}; /*** QR preconditioners (R-SVD) @@ -46,37 +46,31 @@ struct qr_preconditioner_should_do_anything }; }; -template::ret -> struct qr_preconditioner_impl {}; +template ::ret> +struct qr_preconditioner_impl {}; -template -class qr_preconditioner_impl -{ -public: - void allocate(const JacobiSVD&) {} - bool run(JacobiSVD&, const MatrixType&) - { - return false; - } +template +class qr_preconditioner_impl { + public: + void allocate(const JacobiSVD&) {} + bool run(JacobiSVD&, const MatrixType&) { return false; } }; /*** preconditioner using FullPivHouseholderQR ***/ -template -class qr_preconditioner_impl -{ -public: +template +class qr_preconditioner_impl { + public: typedef typename MatrixType::Scalar Scalar; - enum - { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime - }; - typedef Matrix WorkspaceType; + typedef JacobiSVD SVDType; - void allocate(const JacobiSVD& svd) - { + enum { WorkspaceSize = MatrixType::RowsAtCompileTime, MaxWorkspaceSize = MatrixType::MaxRowsAtCompileTime }; + + typedef Matrix WorkspaceType; + + void allocate(const SVDType& svd) { if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { m_qr.~QRType(); @@ -85,8 +79,7 @@ public: if (svd.m_computeFullU) m_workspace.resize(svd.rows()); } - bool run(JacobiSVD& svd, const MatrixType& matrix) - { + bool run(SVDType& svd, const MatrixType& matrix) { if(matrix.rows() > matrix.cols()) { m_qr.compute(matrix); @@ -97,32 +90,33 @@ public: } return false; } + private: typedef FullPivHouseholderQR QRType; QRType m_qr; WorkspaceType m_workspace; }; -template -class qr_preconditioner_impl -{ -public: +template +class qr_preconditioner_impl { + public: typedef typename MatrixType::Scalar Scalar; - enum - { + typedef JacobiSVD SVDType; + + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - Options = MatrixType::Options + MatrixOptions = MatrixType::Options }; - typedef typename internal::make_proper_matrix_type< - Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime - >::type TransposeTypeWithSameStorageOrder; + typedef typename internal::make_proper_matrix_type::type + TransposeTypeWithSameStorageOrder; - void allocate(const JacobiSVD& svd) - { + void allocate(const SVDType& svd) { if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { m_qr.~QRType(); @@ -132,8 +126,7 @@ public: if (svd.m_computeFullV) m_workspace.resize(svd.cols()); } - bool run(JacobiSVD& svd, const MatrixType& matrix) - { + bool run(SVDType& svd, const MatrixType& matrix) { if(matrix.cols() > matrix.rows()) { m_adjoint = matrix.adjoint(); @@ -145,21 +138,31 @@ public: } else return false; } + private: typedef FullPivHouseholderQR QRType; QRType m_qr; TransposeTypeWithSameStorageOrder m_adjoint; - typename internal::plain_row_type::type m_workspace; + typename plain_row_type::type m_workspace; }; /*** preconditioner using ColPivHouseholderQR ***/ -template -class qr_preconditioner_impl -{ -public: - void allocate(const JacobiSVD& svd) - { +template +class qr_preconditioner_impl { + public: + typedef typename MatrixType::Scalar Scalar; + typedef JacobiSVD SVDType; + + enum { + WorkspaceSize = internal::traits::MatrixUColsAtCompileTime, + MaxWorkspaceSize = internal::traits::MatrixUMaxColsAtCompileTime + }; + + typedef Matrix WorkspaceType; + + void allocate(const SVDType& svd) { if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { m_qr.~QRType(); @@ -169,8 +172,7 @@ public: else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); } - bool run(JacobiSVD& svd, const MatrixType& matrix) - { + bool run(SVDType& svd, const MatrixType& matrix) { if(matrix.rows() > matrix.cols()) { m_qr.compute(matrix); @@ -190,29 +192,33 @@ public: private: typedef ColPivHouseholderQR QRType; QRType m_qr; - typename internal::plain_col_type::type m_workspace; + WorkspaceType m_workspace; }; -template -class qr_preconditioner_impl -{ -public: +template +class qr_preconditioner_impl { + public: typedef typename MatrixType::Scalar Scalar; - enum - { + typedef JacobiSVD SVDType; + + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - Options = MatrixType::Options + MatrixOptions = MatrixType::Options, + WorkspaceSize = internal::traits::MatrixVColsAtCompileTime, + MaxWorkspaceSize = internal::traits::MatrixVMaxColsAtCompileTime }; - typedef typename internal::make_proper_matrix_type< - Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime - >::type TransposeTypeWithSameStorageOrder; + typedef Matrix WorkspaceType; - void allocate(const JacobiSVD& svd) - { + typedef typename internal::make_proper_matrix_type::type + TransposeTypeWithSameStorageOrder; + + void allocate(const SVDType& svd) { if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { m_qr.~QRType(); @@ -223,8 +229,7 @@ public: m_adjoint.resize(svd.cols(), svd.rows()); } - bool run(JacobiSVD& svd, const MatrixType& matrix) - { + bool run(SVDType& svd, const MatrixType& matrix) { if(matrix.cols() > matrix.rows()) { m_adjoint = matrix.adjoint(); @@ -247,17 +252,25 @@ private: typedef ColPivHouseholderQR QRType; QRType m_qr; TransposeTypeWithSameStorageOrder m_adjoint; - typename internal::plain_row_type::type m_workspace; + WorkspaceType m_workspace; }; /*** preconditioner using HouseholderQR ***/ -template -class qr_preconditioner_impl -{ -public: - void allocate(const JacobiSVD& svd) - { +template +class qr_preconditioner_impl { + public: + typedef typename MatrixType::Scalar Scalar; + typedef JacobiSVD SVDType; + + enum { + WorkspaceSize = internal::traits::MatrixUColsAtCompileTime, + MaxWorkspaceSize = internal::traits::MatrixUMaxColsAtCompileTime + }; + + typedef Matrix WorkspaceType; + + void allocate(const SVDType& svd) { if (svd.rows() != m_qr.rows() || svd.cols() != m_qr.cols()) { m_qr.~QRType(); @@ -267,8 +280,7 @@ public: else if (svd.m_computeThinU) m_workspace.resize(svd.cols()); } - bool run(JacobiSVD& svd, const MatrixType& matrix) - { + bool run(SVDType& svd, const MatrixType& matrix) { if(matrix.rows() > matrix.cols()) { m_qr.compute(matrix); @@ -284,32 +296,36 @@ public: } return false; } + private: typedef HouseholderQR QRType; QRType m_qr; - typename internal::plain_col_type::type m_workspace; + WorkspaceType m_workspace; }; -template -class qr_preconditioner_impl -{ -public: +template +class qr_preconditioner_impl { + public: typedef typename MatrixType::Scalar Scalar; - enum - { + typedef JacobiSVD SVDType; + + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - Options = MatrixType::Options + MatrixOptions = MatrixType::Options, + WorkspaceSize = internal::traits::MatrixVColsAtCompileTime, + MaxWorkspaceSize = internal::traits::MatrixVMaxColsAtCompileTime }; - typedef typename internal::make_proper_matrix_type< - Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime - >::type TransposeTypeWithSameStorageOrder; + typedef Matrix WorkspaceType; - void allocate(const JacobiSVD& svd) - { + typedef typename internal::make_proper_matrix_type::type + TransposeTypeWithSameStorageOrder; + + void allocate(const SVDType& svd) { if (svd.cols() != m_qr.rows() || svd.rows() != m_qr.cols()) { m_qr.~QRType(); @@ -320,8 +336,7 @@ public: m_adjoint.resize(svd.cols(), svd.rows()); } - bool run(JacobiSVD& svd, const MatrixType& matrix) - { + bool run(SVDType& svd, const MatrixType& matrix) { if(matrix.cols() > matrix.rows()) { m_adjoint = matrix.adjoint(); @@ -344,7 +359,7 @@ private: typedef HouseholderQR QRType; QRType m_qr; TransposeTypeWithSameStorageOrder m_adjoint; - typename internal::plain_row_type::type m_workspace; + WorkspaceType m_workspace; }; /*** 2x2 SVD implementation @@ -352,18 +367,16 @@ private: *** JacobiSVD consists in performing a series of 2x2 SVD subproblems ***/ -template -struct svd_precondition_2x2_block_to_be_real -{ - typedef JacobiSVD SVD; +template +struct svd_precondition_2x2_block_to_be_real { + typedef JacobiSVD SVD; typedef typename MatrixType::RealScalar RealScalar; static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; } }; -template -struct svd_precondition_2x2_block_to_be_real -{ - typedef JacobiSVD SVD; +template +struct svd_precondition_2x2_block_to_be_real { + typedef JacobiSVD SVD; typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; static bool run(typename SVD::WorkMatrixType& work_matrix, SVD& svd, Index p, Index q, RealScalar& maxDiagEntry) @@ -425,249 +438,235 @@ struct svd_precondition_2x2_block_to_be_real } }; -template -struct traits > - : traits -{ +template +struct traits > : svd_traits { typedef MatrixType_ MatrixType; }; } // end namespace internal /** \ingroup SVD_Module - * - * - * \class JacobiSVD - * - * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix - * - * \tparam MatrixType_ the type of the matrix of which we are computing the SVD decomposition - * \tparam QRPreconditioner this optional parameter allows to specify the type of QR decomposition that will be used internally - * for the R-SVD step for non-square matrices. See discussion of possible values below. - * - * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product - * \f[ A = U S V^* \f] - * where \a U is a n-by-n unitary, \a V is a p-by-p unitary, and \a S is a n-by-p real positive matrix which is zero outside of its main diagonal; - * the diagonal entries of S are known as the \em singular \em values of \a A and the columns of \a U and \a V are known as the left - * and right \em singular \em vectors of \a A respectively. - * - * Singular values are always sorted in decreasing order. - * - * This JacobiSVD decomposition computes only the singular values by default. If you want \a U or \a V, you need to ask for them explicitly. - * - * You can ask for only \em thin \a U or \a V to be computed, meaning the following. In case of a rectangular n-by-p matrix, letting \a m be the - * smaller value among \a n and \a p, there are only \a m singular vectors; the remaining columns of \a U and \a V do not correspond to actual - * 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. - * - * Here's an example demonstrating basic usage: - * \include JacobiSVD_basic.cpp - * Output: \verbinclude JacobiSVD_basic.out - * - * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The downside is that it's slower than - * bidiagonalizing SVD algorithms for large square matrices; however its complexity is still \f$ O(n^2p) \f$ where \a n is the smaller dimension and - * \a p is the greater dimension, meaning that it is still of the same order of complexity as the faster bidiagonalizing R-SVD algorithms. - * In particular, like any R-SVD, it takes advantage of non-squareness in that its complexity is only linear in the greater dimension. - * - * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is guaranteed to - * terminate in finite (and reasonable) time. - * - * The possible values for QRPreconditioner are: - * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR. - * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR. - * Contrary to other QRs, it doesn't allow computing thin unitaries. - * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses non-pivoting QR. - * This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing SVD algorithms (since bidiagonalization - * is inherently non-pivoting). However the resulting SVD is still more reliable than bidiagonalizing SVDs because the Jacobi-based iterarive - * process is more reliable than the optimized bidiagonal SVD iterations. - * \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that you will only be computing - * JacobiSVD decompositions of square matrices. Non-square matrices require a QR preconditioner. Using this option will result in - * faster compilation and smaller executable code. It won't significantly speed up computation, since JacobiSVD is always checking - * if QR preconditioning is needed before applying it anyway. - * - * \sa MatrixBase::jacobiSvd() - */ -template class JacobiSVD - : public SVDBase > -{ - typedef SVDBase Base; - public: + * + * + * \class JacobiSVD + * + * \brief Two-sided Jacobi SVD decomposition of a rectangular matrix + * + * \tparam MatrixType_ the type of the matrix of which we are computing the SVD decomposition + * \tparam Options this optional parameter allows one to specify the type of QR decomposition that will be used + * internally for the R-SVD step for non-square matrices. Additionally, it allows one to specify whether to compute thin + * or full unitaries \a U and \a V. See discussion of possible values below. + * + * SVD decomposition consists in decomposing any n-by-p matrix \a A as a product + * \f[ A = U S V^* \f] + * where \a U is a n-by-n unitary, \a V is a p-by-p unitary, and \a S is a n-by-p real positive matrix which is zero + * outside of its main diagonal; the diagonal entries of S are known as the \em singular \em values of \a A and the + * columns of \a U and \a V are known as the left and right \em singular \em vectors of \a A respectively. + * + * Singular values are always sorted in decreasing order. + * + * This JacobiSVD decomposition computes only the singular values by default. If you want \a U or \a V, you need to ask + * for them explicitly. + * + * You can ask for only \em thin \a U or \a V to be computed, meaning the following. In case of a rectangular n-by-p + * matrix, letting \a m be the smaller value among \a n and \a p, there are only \a m singular vectors; the remaining + * columns of \a U and \a V do not correspond to actual 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. + * + * Here's an example demonstrating basic usage: + * \include JacobiSVD_basic.cpp + * Output: \verbinclude JacobiSVD_basic.out + * + * This JacobiSVD class is a two-sided Jacobi R-SVD decomposition, ensuring optimal reliability and accuracy. The + * downside is that it's slower than bidiagonalizing SVD algorithms for large square matrices; however its complexity is + * still \f$ O(n^2p) \f$ where \a n is the smaller dimension and \a p is the greater dimension, meaning that it is still + * of the same order of complexity as the faster bidiagonalizing R-SVD algorithms. In particular, like any R-SVD, it + * takes advantage of non-squareness in that its complexity is only linear in the greater dimension. + * + * If the input matrix has inf or nan coefficients, the result of the computation is undefined, but the computation is + * guaranteed to terminate in finite (and reasonable) time. + * + * The possible QR preconditioners that can be set with Options template parameter are: + * \li ColPivHouseholderQRPreconditioner is the default. In practice it's very safe. It uses column-pivoting QR. + * \li FullPivHouseholderQRPreconditioner, is the safest and slowest. It uses full-pivoting QR. + * Contrary to other QRs, it doesn't allow computing thin unitaries. + * \li HouseholderQRPreconditioner is the fastest, and less safe and accurate than the pivoting variants. It uses + * non-pivoting QR. This is very similar in safety and accuracy to the bidiagonalization process used by bidiagonalizing + * SVD algorithms (since bidiagonalization is inherently non-pivoting). However the resulting SVD is still more reliable + * than bidiagonalizing SVDs because the Jacobi-based iterarive process is more reliable than the optimized bidiagonal + * SVD iterations. \li NoQRPreconditioner allows not to use a QR preconditioner at all. This is useful if you know that + * you will only be computing JacobiSVD decompositions of square matrices. Non-square matrices require a QR + * preconditioner. Using this option will result in faster compilation and smaller executable code. It won't + * significantly speed up computation, since JacobiSVD is always checking if QR preconditioning is needed before + * applying it anyway. + * + * One may also use the Options template parameter to specify how the unitaries should be computed. The options are + * #ComputeThinU, #ComputeThinV, #ComputeFullU, #ComputeFullV. It is not possible to request both the thin and full + * versions of a unitary. By default, unitaries will not be computed. + * + * You can set the QRPreconditioner and unitary options together: JacobiSVD + * + * \sa MatrixBase::jacobiSvd() + */ +template +class JacobiSVD : public SVDBase > { + typedef SVDBase Base; - typedef MatrixType_ MatrixType; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits::Real RealScalar; - enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime, - DiagSizeAtCompileTime = internal::min_size_prefer_dynamic(RowsAtCompileTime,ColsAtCompileTime), - MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - MaxDiagSizeAtCompileTime = internal::min_size_prefer_fixed(MaxRowsAtCompileTime,MaxColsAtCompileTime), - MatrixOptions = MatrixType::Options - }; + public: + typedef MatrixType_ MatrixType; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + enum { + QRPreconditioner = internal::get_qr_preconditioner(Options), + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + DiagSizeAtCompileTime = internal::min_size_prefer_dynamic(RowsAtCompileTime, ColsAtCompileTime), + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + MaxDiagSizeAtCompileTime = internal::min_size_prefer_fixed(MaxRowsAtCompileTime, MaxColsAtCompileTime), + MatrixOptions = MatrixType::Options + }; - typedef typename Base::MatrixUType MatrixUType; - typedef typename Base::MatrixVType MatrixVType; - typedef typename Base::SingularValuesType SingularValuesType; - - typedef typename internal::plain_row_type::type RowType; - typedef typename internal::plain_col_type::type ColType; - typedef Matrix - WorkMatrixType; + typedef typename Base::MatrixUType MatrixUType; + typedef typename Base::MatrixVType MatrixVType; + typedef typename Base::SingularValuesType SingularValuesType; + typedef Matrix + WorkMatrixType; - /** \brief Default Constructor. - * - * The default constructor is useful in cases in which the user intends to - * perform decompositions via JacobiSVD::compute(const MatrixType&). - */ - JacobiSVD() - {} + /** \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via JacobiSVD::compute(const MatrixType&). + */ + JacobiSVD() {} + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem size and \a Options template parameter. + * + * \sa JacobiSVD() + */ + JacobiSVD(Index rows, Index cols) { allocate(rows, cols, internal::get_computation_options(Options)); } - /** \brief Default Constructor with memory preallocation - * - * Like the default constructor but with preallocation of the internal data - * according to the specified problem size. - * \sa JacobiSVD() - */ - JacobiSVD(Index rows, Index cols, unsigned int computationOptions = 0) - { - allocate(rows, cols, computationOptions); - } + /** \brief Default Constructor with memory preallocation + * + * Like the default constructor but with preallocation of the internal data + * according to the specified problem size. + * + * Note: This constructor is deprecated. + * one \b cannot request unitaries using both the \a Options template parameter + * and the constructor. If possible, prefer using the \a Options template parameter. + * + * \param computationOptions specify whether to compute Thin/Full unitaries U/V + * \sa JacobiSVD() + */ + JacobiSVD(Index rows, Index cols, unsigned int computationOptions) { + internal::check_svd_constructor_assertions(computationOptions); + allocate(rows, cols, computationOptions); + } - /** \brief Constructor performing the decomposition of given matrix. - * - * \param matrix the matrix to decompose - * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. - * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, - * #ComputeFullV, #ComputeThinV. - * - * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not - * available with the (non-default) FullPivHouseholderQR preconditioner. - */ - explicit JacobiSVD(const MatrixType& matrix, unsigned int computationOptions = 0) - { - compute(matrix, computationOptions); - } + /** \brief Constructor performing the decomposition of given matrix, using the custom options specified + * with the \a Options template paramter. + * + * \param matrix the matrix to decompose + */ + explicit JacobiSVD(const MatrixType& matrix) { compute_impl(matrix, internal::get_computation_options(Options)); } - /** \brief Method performing the decomposition of given matrix using custom options. - * - * \param matrix the matrix to decompose - * \param computationOptions optional parameter allowing to specify if you want full or thin U or V unitaries to be computed. - * By default, none is computed. This is a bit-field, the possible bits are #ComputeFullU, #ComputeThinU, - * #ComputeFullV, #ComputeThinV. - * - * Thin unitaries are only available if your matrix type has a Dynamic number of columns (for example MatrixXf). They also are not - * available with the (non-default) FullPivHouseholderQR preconditioner. - */ - JacobiSVD& compute(const MatrixType& matrix, unsigned int computationOptions); + /** \brief Constructor performing the decomposition of given matrix using specified options + * for computing unitaries. + * + * Note: This constructor is deprecated. + * One \b cannot request unitiaries using both the \a Options template parameter + * and the constructor. If possible, prefer using the \a Options template parameter. + * + * \param matrix the matrix to decompose + * \param computationOptions specify whether to compute Thin/Full unitaries U/V + */ + JacobiSVD(const MatrixType& matrix, unsigned int computationOptions) { + internal::check_svd_constructor_assertions(computationOptions); + compute_impl(matrix, computationOptions); + } - /** \brief Method performing the decomposition of given matrix using current options. - * - * \param matrix the matrix to decompose - * - * This method uses the current \a computationOptions, as already passed to the constructor or to compute(const MatrixType&, unsigned int). - */ - JacobiSVD& compute(const MatrixType& matrix) - { - return compute(matrix, m_computationOptions); - } + /** \brief Method performing the decomposition of given matrix. Computes Thin/Full unitaries U/V if specified + * using the \a Options template parameter or the class constructor. + * + * \param matrix the matrix to decompose + */ + JacobiSVD& compute(const MatrixType& matrix) { return compute_impl(matrix, m_computationOptions); } - using Base::computeU; - using Base::computeV; - using Base::rows; - using Base::cols; - using Base::rank; + using Base::computeU; + using Base::computeV; + using Base::rows; + using Base::cols; + using Base::rank; - private: - void allocate(Index rows, Index cols, unsigned int computationOptions); + private: + void allocate(Index rows, Index cols, unsigned int computationOptions); + JacobiSVD& compute_impl(const MatrixType& matrix, unsigned int computationOptions); - protected: - 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; - using Base::m_computeFullU; - using Base::m_computeThinU; - using Base::m_computeFullV; - using Base::m_computeThinV; - using Base::m_computationOptions; - using Base::m_nonzeroSingularValues; - using Base::m_rows; - using Base::m_cols; - using Base::m_diagSize; - using Base::m_prescribedThreshold; - WorkMatrixType m_workMatrix; + 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; + using Base::m_matrixU; + 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; + using Base::ShouldComputeThinV; - template - friend struct internal::svd_precondition_2x2_block_to_be_real; - template - friend struct internal::qr_preconditioner_impl; + EIGEN_STATIC_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.") - internal::qr_preconditioner_impl m_qr_precond_morecols; - internal::qr_preconditioner_impl m_qr_precond_morerows; - MatrixType m_scaledMatrix; + template + friend struct internal::svd_precondition_2x2_block_to_be_real; + template + friend struct internal::qr_preconditioner_impl; + + internal::qr_preconditioner_impl + m_qr_precond_morecols; + internal::qr_preconditioner_impl + m_qr_precond_morerows; + WorkMatrixType m_workMatrix; + MatrixType m_scaledMatrix; }; -template -void JacobiSVD::allocate(Eigen::Index rows, Eigen::Index cols, unsigned int computationOptions) -{ - eigen_assert(rows >= 0 && cols >= 0); +template +void JacobiSVD::allocate(Eigen::Index rows, Eigen::Index cols, unsigned int computationOptions) { + if (Base::allocate(rows, cols, computationOptions)) return; - if (m_isAllocated && - rows == m_rows && - cols == m_cols && - computationOptions == m_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_rows = rows; - m_cols = cols; - m_info = Success; - m_isInitialized = false; - m_isAllocated = true; - m_computationOptions = computationOptions; - m_computeFullU = (computationOptions & ComputeFullU) != 0; - m_computeThinU = (computationOptions & ComputeThinU) != 0; - m_computeFullV = (computationOptions & ComputeFullV) != 0; - m_computeThinV = (computationOptions & ComputeThinV) != 0; - eigen_assert(!(m_computeFullU && m_computeThinU) && "JacobiSVD: you can't ask for both full and thin U"); - eigen_assert(!(m_computeFullV && m_computeThinV) && "JacobiSVD: you can't ask for both full and thin V"); - eigen_assert(internal::check_implication(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && - "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns."); - if (QRPreconditioner == FullPivHouseholderQRPreconditioner) - { - eigen_assert(!(m_computeThinU || m_computeThinV) && - "JacobiSVD: can't compute thin U or thin V with the FullPivHouseholderQR preconditioner. " - "Use the ColPivHouseholderQR preconditioner instead."); - } - m_diagSize = (std::min)(m_rows, m_cols); - m_singularValues.resize(m_diagSize); - if(RowsAtCompileTime==Dynamic) - m_matrixU.resize(m_rows, m_computeFullU ? m_rows - : m_computeThinU ? m_diagSize - : 0); - if(ColsAtCompileTime==Dynamic) - m_matrixV.resize(m_cols, m_computeFullV ? m_cols - : m_computeThinV ? m_diagSize - : 0); 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); } -template -JacobiSVD& -JacobiSVD::compute(const MatrixType& matrix, unsigned int computationOptions) -{ +template +JacobiSVD& JacobiSVD::compute_impl(const MatrixType& matrix, + unsigned int computationOptions) { using std::abs; + allocate(matrix.rows(), matrix.cols(), computationOptions); // currently we stop when we reach precision 2*epsilon as the last bit of precision can require an unreasonable number of iterations, @@ -726,8 +725,8 @@ JacobiSVD::compute(const MatrixType& matrix, unsig finished = false; // perform SVD decomposition of 2x2 sub-matrix corresponding to indices p,q to make it diagonal // the complex to real operation returns true if the updated 2x2 block is not already diagonal - if(internal::svd_precondition_2x2_block_to_be_real::run(m_workMatrix, *this, p, q, maxDiagEntry)) - { + if (internal::svd_precondition_2x2_block_to_be_real::run(m_workMatrix, *this, p, q, + maxDiagEntry)) { JacobiRotation j_left, j_right; internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right); @@ -802,13 +801,19 @@ JacobiSVD::compute(const MatrixType& matrix, unsig * * \sa class JacobiSVD */ -template -JacobiSVD::PlainObject> -MatrixBase::jacobiSvd(unsigned int computationOptions) const -{ - return JacobiSVD(*this, computationOptions); +template +template +JacobiSVD::PlainObject, Options> MatrixBase::jacobiSvd() const { + return JacobiSVD(*this); } -} // end namespace Eigen +template +template +JacobiSVD::PlainObject, Options> MatrixBase::jacobiSvd( + unsigned int computationOptions) const { + return JacobiSVD(*this, computationOptions); +} + +} // end namespace Eigen #endif // EIGEN_JACOBISVD_H diff --git a/Eigen/src/SVD/JacobiSVD_LAPACKE.h b/Eigen/src/SVD/JacobiSVD_LAPACKE.h index 611ae8ce4..688df3155 100644 --- a/Eigen/src/SVD/JacobiSVD_LAPACKE.h +++ b/Eigen/src/SVD/JacobiSVD_LAPACKE.h @@ -39,15 +39,15 @@ namespace Eigen { /** \internal Specialization for the data types supported by LAPACKe */ -#define EIGEN_LAPACKE_SVD(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW) \ +#define EIGEN_LAPACKE_SVD(EIGTYPE, LAPACKE_TYPE, LAPACKE_RTYPE, LAPACKE_PREFIX, EIGCOLROW, LAPACKE_COLROW, OPTIONS) \ template<> inline \ -JacobiSVD, ColPivHouseholderQRPreconditioner>& \ -JacobiSVD, ColPivHouseholderQRPreconditioner>::compute(const Matrix& matrix, unsigned int computationOptions) \ +JacobiSVD, OPTIONS>& \ +JacobiSVD, OPTIONS>::compute(const Matrix& matrix) \ { \ typedef Matrix MatrixType; \ /*typedef MatrixType::Scalar Scalar;*/ \ /*typedef MatrixType::RealScalar RealScalar;*/ \ - allocate(matrix.rows(), matrix.cols(), computationOptions); \ + allocate(matrix.rows(), matrix.cols()); \ \ /*const RealScalar precision = RealScalar(2) * NumTraits::epsilon();*/ \ m_nonzeroSingularValues = m_diagSize; \ @@ -56,14 +56,14 @@ JacobiSVD, ColPiv lapack_int matrix_order = LAPACKE_COLROW; \ char jobu, jobvt; \ LAPACKE_TYPE *u, *vt, dummy; \ - jobu = (m_computeFullU) ? 'A' : (m_computeThinU) ? 'S' : 'N'; \ - jobvt = (m_computeFullV) ? 'A' : (m_computeThinV) ? 'S' : 'N'; \ + jobu = (ShouldComputeFullU) ? 'A' : (ShouldComputeThinU) ? 'S' : 'N'; \ + jobvt = (ShouldComputeFullV) ? 'A' : (ShouldComputeThinV) ? 'S' : 'N'; \ if (computeU()) { \ ldu = internal::convert_index(m_matrixU.outerStride()); \ u = (LAPACKE_TYPE*)m_matrixU.data(); \ } else { ldu=1; u=&dummy; }\ MatrixType localV; \ - lapack_int vt_rows = (m_computeFullV) ? internal::convert_index(m_cols) : (m_computeThinV) ? internal::convert_index(m_diagSize) : 1; \ + lapack_int vt_rows = (ShouldComputeFullV) ? internal::convert_index(m_cols) : (ShouldComputeThinV) ? internal::convert_index(m_diagSize) : 1; \ if (computeV()) { \ localV.resize(vt_rows, m_cols); \ ldvt = internal::convert_index(localV.outerStride()); \ @@ -78,15 +78,26 @@ JacobiSVD, ColPiv return *this; \ } -EIGEN_LAPACKE_SVD(double, double, double, d, ColMajor, LAPACK_COL_MAJOR) -EIGEN_LAPACKE_SVD(float, float, float , s, ColMajor, LAPACK_COL_MAJOR) -EIGEN_LAPACKE_SVD(dcomplex, lapack_complex_double, double, z, ColMajor, LAPACK_COL_MAJOR) -EIGEN_LAPACKE_SVD(scomplex, lapack_complex_float, float , c, ColMajor, LAPACK_COL_MAJOR) +#define EIGEN_LAPACK_SVD_OPTIONS(OPTIONS) \ + EIGEN_LAPACKE_SVD(double, double, double, d, ColMajor, LAPACK_COL_MAJOR, OPTIONS) \ + EIGEN_LAPACKE_SVD(float, float, float , s, ColMajor, LAPACK_COL_MAJOR, OPTIONS) \ + EIGEN_LAPACKE_SVD(dcomplex, lapack_complex_double, double, z, ColMajor, LAPACK_COL_MAJOR, OPTIONS) \ + EIGEN_LAPACKE_SVD(scomplex, lapack_complex_float, float , c, ColMajor, LAPACK_COL_MAJOR, OPTIONS) \ +\ + EIGEN_LAPACKE_SVD(double, double, double, d, RowMajor, LAPACK_ROW_MAJOR, OPTIONS) \ + EIGEN_LAPACKE_SVD(float, float, float , s, RowMajor, LAPACK_ROW_MAJOR, OPTIONS) \ + EIGEN_LAPACKE_SVD(dcomplex, lapack_complex_double, double, z, RowMajor, LAPACK_ROW_MAJOR, OPTIONS) \ + EIGEN_LAPACKE_SVD(scomplex, lapack_complex_float, float , c, RowMajor, LAPACK_ROW_MAJOR, OPTIONS) -EIGEN_LAPACKE_SVD(double, double, double, d, RowMajor, LAPACK_ROW_MAJOR) -EIGEN_LAPACKE_SVD(float, float, float , s, RowMajor, LAPACK_ROW_MAJOR) -EIGEN_LAPACKE_SVD(dcomplex, lapack_complex_double, double, z, RowMajor, LAPACK_ROW_MAJOR) -EIGEN_LAPACKE_SVD(scomplex, lapack_complex_float, float , c, RowMajor, LAPACK_ROW_MAJOR) +EIGEN_LAPACK_SVD_OPTIONS(0) +EIGEN_LAPACK_SVD_OPTIONS(ComputeThinU) +EIGEN_LAPACK_SVD_OPTIONS(ComputeThinV) +EIGEN_LAPACK_SVD_OPTIONS(ComputeFullU) +EIGEN_LAPACK_SVD_OPTIONS(ComputeFullV) +EIGEN_LAPACK_SVD_OPTIONS(ComputeThinU | ComputeThinV) +EIGEN_LAPACK_SVD_OPTIONS(ComputeFullU | ComputeFullV) +EIGEN_LAPACK_SVD_OPTIONS(ComputeThinU | ComputeFullV) +EIGEN_LAPACK_SVD_OPTIONS(ComputeFullU | ComputeThinV) } // end namespace Eigen diff --git a/Eigen/src/SVD/SVDBase.h b/Eigen/src/SVD/SVDBase.h index 1e584048f..af716faac 100644 --- a/Eigen/src/SVD/SVDBase.h +++ b/Eigen/src/SVD/SVDBase.h @@ -21,6 +21,36 @@ namespace Eigen { namespace internal { + +enum OptionsMasks { + QRPreconditionerBits = NoQRPreconditioner | HouseholderQRPreconditioner | ColPivHouseholderQRPreconditioner | + FullPivHouseholderQRPreconditioner, + ComputationOptionsBits = ComputeThinU | ComputeFullU | ComputeThinV | ComputeFullV +}; + +constexpr int get_qr_preconditioner(int options) { return options & QRPreconditionerBits; } + +constexpr int get_computation_options(int options) { return options & ComputationOptionsBits; } + +constexpr int should_svd_compute_thin_u(int options) { return options & ComputeThinU; } +constexpr int should_svd_compute_full_u(int options) { return options & ComputeFullU; } +constexpr int should_svd_compute_thin_v(int options) { return options & ComputeThinV; } +constexpr int should_svd_compute_full_v(int options) { return options & ComputeFullV; } + +template +void check_svd_constructor_assertions(unsigned int computationOptions) { + EIGEN_STATIC_ASSERT((Options & ComputationOptionsBits) == 0, + "SVDBase: Cannot request U or V using both static and runtime options, even if they match. " + "Requesting unitaries at runtime through the constructor is DEPRECATED: " + "If possible, prefer requesting unitaries statically, using the Options template parameter."); + eigen_assert( + !(should_svd_compute_thin_u(computationOptions) && MatrixType::ColsAtCompileTime != Dynamic) && + !(should_svd_compute_thin_v(computationOptions) && MatrixType::ColsAtCompileTime != Dynamic) && + "SVDBase: If U or V are requested at runtime through the constructor, then thin U and V are only available when " + "your matrix has a dynamic number of columns."); + (void)computationOptions; +} + template struct traits > : traits { @@ -29,6 +59,32 @@ template struct traits > typedef int StorageIndex; enum { Flags = 0 }; }; + +template +struct svd_traits : traits { + static constexpr bool ShouldComputeFullU = internal::should_svd_compute_full_u(Options); + static constexpr bool ShouldComputeThinU = internal::should_svd_compute_thin_u(Options); + static constexpr bool ShouldComputeFullV = internal::should_svd_compute_full_v(Options); + static constexpr bool ShouldComputeThinV = internal::should_svd_compute_thin_v(Options); + enum { + DiagSizeAtCompileTime = + internal::min_size_prefer_dynamic(MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime), + MaxDiagSizeAtCompileTime = + internal::min_size_prefer_dynamic(MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime), + MatrixUColsAtCompileTime = ShouldComputeFullU ? MatrixType::RowsAtCompileTime + : ShouldComputeThinU ? DiagSizeAtCompileTime + : Dynamic, + MatrixVColsAtCompileTime = ShouldComputeFullV ? MatrixType::ColsAtCompileTime + : ShouldComputeThinV ? DiagSizeAtCompileTime + : Dynamic, + MatrixUMaxColsAtCompileTime = ShouldComputeFullU ? MatrixType::MaxRowsAtCompileTime + : ShouldComputeThinU ? MaxDiagSizeAtCompileTime + : Dynamic, + MatrixVMaxColsAtCompileTime = ShouldComputeFullV ? MatrixType::MaxColsAtCompileTime + : ShouldComputeThinV ? MaxDiagSizeAtCompileTime + : Dynamic + }; +}; } /** \ingroup SVD_Module @@ -74,20 +130,38 @@ public: typedef typename NumTraits::Real RealScalar; typedef typename Eigen::internal::traits::StorageIndex StorageIndex; typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 + + static constexpr bool ShouldComputeFullU = internal::traits::ShouldComputeFullU; + static constexpr bool ShouldComputeThinU = internal::traits::ShouldComputeThinU; + static constexpr bool ShouldComputeFullV = internal::traits::ShouldComputeFullV; + static constexpr bool ShouldComputeThinV = internal::traits::ShouldComputeThinV; + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime, ColsAtCompileTime = MatrixType::ColsAtCompileTime, - DiagSizeAtCompileTime = internal::min_size_prefer_dynamic(RowsAtCompileTime,ColsAtCompileTime), + DiagSizeAtCompileTime = internal::min_size_prefer_dynamic(RowsAtCompileTime, ColsAtCompileTime), MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, - MaxDiagSizeAtCompileTime = internal::min_size_prefer_fixed(MaxRowsAtCompileTime,MaxColsAtCompileTime), - MatrixOptions = MatrixType::Options + MaxDiagSizeAtCompileTime = internal::min_size_prefer_fixed(MaxRowsAtCompileTime, MaxColsAtCompileTime), + MatrixOptions = MatrixType::Options, + MatrixUColsAtCompileTime = internal::traits::MatrixUColsAtCompileTime, + MatrixVColsAtCompileTime = internal::traits::MatrixVColsAtCompileTime, + MatrixUMaxColsAtCompileTime = internal::traits::MatrixUMaxColsAtCompileTime, + MatrixVMaxColsAtCompileTime = internal::traits::MatrixVMaxColsAtCompileTime }; - typedef Matrix MatrixUType; - typedef Matrix MatrixVType; + EIGEN_STATIC_ASSERT(!(ShouldComputeFullU && ShouldComputeThinU), "SVDBase: Cannot request both full and thin U") + EIGEN_STATIC_ASSERT(!(ShouldComputeFullV && ShouldComputeThinV), "SVDBase: Cannot request both full and thin V") + + typedef + typename internal::make_proper_matrix_type::type MatrixUType; + typedef + typename internal::make_proper_matrix_type::type MatrixVType; + typedef typename internal::plain_diag_type::type SingularValuesType; - + Derived& derived() { return *static_cast(this); } const Derived& derived() const { return *static_cast(this); } @@ -266,7 +340,7 @@ protected: } // return true if already allocated - bool allocate(Index rows, Index cols, unsigned int computationOptions) ; + bool allocate(Index rows, Index cols, unsigned int computationOptions); MatrixUType m_matrixU; MatrixVType m_matrixV; @@ -284,19 +358,18 @@ protected: * Default constructor of SVDBase */ SVDBase() - : m_info(Success), - m_isInitialized(false), - m_isAllocated(false), - m_usePrescribedThreshold(false), - m_computeFullU(false), - m_computeThinU(false), - m_computeFullV(false), - m_computeThinV(false), - m_computationOptions(0), - m_rows(-1), m_cols(-1), m_diagSize(0) - { } - - + : m_info(Success), + m_isInitialized(false), + m_isAllocated(false), + m_usePrescribedThreshold(false), + m_computeFullU(false), + m_computeThinU(false), + m_computeFullV(false), + m_computeThinV(false), + m_computationOptions(0), + m_rows(-1), + m_cols(-1), + m_diagSize(0) {} }; #ifndef EIGEN_PARSED_BY_DOXYGEN @@ -330,9 +403,8 @@ void SVDBase::_solve_impl_transposed(const RhsType &rhs, DstType &dst) } #endif -template -bool SVDBase::allocate(Index rows, Index cols, unsigned int computationOptions) -{ +template +bool SVDBase::allocate(Index rows, Index cols, unsigned int computationOptions) { eigen_assert(rows >= 0 && cols >= 0); if (m_isAllocated && @@ -349,14 +421,13 @@ bool SVDBase::allocate(Index rows, Index cols, unsigned int computat m_isInitialized = false; m_isAllocated = true; m_computationOptions = computationOptions; - m_computeFullU = (computationOptions & ComputeFullU) != 0; - m_computeThinU = (computationOptions & ComputeThinU) != 0; - m_computeFullV = (computationOptions & ComputeFullV) != 0; - m_computeThinV = (computationOptions & ComputeThinV) != 0; + m_computeFullU = ShouldComputeFullU || internal::should_svd_compute_full_u(computationOptions); + m_computeThinU = ShouldComputeThinU || internal::should_svd_compute_thin_u(computationOptions); + m_computeFullV = ShouldComputeFullV || internal::should_svd_compute_full_v(computationOptions); + m_computeThinV = ShouldComputeThinV || internal::should_svd_compute_thin_v(computationOptions); + 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"); - eigen_assert(internal::check_implication(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && - "SVDBase: thin U and V are only available when your matrix has a dynamic number of columns."); m_diagSize = (std::min)(m_rows, m_cols); m_singularValues.resize(m_diagSize); diff --git a/bench/dense_solvers.cpp b/bench/dense_solvers.cpp index 24343dcd8..11c755b30 100644 --- a/bench/dense_solvers.cpp +++ b/bench/dense_solvers.cpp @@ -38,8 +38,6 @@ void bench(int id, int rows, int size = Size) A = A*A.adjoint(); BenchTimer t_llt, t_ldlt, t_lu, t_fplu, t_qr, t_cpqr, t_cod, t_fpqr, t_jsvd, t_bdcsvd; - int svd_opt = ComputeThinU|ComputeThinV; - int tries = 5; int rep = 1000/size; if(rep==0) rep = 1; @@ -53,8 +51,8 @@ void bench(int id, int rows, int size = Size) ColPivHouseholderQR cpqr(A.rows(),A.cols()); CompleteOrthogonalDecomposition cod(A.rows(),A.cols()); FullPivHouseholderQR fpqr(A.rows(),A.cols()); - JacobiSVD jsvd(A.rows(),A.cols()); - BDCSVD bdcsvd(A.rows(),A.cols()); + JacobiSVD jsvd(A.rows(),A.cols()); + BDCSVD bdcsvd(A.rows(),A.cols()); BENCH(t_llt, tries, rep, compute_norm_equation(llt,A)); BENCH(t_ldlt, tries, rep, compute_norm_equation(ldlt,A)); @@ -67,9 +65,9 @@ void bench(int id, int rows, int size = Size) if(size*rows<=10000000) BENCH(t_fpqr, tries, rep, compute(fpqr,A)); if(size<500) // JacobiSVD is really too slow for too large matrices - BENCH(t_jsvd, tries, rep, jsvd.compute(A,svd_opt)); + BENCH(t_jsvd, tries, rep, jsvd.compute(A)); // if(size*rows<=20000000) - BENCH(t_bdcsvd, tries, rep, bdcsvd.compute(A,svd_opt)); + BENCH(t_bdcsvd, tries, rep, bdcsvd.compute(A)); results["LLT"][id] = t_llt.best(); results["LDLT"][id] = t_ldlt.best(); diff --git a/doc/UsingBlasLapackBackends.dox b/doc/UsingBlasLapackBackends.dox index caa597122..457282056 100644 --- a/doc/UsingBlasLapackBackends.dox +++ b/doc/UsingBlasLapackBackends.dox @@ -101,8 +101,8 @@ m1.colPivHouseholderQr(); ?geqp3 \endcode Singular value decomposition \n \c EIGEN_USE_LAPACKE \code -JacobiSVD svd; -svd.compute(m1, ComputeThinV); +JacobiSVD svd; +svd.compute(m1); \endcode\code ?gesvd \endcode diff --git a/doc/examples/TutorialLinAlgSVDSolve.cpp b/doc/examples/TutorialLinAlgSVDSolve.cpp index 23ad42295..ca2e78c0c 100644 --- a/doc/examples/TutorialLinAlgSVDSolve.cpp +++ b/doc/examples/TutorialLinAlgSVDSolve.cpp @@ -3,10 +3,10 @@ int main() { - Eigen::MatrixXf A = Eigen::MatrixXf::Random(3, 2); - std::cout << "Here is the matrix A:\n" << A << std::endl; - Eigen::VectorXf b = Eigen::VectorXf::Random(3); - std::cout << "Here is the right hand side b:\n" << b << std::endl; - std::cout << "The least-squares solution is:\n" - << A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b) << std::endl; + MatrixXf A = MatrixXf::Random(3, 2); + cout << "Here is the matrix A:\n" << A << endl; + VectorXf b = VectorXf::Random(3); + cout << "Here is the right hand side b:\n" << b << endl; + cout << "The least-squares solution is:\n" + << A.template bdcSvd().solve(b) << endl; } diff --git a/doc/snippets/JacobiSVD_basic.cpp b/doc/snippets/JacobiSVD_basic.cpp index ab24b9bca..6c21bafc2 100644 --- a/doc/snippets/JacobiSVD_basic.cpp +++ b/doc/snippets/JacobiSVD_basic.cpp @@ -1,6 +1,6 @@ MatrixXf m = MatrixXf::Random(3,2); cout << "Here is the matrix m:" << endl << m << endl; -JacobiSVD svd(m, ComputeThinU | ComputeThinV); +JacobiSVD svd(m); cout << "Its singular values are:" << endl << svd.singularValues() << endl; cout << "Its left singular vectors are the columns of the thin U matrix:" << endl << svd.matrixU() << endl; cout << "Its right singular vectors are the columns of the thin V matrix:" << endl << svd.matrixV() << endl; diff --git a/lapack/svd.cpp b/lapack/svd.cpp index 77b302b6b..83544cf32 100644 --- a/lapack/svd.cpp +++ b/lapack/svd.cpp @@ -135,4 +135,4 @@ EIGEN_LAPACK_FUNC(gesvd,(char *jobu, char *jobv, int *m, int* n, Scalar* a, int else if(*jobv=='O') matrix(a,diag_size,*n,*lda) = svd.matrixV().adjoint(); } return 0; -} +} \ No newline at end of file diff --git a/test/bdcsvd.cpp b/test/bdcsvd.cpp index f98cdcaf6..6eb82f09f 100644 --- a/test/bdcsvd.cpp +++ b/test/bdcsvd.cpp @@ -16,29 +16,12 @@ #include "main.h" #include -#include -#include - #define SVD_DEFAULT(M) BDCSVD #define SVD_FOR_MIN_NORM(M) BDCSVD +#define SVD_STATIC_OPTIONS(M, O) BDCSVD #include "svd_common.h" -// Check all variants of JacobiSVD -template -void bdcsvd(const MatrixType& a = MatrixType(), bool pickrandom = true) -{ - MatrixType m; - if(pickrandom) { - m.resizeLike(a); - svd_fill_random(m); - } - else - m = a; - - CALL_SUBTEST(( svd_test_all_computation_options >(m, false) )); -} - template void bdcsvd_method() { @@ -52,25 +35,22 @@ void bdcsvd_method() VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).solve(m), m); VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).transpose().solve(m), m); VERIFY_IS_APPROX(m.bdcSvd(ComputeFullU|ComputeFullV).adjoint().solve(m), m); + VERIFY_IS_APPROX(m.template bdcSvd().solve(m), m); + VERIFY_IS_APPROX(m.template bdcSvd().transpose().solve(m), m); + VERIFY_IS_APPROX(m.template bdcSvd().adjoint().solve(m), m); } -// Compare the Singular values returned with Jacobi and Bdc. -template -void compare_bdc_jacobi(const MatrixType& a = MatrixType(), unsigned int computationOptions = 0, int algoswap = 16, bool random = true) -{ +// compare the Singular values returned with Jacobi and Bdc +template +void compare_bdc_jacobi(const MatrixType& a = MatrixType(), int algoswap = 16, bool random = true) { MatrixType m = random ? MatrixType::Random(a.rows(), a.cols()) : a; - BDCSVD bdc_svd(m.rows(), m.cols(), computationOptions); + BDCSVD bdc_svd(m.rows(), m.cols()); bdc_svd.setSwitchSize(algoswap); bdc_svd.compute(m); JacobiSVD jacobi_svd(m); VERIFY_IS_APPROX(bdc_svd.singularValues(), jacobi_svd.singularValues()); - - if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); - if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(bdc_svd.matrixU(), jacobi_svd.matrixU()); - if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); - if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(bdc_svd.matrixV(), jacobi_svd.matrixV()); } // Verifies total deflation is **not** triggered. @@ -91,41 +71,72 @@ void compare_bdc_jacobi_instance(bool structure_as_m, int algoswap = 16) -20.794, 8.68496, -4.83103, -8.4981, -10.5451, 23.9072; } - compare_bdc_jacobi(m, 0, algoswap, false); + compare_bdc_jacobi(m, algoswap, false); +} + +template +void bdcsvd_all_options(const MatrixType& input = MatrixType()) { + MatrixType m = input; + svd_fill_random(m); + svd_option_checks(m); +} + +template +void bdcsvd_verify_assert(const MatrixType& input = MatrixType()) { + svd_verify_assert(input); + svd_verify_constructor_options_assert>(input); } EIGEN_DECLARE_TEST(bdcsvd) { - CALL_SUBTEST_3(( svd_verify_assert >(Matrix3f()) )); - CALL_SUBTEST_4(( svd_verify_assert >(Matrix4d()) )); - CALL_SUBTEST_7(( svd_verify_assert >(MatrixXf(10,12)) )); - CALL_SUBTEST_8(( svd_verify_assert >(MatrixXcd(7,5)) )); + CALL_SUBTEST_3((bdcsvd_verify_assert())); + CALL_SUBTEST_4((bdcsvd_verify_assert())); + CALL_SUBTEST_7((bdcsvd_verify_assert>())); + CALL_SUBTEST_7((bdcsvd_verify_assert>())); + CALL_SUBTEST_9((bdcsvd_verify_assert, 20, 27>>())); - CALL_SUBTEST_101(( svd_all_trivial_2x2(bdcsvd) )); - CALL_SUBTEST_102(( svd_all_trivial_2x2(bdcsvd) )); - - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_3(( bdcsvd() )); - CALL_SUBTEST_4(( bdcsvd() )); - CALL_SUBTEST_5(( bdcsvd >() )); + CALL_SUBTEST_101((svd_all_trivial_2x2(bdcsvd_all_options))); + CALL_SUBTEST_102((svd_all_trivial_2x2(bdcsvd_all_options))); + for (int i = 0; i < g_repeat; i++) { int r = internal::random(1, EIGEN_TEST_MAX_SIZE/2), c = internal::random(1, EIGEN_TEST_MAX_SIZE/2); TEST_SET_BUT_UNUSED_VARIABLE(r) TEST_SET_BUT_UNUSED_VARIABLE(c) - CALL_SUBTEST_6(( bdcsvd(Matrix(r,2)) )); - CALL_SUBTEST_7(( bdcsvd(MatrixXf(r,c)) )); - CALL_SUBTEST_7(( compare_bdc_jacobi(MatrixXf(r,c)) )); - CALL_SUBTEST_10(( bdcsvd(MatrixXd(r,c)) )); - CALL_SUBTEST_10(( compare_bdc_jacobi(MatrixXd(r,c)) )); - CALL_SUBTEST_8(( bdcsvd(MatrixXcd(r,c)) )); - CALL_SUBTEST_8(( compare_bdc_jacobi(MatrixXcd(r,c)) )); - + CALL_SUBTEST_7((compare_bdc_jacobi(MatrixXf(r, c)))); + CALL_SUBTEST_10((compare_bdc_jacobi(MatrixXd(r, c)))); + CALL_SUBTEST_8((compare_bdc_jacobi(MatrixXcd(r, c)))); // Test on inf/nan matrix - CALL_SUBTEST_7( (svd_inf_nan, MatrixXf>()) ); - CALL_SUBTEST_10( (svd_inf_nan, MatrixXd>()) ); + CALL_SUBTEST_7((svd_inf_nan())); + CALL_SUBTEST_10((svd_inf_nan())); + + // Verify some computations using all combinations of the Options template parameter. + CALL_SUBTEST_3((bdcsvd_all_options())); + CALL_SUBTEST_3((bdcsvd_all_options>())); + CALL_SUBTEST_4((bdcsvd_all_options>())); + CALL_SUBTEST_4((bdcsvd_all_options>())); + CALL_SUBTEST_5((bdcsvd_all_options>(Matrix(r, 30)))); + CALL_SUBTEST_5((bdcsvd_all_options>(Matrix(20, c)))); + CALL_SUBTEST_7((bdcsvd_all_options(MatrixXf(r, c)))); + CALL_SUBTEST_8((bdcsvd_all_options(MatrixXcd(r, c)))); + CALL_SUBTEST_10((bdcsvd_all_options(MatrixXd(r, c)))); + CALL_SUBTEST_14((bdcsvd_all_options>())); + CALL_SUBTEST_14((bdcsvd_all_options>())); + + CALL_SUBTEST_15(( + svd_check_max_size_matrix, ColPivHouseholderQRPreconditioner>( + r, c))); + CALL_SUBTEST_15( + (svd_check_max_size_matrix, HouseholderQRPreconditioner>(r, + c))); + CALL_SUBTEST_15(( + svd_check_max_size_matrix, ColPivHouseholderQRPreconditioner>( + r, c))); + CALL_SUBTEST_15( + (svd_check_max_size_matrix, HouseholderQRPreconditioner>(r, + c))); } // test matrixbase method diff --git a/test/boostmultiprec.cpp b/test/boostmultiprec.cpp index e83e97044..e2fc9a811 100644 --- a/test/boostmultiprec.cpp +++ b/test/boostmultiprec.cpp @@ -200,8 +200,8 @@ EIGEN_DECLARE_TEST(boostmultiprec) TEST_SET_BUT_UNUSED_VARIABLE(s) } - CALL_SUBTEST_9(( jacobisvd(Mat(internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); - CALL_SUBTEST_10(( bdcsvd(Mat(internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); + CALL_SUBTEST_9(( jacobisvd_all_options(Mat(internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); + CALL_SUBTEST_10(( bdcsvd_all_options(Mat(internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); CALL_SUBTEST_11(( test_simplicial_cholesky_T() )); } diff --git a/test/diagonalmatrices.cpp b/test/diagonalmatrices.cpp index a88c3c124..76c8886e6 100644 --- a/test/diagonalmatrices.cpp +++ b/test/diagonalmatrices.cpp @@ -96,9 +96,9 @@ template void diagonalmatrices(const MatrixType& m) // products do not allocate memory MatrixType res(rows, cols); internal::set_is_malloc_allowed(false); - res.noalias() = ldm1 * m; - res.noalias() = m * rdm1; - res.noalias() = ldm1 * m * rdm1; + res.noalias() = ldm1 * m1; + res.noalias() = m1 * rdm1; + res.noalias() = ldm1 * m1 * rdm1; internal::set_is_malloc_allowed(true); // scalar multiple diff --git a/test/geo_alignedbox.cpp b/test/geo_alignedbox.cpp index 7b1684f29..e4dab32e0 100644 --- a/test/geo_alignedbox.cpp +++ b/test/geo_alignedbox.cpp @@ -211,7 +211,7 @@ MatrixType randomRotationMatrix() // https://www.isprs-ann-photogramm-remote-sens-spatial-inf-sci.net/III-7/103/2016/isprs-annals-III-7-103-2016.pdf const MatrixType rand = MatrixType::Random(); const MatrixType q = rand.householderQr().householderQ(); - const JacobiSVD svd = q.jacobiSvd(ComputeFullU | ComputeFullV); + const JacobiSVD svd(q); const typename MatrixType::Scalar det = (svd.matrixU() * svd.matrixV().transpose()).determinant(); MatrixType diag = rand.Identity(); diag(MatrixType::RowsAtCompileTime - 1, MatrixType::ColsAtCompileTime - 1) = det; diff --git a/test/jacobisvd.cpp b/test/jacobisvd.cpp index 5b15c5a27..e9d0cc835 100644 --- a/test/jacobisvd.cpp +++ b/test/jacobisvd.cpp @@ -16,49 +16,9 @@ #define SVD_DEFAULT(M) JacobiSVD #define SVD_FOR_MIN_NORM(M) JacobiSVD +#define SVD_STATIC_OPTIONS(M, O) JacobiSVD #include "svd_common.h" -// Check all variants of JacobiSVD -template -void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) -{ - MatrixType m = a; - if(pickrandom) - svd_fill_random(m); - - CALL_SUBTEST(( svd_test_all_computation_options >(m, true) )); // check full only - CALL_SUBTEST(( svd_test_all_computation_options >(m, false) )); - CALL_SUBTEST(( svd_test_all_computation_options >(m, false) )); - if(m.rows()==m.cols()) - CALL_SUBTEST(( svd_test_all_computation_options >(m, false) )); -} - -template void jacobisvd_verify_assert(const MatrixType& m) -{ - svd_verify_assert >(m); - svd_verify_assert >(m, true); - svd_verify_assert >(m); - svd_verify_assert >(m); - Index rows = m.rows(); - Index cols = m.cols(); - - enum { - ColsAtCompileTime = MatrixType::ColsAtCompileTime - }; - - - MatrixType a = MatrixType::Zero(rows, cols); - a.setZero(); - - if (ColsAtCompileTime == Dynamic) - { - JacobiSVD svd_fullqr; - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeFullU|ComputeThinV)) - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeThinV)) - VERIFY_RAISES_ASSERT(svd_fullqr.compute(a, ComputeThinU|ComputeFullV)) - } -} - template void jacobisvd_method() { @@ -72,6 +32,55 @@ void jacobisvd_method() VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).solve(m), m); VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).transpose().solve(m), m); VERIFY_IS_APPROX(m.jacobiSvd(ComputeFullU|ComputeFullV).adjoint().solve(m), m); + VERIFY_IS_APPROX(m.template jacobiSvd().solve(m), m); + VERIFY_IS_APPROX(m.template jacobiSvd().transpose().solve(m), m); + VERIFY_IS_APPROX(m.template jacobiSvd().adjoint().solve(m), m); +} + +template +void jacobisvd_all_options(const MatrixType& input = MatrixType()) { + MatrixType m = input; + svd_fill_random(m); + svd_option_checks(m); + svd_option_checks(m); + svd_option_checks(m); + svd_option_checks_full_only( + m); // FullPiv only used when computing full unitaries +} + +template +void jacobisvd_verify_assert(const MatrixType& m = MatrixType()) { + svd_verify_assert(m); + svd_verify_assert(m); + svd_verify_assert(m); + svd_verify_assert_full_only(m); + + svd_verify_constructor_options_assert>(m); + svd_verify_constructor_options_assert>(m); + svd_verify_constructor_options_assert>(m); + svd_verify_constructor_options_assert>(m, true); +} + +template +void jacobisvd_verify_inputs(const MatrixType& m = MatrixType()) { + // check defaults + typedef JacobiSVD DefaultSVD; + DefaultSVD defaultSvd(m); + VERIFY((int)DefaultSVD::QRPreconditioner == (int)ColPivHouseholderQRPreconditioner); + VERIFY(!defaultSvd.computeU()); + VERIFY(!defaultSvd.computeV()); + + // ColPivHouseholderQR is always default in presence of other options. + VERIFY(((int)JacobiSVD::QRPreconditioner == (int)ColPivHouseholderQRPreconditioner)); + VERIFY(((int)JacobiSVD::QRPreconditioner == (int)ColPivHouseholderQRPreconditioner)); + VERIFY(((int)JacobiSVD::QRPreconditioner == + (int)ColPivHouseholderQRPreconditioner)); + VERIFY(((int)JacobiSVD::QRPreconditioner == + (int)ColPivHouseholderQRPreconditioner)); + VERIFY(((int)JacobiSVD::QRPreconditioner == + (int)ColPivHouseholderQRPreconditioner)); + VERIFY(((int)JacobiSVD::QRPreconditioner == + (int)ColPivHouseholderQRPreconditioner)); } namespace Foo { @@ -91,45 +100,73 @@ void msvc_workaround() EIGEN_DECLARE_TEST(jacobisvd) { - CALL_SUBTEST_3(( jacobisvd_verify_assert(Matrix3f()) )); - CALL_SUBTEST_4(( jacobisvd_verify_assert(Matrix4d()) )); - CALL_SUBTEST_7(( jacobisvd_verify_assert(MatrixXf(10,12)) )); - CALL_SUBTEST_8(( jacobisvd_verify_assert(MatrixXcd(7,5)) )); - - CALL_SUBTEST_11(svd_all_trivial_2x2(jacobisvd)); - CALL_SUBTEST_12(svd_all_trivial_2x2(jacobisvd)); + CALL_SUBTEST_4((jacobisvd_verify_inputs())); + CALL_SUBTEST_7((jacobisvd_verify_inputs(Matrix(10, 12)))); + CALL_SUBTEST_8((jacobisvd_verify_inputs, 7, 5>>())); - for(int i = 0; i < g_repeat; i++) { - CALL_SUBTEST_3(( jacobisvd() )); - CALL_SUBTEST_4(( jacobisvd() )); - CALL_SUBTEST_5(( jacobisvd >() )); - CALL_SUBTEST_6(( jacobisvd >(Matrix(10,2)) )); + CALL_SUBTEST_3((jacobisvd_verify_assert())); + CALL_SUBTEST_4((jacobisvd_verify_assert())); + CALL_SUBTEST_7((jacobisvd_verify_assert>())); + CALL_SUBTEST_7((jacobisvd_verify_assert>())); + CALL_SUBTEST_7((jacobisvd_verify_assert(MatrixXf(10, 12)))); + CALL_SUBTEST_8((jacobisvd_verify_assert(MatrixXcd(7, 5)))); + CALL_SUBTEST_11(svd_all_trivial_2x2(jacobisvd_all_options)); + CALL_SUBTEST_12(svd_all_trivial_2x2(jacobisvd_all_options)); + + for (int i = 0; i < g_repeat; i++) { int r = internal::random(1, 30), c = internal::random(1, 30); TEST_SET_BUT_UNUSED_VARIABLE(r) TEST_SET_BUT_UNUSED_VARIABLE(c) - - CALL_SUBTEST_10(( jacobisvd(MatrixXd(r,c)) )); - CALL_SUBTEST_7(( jacobisvd(MatrixXf(r,c)) )); - CALL_SUBTEST_8(( jacobisvd(MatrixXcd(r,c)) )); - (void) r; - (void) c; + + CALL_SUBTEST_3((jacobisvd_all_options())); + CALL_SUBTEST_3((jacobisvd_all_options>())); + CALL_SUBTEST_4((jacobisvd_all_options())); + CALL_SUBTEST_4((jacobisvd_all_options>())); + CALL_SUBTEST_4((jacobisvd_all_options>())); + CALL_SUBTEST_5((jacobisvd_all_options>(Matrix(r, 16)))); + CALL_SUBTEST_5((jacobisvd_all_options>(Matrix(10, c)))); + CALL_SUBTEST_7((jacobisvd_all_options(MatrixXf(r, c)))); + CALL_SUBTEST_8((jacobisvd_all_options(MatrixXcd(r, c)))); + CALL_SUBTEST_10((jacobisvd_all_options(MatrixXd(r, c)))); + CALL_SUBTEST_14((jacobisvd_all_options>())); + CALL_SUBTEST_14((jacobisvd_all_options>())); + + MatrixXcd noQRTest = MatrixXcd(r, r); + svd_fill_random(noQRTest); + CALL_SUBTEST_16((svd_option_checks(noQRTest))); + + CALL_SUBTEST_15(( + svd_check_max_size_matrix, ColPivHouseholderQRPreconditioner>( + r, c))); + CALL_SUBTEST_15( + (svd_check_max_size_matrix, HouseholderQRPreconditioner>(r, + c))); + CALL_SUBTEST_15(( + svd_check_max_size_matrix, ColPivHouseholderQRPreconditioner>( + r, c))); + CALL_SUBTEST_15( + (svd_check_max_size_matrix, HouseholderQRPreconditioner>(r, + c))); // Test on inf/nan matrix - CALL_SUBTEST_7( (svd_inf_nan, MatrixXf>()) ); - CALL_SUBTEST_10( (svd_inf_nan, MatrixXd>()) ); + CALL_SUBTEST_7((svd_inf_nan())); + CALL_SUBTEST_10((svd_inf_nan())); - // bug1395 test compile-time vectors as input - CALL_SUBTEST_13(( jacobisvd_verify_assert(Matrix()) )); - CALL_SUBTEST_13(( jacobisvd_verify_assert(Matrix()) )); - CALL_SUBTEST_13(( jacobisvd_verify_assert(Matrix(r)) )); - CALL_SUBTEST_13(( jacobisvd_verify_assert(Matrix(c)) )); + CALL_SUBTEST_13((jacobisvd_verify_assert>())); + CALL_SUBTEST_13((jacobisvd_verify_assert>())); + CALL_SUBTEST_13((jacobisvd_verify_assert>(Matrix(r)))); + CALL_SUBTEST_13((jacobisvd_verify_assert>(Matrix(c)))); } - CALL_SUBTEST_7(( jacobisvd(MatrixXf(internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) )); - CALL_SUBTEST_8(( jacobisvd(MatrixXcd(internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3), internal::random(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3))) )); + CALL_SUBTEST_7((jacobisvd_all_options( + MatrixXd(internal::random(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 2), + internal::random(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 2))))); + CALL_SUBTEST_8((jacobisvd_all_options( + MatrixXcd(internal::random(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 3), + internal::random(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 3))))); // test matrixbase method CALL_SUBTEST_1(( jacobisvd_method() )); diff --git a/test/nomalloc.cpp b/test/nomalloc.cpp index cb4c073e9..689a4ccf4 100644 --- a/test/nomalloc.cpp +++ b/test/nomalloc.cpp @@ -152,7 +152,7 @@ void ctms_decompositions() x = fpQR.solve(b); // SVD module - Eigen::JacobiSVD jSVD; jSVD.compute(A, ComputeFullU | ComputeFullV); + Eigen::JacobiSVD jSVD; jSVD.compute(A); } void test_zerosized() { diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp index 06f16438f..e25167138 100644 --- a/test/qr_colpivoting.cpp +++ b/test/qr_colpivoting.cpp @@ -55,7 +55,7 @@ void cod() { MatrixType exact_solution = MatrixType::Random(cols, cols2); MatrixType rhs = matrix * exact_solution; MatrixType cod_solution = cod.solve(rhs); - JacobiSVD svd(matrix, ComputeThinU | ComputeThinV); + JacobiSVD svd(matrix); MatrixType svd_solution = svd.solve(rhs); VERIFY_IS_APPROX(cod_solution, svd_solution); @@ -88,7 +88,7 @@ void cod_fixedsize() { exact_solution.setRandom(Cols, Cols2); Matrix rhs = matrix * exact_solution; Matrix cod_solution = cod.solve(rhs); - JacobiSVD svd(matrix, ComputeFullU | ComputeFullV); + JacobiSVD svd(matrix); Matrix svd_solution = svd.solve(rhs); VERIFY_IS_APPROX(cod_solution, svd_solution); diff --git a/test/svd_common.h b/test/svd_common.h index eae4c0bfe..1514f674f 100644 --- a/test/svd_common.h +++ b/test/svd_common.h @@ -16,6 +16,10 @@ #error a macro SVD_FOR_MIN_NORM(MatrixType) must be defined prior to including svd_common.h #endif +#ifndef SVD_STATIC_OPTIONS +#error a macro SVD_STATIC_OPTIONS(MatrixType, Options) must be defined prior to including svd_common.h +#endif + #include "svd_fill.h" #include "solverbase.h" @@ -55,50 +59,44 @@ void svd_check_full(const MatrixType& m, const SvdType& svd) } // Compare partial SVD defined by computationOptions to a full SVD referenceSvd -template -void svd_compare_to_full(const MatrixType& m, - unsigned int computationOptions, - const SvdType& referenceSvd) -{ +template +void svd_compare_to_full(const MatrixType& m, const SvdType& referenceSvd) { typedef typename MatrixType::RealScalar RealScalar; Index rows = m.rows(); Index cols = m.cols(); Index diagSize = (std::min)(rows, cols); RealScalar prec = test_precision(); - SvdType svd(m, computationOptions); + SVD_STATIC_OPTIONS(MatrixType, Options) svd(m); VERIFY_IS_APPROX(svd.singularValues(), referenceSvd.singularValues()); - - if(computationOptions & (ComputeFullV|ComputeThinV)) - { + + if (Options & (ComputeFullV | ComputeThinV)) { VERIFY( (svd.matrixV().adjoint()*svd.matrixV()).isIdentity(prec) ); VERIFY_IS_APPROX( svd.matrixV().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint(), referenceSvd.matrixV().leftCols(diagSize) * referenceSvd.singularValues().asDiagonal() * referenceSvd.matrixV().leftCols(diagSize).adjoint()); } - - if(computationOptions & (ComputeFullU|ComputeThinU)) - { + + if (Options & (ComputeFullU | ComputeThinU)) { VERIFY( (svd.matrixU().adjoint()*svd.matrixU()).isIdentity(prec) ); VERIFY_IS_APPROX( svd.matrixU().leftCols(diagSize) * svd.singularValues().cwiseAbs2().asDiagonal() * svd.matrixU().leftCols(diagSize).adjoint(), referenceSvd.matrixU().leftCols(diagSize) * referenceSvd.singularValues().cwiseAbs2().asDiagonal() * referenceSvd.matrixU().leftCols(diagSize).adjoint()); } - + // The following checks are not critical. - // For instance, with Dived&Conquer SVD, if only the factor 'V' is computedt then different matrix-matrix product implementation will be used - // and the resulting 'V' factor might be significantly different when the SVD decomposition is not unique, especially with single precision float. + // For instance, with Dived&Conquer SVD, if only the factor 'V' is computed then different matrix-matrix product + // implementation will be used and the resulting 'V' factor might be significantly different when the SVD + // decomposition is not unique, especially with single precision float. ++g_test_level; - if(computationOptions & ComputeFullU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU()); - if(computationOptions & ComputeThinU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); - if(computationOptions & ComputeFullV) VERIFY_IS_APPROX(svd.matrixV().cwiseAbs(), referenceSvd.matrixV().cwiseAbs()); - if(computationOptions & ComputeThinV) VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); + if (Options & ComputeFullU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU()); + if (Options & ComputeThinU) VERIFY_IS_APPROX(svd.matrixU(), referenceSvd.matrixU().leftCols(diagSize)); + if (Options & ComputeFullV) VERIFY_IS_APPROX(svd.matrixV().cwiseAbs(), referenceSvd.matrixV().cwiseAbs()); + if (Options & ComputeThinV) VERIFY_IS_APPROX(svd.matrixV(), referenceSvd.matrixV().leftCols(diagSize)); --g_test_level; } -// -template -void svd_least_square(const MatrixType& m, unsigned int computationOptions) -{ +template +void svd_least_square(const MatrixType& m) { typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::RealScalar RealScalar; Index rows = m.rows(); @@ -113,9 +111,10 @@ void svd_least_square(const MatrixType& m, unsigned int computationOptions) typedef Matrix SolutionType; RhsType rhs = RhsType::Random(rows, internal::random(1, cols)); - SvdType svd(m, computationOptions); + SvdType svd(m); - if(internal::is_same::value) svd.setThreshold(1e-8); + if (internal::is_same::value) + svd.setThreshold(1e-8); else if(internal::is_same::value) svd.setThreshold(2e-4); SolutionType x = svd.solve(rhs); @@ -162,10 +161,9 @@ void svd_least_square(const MatrixType& m, unsigned int computationOptions) } } -// check minimal norm solutions, the inoput matrix m is only used to recover problem size -template -void svd_min_norm(const MatrixType& m, unsigned int computationOptions) -{ +// check minimal norm solutions, the input matrix m is only used to recover problem size +template +void svd_min_norm(const MatrixType& m) { typedef typename MatrixType::Scalar Scalar; Index cols = m.cols(); @@ -199,7 +197,7 @@ void svd_min_norm(const MatrixType& m, unsigned int computationOptions) tmp.tail(cols-rank).setZero(); SolutionType x21 = qr.householderQ() * tmp; // now check with SVD - SVD_FOR_MIN_NORM(MatrixType2) svd2(m2, computationOptions); + SVD_STATIC_OPTIONS(MatrixType2, Options) svd2(m2); SolutionType x22 = svd2.solve(rhs2); VERIFY_IS_APPROX(m2*x21, rhs2); VERIFY_IS_APPROX(m2*x22, rhs2); @@ -212,7 +210,7 @@ void svd_min_norm(const MatrixType& m, unsigned int computationOptions) Matrix C = Matrix::Random(rows3,rank); MatrixType3 m3 = C * m2; RhsType3 rhs3 = C * rhs2; - SVD_FOR_MIN_NORM(MatrixType3) svd3(m3, computationOptions); + SVD_STATIC_OPTIONS(MatrixType3, Options) svd3(m3); SolutionType x3 = svd3.solve(rhs3); VERIFY_IS_APPROX(m3*x3, rhs3); VERIFY_IS_APPROX(m3*x21, rhs3); @@ -239,57 +237,6 @@ void svd_test_solvers(const MatrixType& m, const SolverType& solver) { check_solverbase(m, solver, rows, cols, cols2); } -// Check full, compare_to_full, least_square, and min_norm for all possible compute-options -template -void svd_test_all_computation_options(const MatrixType& m, bool full_only) -{ -// if (QRPreconditioner == NoQRPreconditioner && m.rows() != m.cols()) -// return; - STATIC_CHECK(( internal::is_same::value )); - - SvdType fullSvd(m, ComputeFullU|ComputeFullV); - CALL_SUBTEST(( svd_check_full(m, fullSvd) )); - CALL_SUBTEST(( svd_least_square(m, ComputeFullU | ComputeFullV) )); - CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeFullV) )); - - #if defined __INTEL_COMPILER - // remark #111: statement is unreachable - #pragma warning disable 111 - #endif - - svd_test_solvers(m, fullSvd); - - if(full_only) - return; - - CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullU, fullSvd) )); - CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullV, fullSvd) )); - CALL_SUBTEST(( svd_compare_to_full(m, 0, fullSvd) )); - - if (MatrixType::ColsAtCompileTime == Dynamic) { - // thin U/V are only available with dynamic number of columns - CALL_SUBTEST(( svd_compare_to_full(m, ComputeFullU|ComputeThinV, fullSvd) )); - CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinV, fullSvd) )); - CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU|ComputeFullV, fullSvd) )); - CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU , fullSvd) )); - CALL_SUBTEST(( svd_compare_to_full(m, ComputeThinU|ComputeThinV, fullSvd) )); - - CALL_SUBTEST(( svd_least_square(m, ComputeFullU | ComputeThinV) )); - CALL_SUBTEST(( svd_least_square(m, ComputeThinU | ComputeFullV) )); - CALL_SUBTEST(( svd_least_square(m, ComputeThinU | ComputeThinV) )); - - CALL_SUBTEST(( svd_min_norm(m, ComputeFullU | ComputeThinV) )); - CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeFullV) )); - CALL_SUBTEST(( svd_min_norm(m, ComputeThinU | ComputeThinV) )); - - // test reconstruction - Index diagSize = (std::min)(m.rows(), m.cols()); - SvdType svd(m, ComputeThinU | ComputeThinV); - VERIFY_IS_APPROX(m, svd.matrixU().leftCols(diagSize) * svd.singularValues().asDiagonal() * svd.matrixV().leftCols(diagSize).adjoint()); - } -} - - // work around stupid msvc error when constructing at compile time an expression that involves // a division by zero, even if the numeric type has floating point template @@ -300,29 +247,28 @@ template EIGEN_DONT_INLINE T sub(T a, T b) { return a - b; } // This function verifies we don't iterate infinitely on nan/inf values, // and that info() returns InvalidInput. -template -void svd_inf_nan() -{ - SvdType svd; +template +void svd_inf_nan() { + SVD_STATIC_OPTIONS(MatrixType, ComputeFullU | ComputeFullV) svd; typedef typename MatrixType::Scalar Scalar; Scalar some_inf = Scalar(1) / zero(); VERIFY(sub(some_inf, some_inf) != sub(some_inf, some_inf)); - svd.compute(MatrixType::Constant(10,10,some_inf), ComputeFullU | ComputeFullV); + svd.compute(MatrixType::Constant(10, 10, some_inf)); VERIFY(svd.info() == InvalidInput); Scalar nan = std::numeric_limits::quiet_NaN(); VERIFY(nan != nan); - svd.compute(MatrixType::Constant(10,10,nan), ComputeFullU | ComputeFullV); + svd.compute(MatrixType::Constant(10, 10, nan)); VERIFY(svd.info() == InvalidInput); MatrixType m = MatrixType::Zero(10,10); m(internal::random(0,9), internal::random(0,9)) = some_inf; - svd.compute(m, ComputeFullU | ComputeFullV); + svd.compute(m); VERIFY(svd.info() == InvalidInput); m = MatrixType::Zero(10,10); m(internal::random(0,9), internal::random(0,9)) = nan; - svd.compute(m, ComputeFullU | ComputeFullV); + svd.compute(m); VERIFY(svd.info() == InvalidInput); // regression test for bug 791 @@ -330,7 +276,7 @@ void svd_inf_nan() m << 0, 2*NumTraits::epsilon(), 0.5, 0, -0.5, 0, nan, 0, 0; - svd.compute(m, ComputeFullU | ComputeFullV); + svd.compute(m); VERIFY(svd.info() == InvalidInput); m.resize(4,4); @@ -338,7 +284,7 @@ void svd_inf_nan() 0, 3, 1, 2e-308, 1, 0, 1, nan, 0, nan, nan, 0; - svd.compute(m, ComputeFullU | ComputeFullV); + svd.compute(m); VERIFY(svd.info() == InvalidInput); } @@ -355,8 +301,8 @@ void svd_underoverflow() Matrix2d M; M << -7.90884e-313, -4.94e-324, 0, 5.60844e-313; - SVD_DEFAULT(Matrix2d) svd; - svd.compute(M,ComputeFullU|ComputeFullV); + SVD_STATIC_OPTIONS(Matrix2d, ComputeFullU | ComputeFullV) svd; + svd.compute(M); CALL_SUBTEST( svd_check_full(M,svd) ); // Check all 2x2 matrices made with the following coefficients: @@ -367,7 +313,7 @@ void svd_underoverflow() do { M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3)); - svd.compute(M,ComputeFullU|ComputeFullV); + svd.compute(M); CALL_SUBTEST( svd_check_full(M,svd) ); id(k)++; @@ -390,16 +336,13 @@ void svd_underoverflow() 3.7841695601406358e+307, 2.4331702789740617e+306, -3.5235707140272905e+307, -8.7190887618028355e+307, -7.3453213709232193e+307, -2.4367363684472105e+307; - SVD_DEFAULT(Matrix3d) svd3; - svd3.compute(M3,ComputeFullU|ComputeFullV); // just check we don't loop indefinitely + SVD_STATIC_OPTIONS(Matrix3d, ComputeFullU | ComputeFullV) svd3; + svd3.compute(M3); // just check we don't loop indefinitely CALL_SUBTEST( svd_check_full(M3,svd3) ); } -// void jacobisvd(const MatrixType& a = MatrixType(), bool pickrandom = true) - -template -void svd_all_trivial_2x2( void (*cb)(const MatrixType&,bool) ) -{ +template +void svd_all_trivial_2x2(void (*cb)(const MatrixType&)) { MatrixType M; VectorXd value_set(3); value_set << 0, 1, -1; @@ -408,9 +351,9 @@ void svd_all_trivial_2x2( void (*cb)(const MatrixType&,bool) ) do { M << value_set(id(0)), value_set(id(1)), value_set(id(2)), value_set(id(3)); - - cb(M,false); - + + cb(M); + id(k)++; if(id(k)>=value_set.size()) { @@ -434,22 +377,10 @@ void svd_preallocate() internal::set_is_malloc_allowed(true); svd.compute(m); VERIFY_IS_APPROX(svd.singularValues(), v); + VERIFY_RAISES_ASSERT(svd.matrixU()); + VERIFY_RAISES_ASSERT(svd.matrixV()); - SVD_DEFAULT(MatrixXf) svd2(3,3); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - VERIFY_IS_APPROX(svd2.singularValues(), v); - VERIFY_RAISES_ASSERT(svd2.matrixU()); - VERIFY_RAISES_ASSERT(svd2.matrixV()); - svd2.compute(m, ComputeFullU | ComputeFullV); - VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); - VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); - internal::set_is_malloc_allowed(false); - svd2.compute(m); - internal::set_is_malloc_allowed(true); - - SVD_DEFAULT(MatrixXf) svd3(3,3,ComputeFullU|ComputeFullV); + SVD_STATIC_OPTIONS(MatrixXf, ComputeFullU | ComputeFullV) svd2(3, 3); internal::set_is_malloc_allowed(false); svd2.compute(m); internal::set_is_malloc_allowed(true); @@ -457,13 +388,201 @@ void svd_preallocate() VERIFY_IS_APPROX(svd2.matrixU(), Matrix3f::Identity()); VERIFY_IS_APPROX(svd2.matrixV(), Matrix3f::Identity()); internal::set_is_malloc_allowed(false); - svd2.compute(m, ComputeFullU|ComputeFullV); + svd2.compute(m); internal::set_is_malloc_allowed(true); } -template -void svd_verify_assert(const MatrixType& m, bool fullOnly = false) -{ +template +void svd_verify_assert_full_only(const MatrixType& m = MatrixType()) { + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime }; + + typedef Matrix RhsType; + RhsType rhs = RhsType::Zero(m.rows()); + + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner) svd0; + VERIFY_RAISES_ASSERT((svd0.matrixU())); + VERIFY_RAISES_ASSERT((svd0.singularValues())); + VERIFY_RAISES_ASSERT((svd0.matrixV())); + VERIFY_RAISES_ASSERT((svd0.solve(rhs))); + VERIFY_RAISES_ASSERT((svd0.transpose().solve(rhs))); + VERIFY_RAISES_ASSERT((svd0.adjoint().solve(rhs))); + + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner) svd1(m); + VERIFY_RAISES_ASSERT((svd1.matrixU())); + VERIFY_RAISES_ASSERT((svd1.matrixV())); + VERIFY_RAISES_ASSERT((svd1.solve(rhs))); + + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeFullU) svdFullU(m); + VERIFY_RAISES_ASSERT((svdFullU.matrixV())); + VERIFY_RAISES_ASSERT((svdFullU.solve(rhs))); + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeFullV) svdFullV(m); + VERIFY_RAISES_ASSERT((svdFullV.matrixU())); + VERIFY_RAISES_ASSERT((svdFullV.solve(rhs))); +} + +template +void svd_verify_assert(const MatrixType& m = MatrixType()) { + enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime }; + + typedef Matrix RhsType; + RhsType rhs = RhsType::Zero(m.rows()); + + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeThinU) svdThinU(m); + VERIFY_RAISES_ASSERT((svdThinU.matrixV())); + VERIFY_RAISES_ASSERT((svdThinU.solve(rhs))); + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeThinV) svdThinV(m); + VERIFY_RAISES_ASSERT((svdThinV.matrixU())); + VERIFY_RAISES_ASSERT((svdThinV.solve(rhs))); + + svd_verify_assert_full_only(m); +} + +template +void svd_compute_checks(const MatrixType& m) { + typedef SVD_STATIC_OPTIONS(MatrixType, Options) SVDType; + + enum { + RowsAtCompileTime = MatrixType::RowsAtCompileTime, + ColsAtCompileTime = MatrixType::ColsAtCompileTime, + DiagAtCompileTime = internal::min_size_prefer_dynamic(RowsAtCompileTime, ColsAtCompileTime), + MatrixURowsAtCompileTime = SVDType::MatrixUType::RowsAtCompileTime, + MatrixUColsAtCompileTime = SVDType::MatrixUType::ColsAtCompileTime, + MatrixVRowsAtCompileTime = SVDType::MatrixVType::RowsAtCompileTime, + MatrixVColsAtCompileTime = SVDType::MatrixVType::ColsAtCompileTime + }; + + SVDType staticSvd(m); + + VERIFY(MatrixURowsAtCompileTime == RowsAtCompileTime); + VERIFY(MatrixVRowsAtCompileTime == ColsAtCompileTime); + if (Options & ComputeThinU) VERIFY(MatrixUColsAtCompileTime == DiagAtCompileTime); + if (Options & ComputeFullU) VERIFY(MatrixUColsAtCompileTime == RowsAtCompileTime); + if (Options & ComputeThinV) VERIFY(MatrixVColsAtCompileTime == DiagAtCompileTime); + if (Options & ComputeFullV) VERIFY(MatrixVColsAtCompileTime == ColsAtCompileTime); + + if (Options & (ComputeThinU | ComputeFullU)) + VERIFY(staticSvd.computeU()); + else + VERIFY(!staticSvd.computeU()); + if (Options & (ComputeThinV | ComputeFullV)) + VERIFY(staticSvd.computeV()); + else + VERIFY(!staticSvd.computeV()); + + if (staticSvd.computeU()) VERIFY(staticSvd.matrixU().isUnitary()); + if (staticSvd.computeV()) VERIFY(staticSvd.matrixV().isUnitary()); + + if (staticSvd.computeU() && staticSvd.computeV()) { + svd_test_solvers(m, staticSvd); + svd_least_square(m); + // svd_min_norm generates non-square matrices so it can't be used with NoQRPreconditioner + if ((Options & internal::QRPreconditionerBits) != NoQRPreconditioner) svd_min_norm(m); + } +} + +template +void svd_check_constructor_options(const MatrixType& m, unsigned int computationOptions) { + const bool thinUnitary = (computationOptions & ComputeThinU) != 0 || (computationOptions & ComputeThinV) != 0; + + if (SvdType::ColsAtCompileTime != Dynamic && thinUnitary) { + VERIFY_RAISES_ASSERT(SvdType svd(m, computationOptions)); + return; + } + + Index diagSize = (std::min)(m.rows(), m.cols()); + + SvdType svd(m, computationOptions); + if (svd.computeU()) { + VERIFY(svd.matrixU().isUnitary()); + if (computationOptions & ComputeThinU) VERIFY(svd.matrixU().cols() == diagSize); + } + + if (svd.computeV()) { + VERIFY(svd.matrixV().isUnitary()); + if (computationOptions & ComputeThinV) VERIFY(svd.matrixV().cols() == diagSize); + } + if (svd.computeU() && svd.computeV()) { + svd_test_solvers(m, svd); + svd.matrixU().isUnitary(); + svd.matrixV().isUnitary(); + } +} + +template +void svd_option_checks(const MatrixType& m) { + svd_compute_checks(m); + svd_compute_checks(m); + svd_compute_checks(m); + svd_compute_checks(m); + + svd_compute_checks(m); + svd_compute_checks(m); + svd_compute_checks(m); + + svd_compute_checks(m); + svd_compute_checks(m); + + typedef SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeFullU | ComputeFullV) FullSvdType; + FullSvdType fullSvd(m); + svd_check_full(m, fullSvd); + svd_compare_to_full(m, fullSvd); + + typedef SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner) DynamicSvd; + svd_check_constructor_options(m, 0); + svd_check_constructor_options(m, ComputeThinU); + svd_check_constructor_options(m, ComputeThinV); + svd_check_constructor_options(m, ComputeThinU | ComputeThinV); + + svd_check_constructor_options(m, ComputeFullU); + svd_check_constructor_options(m, ComputeFullV); + svd_check_constructor_options(m, ComputeFullU | ComputeFullV); + + svd_check_constructor_options(m, ComputeThinU | ComputeFullV); + svd_check_constructor_options(m, ComputeFullU | ComputeThinV); +} + +template +void svd_option_checks_full_only(const MatrixType& m) { + svd_compute_checks(m); + svd_compute_checks(m); + svd_compute_checks(m); + + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeFullU | ComputeFullV) fullSvd(m); + svd_check_full(m, fullSvd); +} + +template +void svd_check_max_size_matrix(int initialRows, int initialCols) { + enum { + MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + }; + + int rows = MaxRowsAtCompileTime == Dynamic ? initialRows : (std::min)(initialRows, (int)MaxRowsAtCompileTime); + int cols = MaxColsAtCompileTime == Dynamic ? initialCols : (std::min)(initialCols, (int)MaxColsAtCompileTime); + + MatrixType m(rows, cols); + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeThinU | ComputeThinV) thinSvd(m); + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeThinU | ComputeFullV) mixedSvd1(m); + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeFullU | ComputeThinV) mixedSvd2(m); + SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeFullU | ComputeFullV) fullSvd(m); + + MatrixType n(MaxRowsAtCompileTime, MaxColsAtCompileTime); + thinSvd.compute(n); + mixedSvd1.compute(n); + mixedSvd2.compute(n); + fullSvd.compute(n); + + MatrixX dynamicMatrix(MaxRowsAtCompileTime + 1, MaxColsAtCompileTime + 1); + + VERIFY_RAISES_ASSERT(thinSvd.compute(dynamicMatrix)); + VERIFY_RAISES_ASSERT(mixedSvd1.compute(dynamicMatrix)); + VERIFY_RAISES_ASSERT(mixedSvd2.compute(dynamicMatrix)); + VERIFY_RAISES_ASSERT(fullSvd.compute(dynamicMatrix)); +} + +template +void svd_verify_constructor_options_assert(const MatrixType& m, bool fullOnly = false) { typedef typename MatrixType::Scalar Scalar; Index rows = m.rows(); Index cols = m.cols(); @@ -482,40 +601,38 @@ void svd_verify_assert(const MatrixType& m, bool fullOnly = false) VERIFY_RAISES_ASSERT(svd.solve(rhs)) VERIFY_RAISES_ASSERT(svd.transpose().solve(rhs)) VERIFY_RAISES_ASSERT(svd.adjoint().solve(rhs)) - MatrixType a = MatrixType::Zero(rows, cols); - a.setZero(); - svd.compute(a, 0); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.matrixV()) - svd.singularValues(); - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - svd.compute(a, ComputeFullU); - svd.matrixU(); - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - svd.compute(a, ComputeFullV); - svd.matrixV(); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) + MatrixType a = MatrixType::Zero(rows, cols); + SvdType svd2(a, 0); + VERIFY_RAISES_ASSERT(svd2.matrixU()) + VERIFY_RAISES_ASSERT(svd2.matrixV()) + svd2.singularValues(); + VERIFY_RAISES_ASSERT(svd2.solve(rhs)) + + SvdType svd3(a, ComputeFullU); + svd3.matrixU(); + VERIFY_RAISES_ASSERT(svd3.matrixV()) + VERIFY_RAISES_ASSERT(svd3.solve(rhs)) + + SvdType svd4(a, ComputeFullV); + svd4.matrixV(); + VERIFY_RAISES_ASSERT(svd4.matrixU()) + VERIFY_RAISES_ASSERT(svd4.solve(rhs)) if (!fullOnly && ColsAtCompileTime == Dynamic) { - svd.compute(a, ComputeThinU); - svd.matrixU(); - VERIFY_RAISES_ASSERT(svd.matrixV()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - svd.compute(a, ComputeThinV); - svd.matrixV(); - VERIFY_RAISES_ASSERT(svd.matrixU()) - VERIFY_RAISES_ASSERT(svd.solve(rhs)) - } - else - { - VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinU)) - VERIFY_RAISES_ASSERT(svd.compute(a, ComputeThinV)) + SvdType svd5(a, ComputeThinU); + svd5.matrixU(); + VERIFY_RAISES_ASSERT(svd5.matrixV()) + VERIFY_RAISES_ASSERT(svd5.solve(rhs)) + + SvdType svd6(a, ComputeThinV); + svd6.matrixV(); + VERIFY_RAISES_ASSERT(svd6.matrixU()) + VERIFY_RAISES_ASSERT(svd6.solve(rhs)) } } #undef SVD_DEFAULT #undef SVD_FOR_MIN_NORM +#undef SVD_STATIC_OPTIONS