mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-06 19:10:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Update SVD Module with Options template parameter

Browse Source
This commit is contained in:
Arthur 2022-02-02 00:15:44 +00:00 committed by Rasmus Munk Larsen
parent 89c6ab2385
commit 18b50458b6
24 changed files with 1112 additions and 815 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Eigen/src/Core/MatrixBase.h
View File

@ -375,8 +375,15 @@ template<typename Derived> class MatrixBase
/////////// SVD module ///////////
inline JacobiSVD<PlainObject> jacobiSvd(unsigned int computationOptions = 0) const;
inline BDCSVD<PlainObject> bdcSvd(unsigned int computationOptions = 0) const;
template<int Options = 0>
inline JacobiSVD<PlainObject, Options> jacobiSvd() const;
template<int Options = 0>
inline JacobiSVD<PlainObject, Options> jacobiSvd(unsigned int computationOptions) const;
template<int Options = 0>
inline BDCSVD<PlainObject, Options> bdcSvd() const;
template<int Options = 0>
inline BDCSVD<PlainObject, Options> bdcSvd(unsigned int computationOptions) const;
/////////// Geometry module ///////////

12
Eigen/src/Core/util/Constants.h
View File

@ -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

4
Eigen/src/Core/util/ForwardDeclarations.h
View File

@ -260,8 +260,8 @@ template<typename MatrixType> class ColPivHouseholderQR;
template<typename MatrixType> class FullPivHouseholderQR;
template<typename MatrixType> class CompleteOrthogonalDecomposition;
template<typename MatrixType> class SVDBase;
template<typename MatrixType, int QRPreconditioner = ColPivHouseholderQRPreconditioner> class JacobiSVD;
template<typename MatrixType> class BDCSVD;
template<typename MatrixType, int Options = 0> class JacobiSVD;
template<typename MatrixType, int Options = 0> class BDCSVD;
template<typename MatrixType, int UpLo = Lower> class LLT;
template<typename MatrixType, int UpLo = Lower> class LDLT;
template<typename VectorsType, typename CoeffsType, int Side=OnTheLeft> class HouseholderSequence;

2
Eigen/src/Geometry/Hyperplane.h
View File

@ -108,7 +108,7 @@ public:
if(norm <= v0.norm() * v1.norm() * NumTraits<RealScalar>::epsilon())
{
Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
JacobiSVD<Matrix<Scalar,2,3>, ComputeFullV> svd(m);
result.normal() = svd.matrixV().col(2);
}
else

2
Eigen/src/Geometry/Quaternion.h
View File

@ -651,7 +651,7 @@ EIGEN_DEVICE_FUNC inline Derived& QuaternionBase<Derived>::setFromTwoVectors(con
{
c = numext::maxi(c,Scalar(-1));
Matrix<Scalar,2,3> m; m << v0.transpose(), v1.transpose();
JacobiSVD<Matrix<Scalar,2,3> > svd(m, ComputeFullV);
JacobiSVD<Matrix<Scalar,2,3>, ComputeFullV> svd(m);
Vector3 axis = svd.matrixV().col(2);
Scalar w2 = (Scalar(1)+c)*Scalar(0.5);

4
Eigen/src/Geometry/Transform.h
View File

@ -1105,7 +1105,7 @@ template<typename RotationMatrixType, typename ScalingMatrixType>
EIGEN_DEVICE_FUNC void Transform<Scalar,Dim,Mode,Options>::computeRotationScaling(RotationMatrixType *rotation, ScalingMatrixType *scaling) const
{
// Note that JacobiSVD is faster than BDCSVD for small matrices.
JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU | ComputeFullV);
JacobiSVD<LinearMatrixType, ComputeFullU | ComputeFullV> 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<typename ScalingMatrixType, typename RotationMatrixType>
EIGEN_DEVICE_FUNC void Transform<Scalar,Dim,Mode,Options>::computeScalingRotation(ScalingMatrixType *scaling, RotationMatrixType *rotation) const
{
// Note that JacobiSVD is faster than BDCSVD for small matrices.
JacobiSVD<LinearMatrixType> svd(linear(), ComputeFullU | ComputeFullV);
JacobiSVD<LinearMatrixType, ComputeFullU | ComputeFullV> 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());

2
Eigen/src/Geometry/Umeyama.h
View File

@ -127,7 +127,7 @@ umeyama(const MatrixBase<Derived>& src, const MatrixBase<OtherDerived>& dst, boo
// Eq. (38)
const MatrixType sigma = one_over_n * dst_demean * src_demean.transpose();
JacobiSVD<MatrixType> svd(sigma, ComputeFullU | ComputeFullV);
JacobiSVD<MatrixType, ComputeFullU | ComputeFullV> svd(sigma);
// Initialize the resulting transformation with an identity matrix...
TransformationMatrixType Rt = TransformationMatrixType::Identity(m+1,m+1);

320
Eigen/src/SVD/BDCSVD.h
View File

@ -39,19 +39,32 @@ namespace Eigen {
IOFormat bdcsvdfmt(8, 0, ", ", "\n", " [", "]");
#endif
template<typename MatrixType_> class BDCSVD;
template <typename MatrixType_, int Options>
class BDCSVD;
namespace internal {
template<typename MatrixType_>
struct traits<BDCSVD<MatrixType_> >
: traits<MatrixType_>
{
template <typename MatrixType_, int Options>
struct traits<BDCSVD<MatrixType_, Options> > : svd_traits<MatrixType_, Options> {
typedef MatrixType_ MatrixType;
};
} // end namespace internal
template <typename MatrixType, int Options>
struct allocate_small_svd {
static void run(JacobiSVD<MatrixType, Options>& smallSvd, Index rows, Index cols, unsigned int computationOptions) {
(void)computationOptions;
smallSvd = JacobiSVD<MatrixType, Options>(rows, cols);
}
};
template <typename MatrixType>
struct allocate_small_svd<MatrixType, 0> {
static void run(JacobiSVD<MatrixType>& smallSvd, Index rows, Index cols, unsigned int computationOptions) {
smallSvd = JacobiSVD<MatrixType>(rows, cols, computationOptions);
}
};
} // end namespace internal
/** \ingroup SVD_Module
*
@ -62,6 +75,11 @@ struct traits<BDCSVD<MatrixType_> >
*
* \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<BDCSVD<MatrixType_> >
*
* \sa class JacobiSVD
*/
template<typename MatrixType_>
class BDCSVD : public SVDBase<BDCSVD<MatrixType_> >
{
template <typename MatrixType_, int Options>
class BDCSVD : public SVDBase<BDCSVD<MatrixType_, Options> > {
typedef SVDBase<BDCSVD> 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.
*
* <b>Note: This constructor is deprecated.</b>
* 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<MatrixType, Options>(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.
*
* <b>Note: This constructor is deprecated.</b>
* 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<MatrixType, Options>(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<MatrixXr,Dynamic,Dynamic> 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 <typename SVDType>
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<MatrixType, Options> 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<typename MatrixType>
void BDCSVD<MatrixType>::allocate(Eigen::Index rows, Eigen::Index cols, unsigned int computationOptions)
{
m_isTranspose = (cols > rows);
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::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<MatrixType, Options>::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<MatrixType>::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<typename MatrixType>
BDCSVD<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsigned int computationOptions)
{
template <typename MatrixType, int Options>
BDCSVD<MatrixType, Options>& BDCSVD<MatrixType, Options>::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<RealScalar>::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<MatrixType> 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<MatrixType>& BDCSVD<MatrixType>::compute(const MatrixType& matrix, unsign
m_isInitialized = true;
return *this;
}// end compute
} // end compute
template<typename MatrixType>
template<typename HouseholderU, typename HouseholderV, typename NaiveU, typename NaiveV>
void BDCSVD<MatrixType>::copyUV(const HouseholderU &householderU, const HouseholderV &householderV, const NaiveU &naiveU, const NaiveV &naiveV)
{
template <typename MatrixType, int Options>
template <typename HouseholderU, typename HouseholderV, typename NaiveU, typename NaiveV>
void BDCSVD<MatrixType, Options>::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<MatrixType>::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<typename MatrixType>
void BDCSVD<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> A, const MatrixXr &B, Index n1)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::structured_update(Block<MatrixXr, Dynamic, Dynamic> A, const MatrixXr& B, Index n1) {
Index n = A.rows();
if(n>100)
{
@ -403,8 +425,26 @@ void BDCSVD<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> 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 <typename MatrixType, int Options>
template <typename SVDType>
void BDCSVD<MatrixType, Options>::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<MatrixType>::structured_update(Block<MatrixXr,Dynamic,Dynamic> 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<typename MatrixType>
void BDCSVD<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::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<MatrixType>::divide(Eigen::Index firstCol, Eigen::Index lastCol, Eig
// matrices.
if (n < m_algoswap)
{
// FIXME this line involves temporaries
JacobiSVD<MatrixXr> b(m_computed.block(firstCol, firstCol, n + 1, n), ComputeFullU | (m_compV ? ComputeFullV : 0));
m_info = b.info();
if (m_info != Success && m_info != NoConvergence) return;
if (m_compU)
m_naiveU.block(firstCol, firstCol, n + 1, n + 1).real() = b.matrixU();
else
{
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<MatrixXr, ComputeFullU | ComputeFullV> baseSvd;
computeBaseCase(baseSvd, n, firstCol, firstRowW, firstColW, shift);
} else {
JacobiSVD<MatrixXr, ComputeFullU> 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<MatrixType>::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<MatrixType>::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 <typename MatrixType>
void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, MatrixXr& U, VectorType& singVals, MatrixXr& V)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, MatrixXr& U,
VectorType& singVals, MatrixXr& V) {
const RealScalar considerZero = (std::numeric_limits<RealScalar>::min)();
using std::abs;
ArrayRef col0 = m_computed.col(firstCol).segment(firstCol, n);
@ -728,9 +762,10 @@ void BDCSVD<MatrixType>::computeSVDofM(Eigen::Index firstCol, Eigen::Index n, Ma
#endif
}
template <typename MatrixType>
typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar mu, const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const ArrayRef& diagShifted, RealScalar shift)
{
template <typename MatrixType, int Options>
typename BDCSVD<MatrixType, Options>::RealScalar BDCSVD<MatrixType, Options>::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<m; ++i)
@ -741,13 +776,11 @@ typename BDCSVD<MatrixType>::RealScalar BDCSVD<MatrixType>::secularEq(RealScalar
res += (col0(j) / (diagShifted(j) - mu)) * (col0(j) / (diag(j) + shift + mu));
}
return res;
}
template <typename MatrixType>
void BDCSVD<MatrixType>::computeSingVals(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm,
VectorType& singVals, ArrayRef shifts, ArrayRef mus)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::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<MatrixType>::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 <typename MatrixType>
void BDCSVD<MatrixType>::perturbCol0
(const ArrayRef& col0, const ArrayRef& diag, const IndicesRef &perm, const VectorType& singVals,
const ArrayRef& shifts, const ArrayRef& mus, ArrayRef zhat)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::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<MatrixType>::perturbCol0
}
// compute singular vectors
template <typename MatrixType>
void BDCSVD<MatrixType>::computeSingVecs
(const ArrayRef& zhat, const ArrayRef& diag, const IndicesRef &perm, const VectorType& singVals,
const ArrayRef& shifts, const ArrayRef& mus, MatrixXr& U, MatrixXr& V)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::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<MatrixType>::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 <typename MatrixType>
void BDCSVD<MatrixType>::deflation43(Eigen::Index firstCol, Eigen::Index shift, Eigen::Index i, Eigen::Index size)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::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<MatrixType>::deflation43(Eigen::Index firstCol, Eigen::Index shift,
JacobiRotation<RealScalar> 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 <typename MatrixType>
void BDCSVD<MatrixType>::deflation44(Eigen::Index firstColu , Eigen::Index firstColm, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index i, Eigen::Index j, Eigen::Index size)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::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<MatrixType>::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 <typename MatrixType>
void BDCSVD<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol, Eigen::Index k, Eigen::Index firstRowW, Eigen::Index firstColW, Eigen::Index shift)
{
template <typename MatrixType, int Options>
void BDCSVD<MatrixType, Options>::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<MatrixType>::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<MatrixType>::deflation(Eigen::Index firstCol, Eigen::Index lastCol,
*
* \sa class BDCSVD
*/
template<typename Derived>
BDCSVD<typename MatrixBase<Derived>::PlainObject>
MatrixBase<Derived>::bdcSvd(unsigned int computationOptions) const
{
return BDCSVD<PlainObject>(*this, computationOptions);
template <typename Derived>
template <int Options>
BDCSVD<typename MatrixBase<Derived>::PlainObject, Options> MatrixBase<Derived>::bdcSvd() const {
return BDCSVD<PlainObject, Options>(*this);
}
/** \svd_module
*
* \return the singular value decomposition of \c *this computed by Divide & Conquer algorithm
*
* \sa class BDCSVD
*/
template <typename Derived>
template <int Options>
BDCSVD<typename MatrixBase<Derived>::PlainObject, Options> MatrixBase<Derived>::bdcSvd(
unsigned int computationOptions) const {
return BDCSVD<PlainObject, Options>(*this, computationOptions);
}
} // end namespace Eigen

647
Eigen/src/SVD/JacobiSVD.h
View File

@ -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<typename MatrixType, int QRPreconditioner,
bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
template <typename MatrixType, int Options, bool IsComplex = NumTraits<typename MatrixType::Scalar>::IsComplex>
struct svd_precondition_2x2_block_to_be_real {};
/*** QR preconditioners (R-SVD)
@ -46,37 +46,31 @@ struct qr_preconditioner_should_do_anything
};
};
template<typename MatrixType, int QRPreconditioner, int Case,
bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret
> struct qr_preconditioner_impl {};
template <typename MatrixType, int Options, int QRPreconditioner, int Case,
bool DoAnything = qr_preconditioner_should_do_anything<MatrixType, QRPreconditioner, Case>::ret>
struct qr_preconditioner_impl {};
template<typename MatrixType, int QRPreconditioner, int Case>
class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false>
{
public:
void allocate(const JacobiSVD<MatrixType, QRPreconditioner>&) {}
bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&)
{
return false;
}
template <typename MatrixType, int Options, int QRPreconditioner, int Case>
class qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, Case, false> {
public:
void allocate(const JacobiSVD<MatrixType, Options>&) {}
bool run(JacobiSVD<MatrixType, Options>&, const MatrixType&) { return false; }
};
/*** preconditioner using FullPivHouseholderQR ***/
template<typename MatrixType>
class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
{
public:
template <typename MatrixType, int Options>
class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols,
true> {
public:
typedef typename MatrixType::Scalar Scalar;
enum
{
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime
};
typedef Matrix<Scalar, 1, RowsAtCompileTime, RowMajor, 1, MaxRowsAtCompileTime> WorkspaceType;
typedef JacobiSVD<MatrixType, Options> SVDType;
void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd)
{
enum { WorkspaceSize = MatrixType::RowsAtCompileTime, MaxWorkspaceSize = MatrixType::MaxRowsAtCompileTime };
typedef Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize> 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<MatrixType, FullPivHouseholderQRPreconditioner>& 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<MatrixType> QRType;
QRType m_qr;
WorkspaceType m_workspace;
};
template<typename MatrixType>
class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
{
public:
template <typename MatrixType, int Options>
class qr_preconditioner_impl<MatrixType, Options, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows,
true> {
public:
typedef typename MatrixType::Scalar Scalar;
enum
{
typedef JacobiSVD<MatrixType, Options> 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<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions,
MaxColsAtCompileTime, MaxRowsAtCompileTime>::type
TransposeTypeWithSameStorageOrder;
void allocate(const JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& 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<MatrixType, FullPivHouseholderQRPreconditioner>& 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<TransposeTypeWithSameStorageOrder> QRType;
QRType m_qr;
TransposeTypeWithSameStorageOrder m_adjoint;
typename internal::plain_row_type<MatrixType>::type m_workspace;
typename plain_row_type<MatrixType>::type m_workspace;
};
/*** preconditioner using ColPivHouseholderQR ***/
template<typename MatrixType>
class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
{
public:
void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
{
template <typename MatrixType, int Options>
class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols,
true> {
public:
typedef typename MatrixType::Scalar Scalar;
typedef JacobiSVD<MatrixType, Options> SVDType;
enum {
WorkspaceSize = internal::traits<SVDType>::MatrixUColsAtCompileTime,
MaxWorkspaceSize = internal::traits<SVDType>::MatrixUMaxColsAtCompileTime
};
typedef Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize> 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<MatrixType, ColPivHouseholderQRPreconditioner>& 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<MatrixType> QRType;
QRType m_qr;
typename internal::plain_col_type<MatrixType>::type m_workspace;
WorkspaceType m_workspace;
};
template<typename MatrixType>
class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
{
public:
template <typename MatrixType, int Options>
class qr_preconditioner_impl<MatrixType, Options, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows,
true> {
public:
typedef typename MatrixType::Scalar Scalar;
enum
{
typedef JacobiSVD<MatrixType, Options> SVDType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
Options = MatrixType::Options
MatrixOptions = MatrixType::Options,
WorkspaceSize = internal::traits<SVDType>::MatrixVColsAtCompileTime,
MaxWorkspaceSize = internal::traits<SVDType>::MatrixVMaxColsAtCompileTime
};
typedef typename internal::make_proper_matrix_type<
Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
>::type TransposeTypeWithSameStorageOrder;
typedef Matrix<Scalar, WorkspaceSize, 1, ColMajor, MaxWorkspaceSize, 1> WorkspaceType;
void allocate(const JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd)
{
typedef typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions,
MaxColsAtCompileTime, MaxRowsAtCompileTime>::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<MatrixType, ColPivHouseholderQRPreconditioner>& 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<TransposeTypeWithSameStorageOrder> QRType;
QRType m_qr;
TransposeTypeWithSameStorageOrder m_adjoint;
typename internal::plain_row_type<MatrixType>::type m_workspace;
WorkspaceType m_workspace;
};
/*** preconditioner using HouseholderQR ***/
template<typename MatrixType>
class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true>
{
public:
void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
{
template <typename MatrixType, int Options>
class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> {
public:
typedef typename MatrixType::Scalar Scalar;
typedef JacobiSVD<MatrixType, Options> SVDType;
enum {
WorkspaceSize = internal::traits<SVDType>::MatrixUColsAtCompileTime,
MaxWorkspaceSize = internal::traits<SVDType>::MatrixUMaxColsAtCompileTime
};
typedef Matrix<Scalar, 1, WorkspaceSize, RowMajor, 1, MaxWorkspaceSize> 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<MatrixType, HouseholderQRPreconditioner>& 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<MatrixType> QRType;
QRType m_qr;
typename internal::plain_col_type<MatrixType>::type m_workspace;
WorkspaceType m_workspace;
};
template<typename MatrixType>
class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true>
{
public:
template <typename MatrixType, int Options>
class qr_preconditioner_impl<MatrixType, Options, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> {
public:
typedef typename MatrixType::Scalar Scalar;
enum
{
typedef JacobiSVD<MatrixType, Options> SVDType;
enum {
RowsAtCompileTime = MatrixType::RowsAtCompileTime,
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime,
Options = MatrixType::Options
MatrixOptions = MatrixType::Options,
WorkspaceSize = internal::traits<SVDType>::MatrixVColsAtCompileTime,
MaxWorkspaceSize = internal::traits<SVDType>::MatrixVMaxColsAtCompileTime
};
typedef typename internal::make_proper_matrix_type<
Scalar, ColsAtCompileTime, RowsAtCompileTime, Options, MaxColsAtCompileTime, MaxRowsAtCompileTime
>::type TransposeTypeWithSameStorageOrder;
typedef Matrix<Scalar, WorkspaceSize, 1, ColMajor, MaxWorkspaceSize, 1> WorkspaceType;
void allocate(const JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd)
{
typedef typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, RowsAtCompileTime, MatrixOptions,
MaxColsAtCompileTime, MaxRowsAtCompileTime>::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<MatrixType, HouseholderQRPreconditioner>& 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<TransposeTypeWithSameStorageOrder> QRType;
QRType m_qr;
TransposeTypeWithSameStorageOrder m_adjoint;
typename internal::plain_row_type<MatrixType>::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<typename MatrixType, int QRPreconditioner>
struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, false>
{
typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
template <typename MatrixType, int Options>
struct svd_precondition_2x2_block_to_be_real<MatrixType, Options, false> {
typedef JacobiSVD<MatrixType, Options> SVD;
typedef typename MatrixType::RealScalar RealScalar;
static bool run(typename SVD::WorkMatrixType&, SVD&, Index, Index, RealScalar&) { return true; }
};
template<typename MatrixType, int QRPreconditioner>
struct svd_precondition_2x2_block_to_be_real<MatrixType, QRPreconditioner, true>
{
typedef JacobiSVD<MatrixType, QRPreconditioner> SVD;
template <typename MatrixType, int Options>
struct svd_precondition_2x2_block_to_be_real<MatrixType, Options, true> {
typedef JacobiSVD<MatrixType, Options> 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<MatrixType, QRPreconditioner, true>
}
};
template<typename MatrixType_, int QRPreconditioner>
struct traits<JacobiSVD<MatrixType_,QRPreconditioner> >
: traits<MatrixType_>
{
template <typename MatrixType_, int Options>
struct traits<JacobiSVD<MatrixType_, Options> > : svd_traits<MatrixType_, Options> {
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<typename MatrixType_, int QRPreconditioner> class JacobiSVD
: public SVDBase<JacobiSVD<MatrixType_,QRPreconditioner> >
{
typedef SVDBase<JacobiSVD> 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<MatrixType,
* ColPivHouseholderQRPreconditioner | ComputeThinU | ComputeFullV>
*
* \sa MatrixBase::jacobiSvd()
*/
template <typename MatrixType_, int Options>
class JacobiSVD : public SVDBase<JacobiSVD<MatrixType_, Options> > {
typedef SVDBase<JacobiSVD> Base;
typedef MatrixType_ MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<typename MatrixType::Scalar>::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<typename MatrixType::Scalar>::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<MatrixType>::type RowType;
typedef typename internal::plain_col_type<MatrixType>::type ColType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime,
MatrixOptions, MaxDiagSizeAtCompileTime, MaxDiagSizeAtCompileTime>
WorkMatrixType;
typedef typename Base::MatrixUType MatrixUType;
typedef typename Base::MatrixVType MatrixVType;
typedef typename Base::SingularValuesType SingularValuesType;
typedef Matrix<Scalar, DiagSizeAtCompileTime, DiagSizeAtCompileTime, MatrixOptions, MaxDiagSizeAtCompileTime,
MaxDiagSizeAtCompileTime>
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.
*
* <b>Note: This constructor is deprecated.</b>
* 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<MatrixType, Options>(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.
*
* <b>Note: This constructor is deprecated.</b>
* 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<MatrixType, Options>(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<typename MatrixType__, int QRPreconditioner_, bool IsComplex_>
friend struct internal::svd_precondition_2x2_block_to_be_real;
template<typename MatrixType__, int QRPreconditioner_, int Case_, bool DoAnything_>
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<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols;
internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows;
MatrixType m_scaledMatrix;
template <typename MatrixType__, int Options_, bool IsComplex_>
friend struct internal::svd_precondition_2x2_block_to_be_real;
template <typename MatrixType__, int Options_, int QRPreconditioner_, int Case_, bool DoAnything_>
friend struct internal::qr_preconditioner_impl;
internal::qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>
m_qr_precond_morecols;
internal::qr_preconditioner_impl<MatrixType, Options, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>
m_qr_precond_morerows;
WorkMatrixType m_workMatrix;
MatrixType m_scaledMatrix;
};
template<typename MatrixType, int QRPreconditioner>
void JacobiSVD<MatrixType, QRPreconditioner>::allocate(Eigen::Index rows, Eigen::Index cols, unsigned int computationOptions)
{
eigen_assert(rows >= 0 && cols >= 0);
template <typename MatrixType, int Options>
void JacobiSVD<MatrixType, Options>::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<typename MatrixType, int QRPreconditioner>
JacobiSVD<MatrixType, QRPreconditioner>&
JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsigned int computationOptions)
{
template <typename MatrixType, int Options>
JacobiSVD<MatrixType, Options>& JacobiSVD<MatrixType, Options>::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<MatrixType, QRPreconditioner>::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<MatrixType, QRPreconditioner>::run(m_workMatrix, *this, p, q, maxDiagEntry))
{
if (internal::svd_precondition_2x2_block_to_be_real<MatrixType, Options>::run(m_workMatrix, *this, p, q,
maxDiagEntry)) {
JacobiRotation<RealScalar> j_left, j_right;
internal::real_2x2_jacobi_svd(m_workMatrix, p, q, &j_left, &j_right);
@ -802,13 +801,19 @@ JacobiSVD<MatrixType, QRPreconditioner>::compute(const MatrixType& matrix, unsig
*
* \sa class JacobiSVD
*/
template<typename Derived>
JacobiSVD<typename MatrixBase<Derived>::PlainObject>
MatrixBase<Derived>::jacobiSvd(unsigned int computationOptions) const
{
return JacobiSVD<PlainObject>(*this, computationOptions);
template <typename Derived>
template <int Options>
JacobiSVD<typename MatrixBase<Derived>::PlainObject, Options> MatrixBase<Derived>::jacobiSvd() const {
return JacobiSVD<PlainObject, Options>(*this);
}
} // end namespace Eigen
template <typename Derived>
template <int Options>
JacobiSVD<typename MatrixBase<Derived>::PlainObject, Options> MatrixBase<Derived>::jacobiSvd(
unsigned int computationOptions) const {
return JacobiSVD<PlainObject, Options>(*this, computationOptions);
}
} // end namespace Eigen
#endif // EIGEN_JACOBISVD_H

41
Eigen/src/SVD/JacobiSVD_LAPACKE.h
View File

@ -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<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, ColPivHouseholderQRPreconditioner>& \
JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, ColPivHouseholderQRPreconditioner>::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>& matrix, unsigned int computationOptions) \
JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>& \
JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, OPTIONS>::compute(const Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>& matrix) \
{ \
typedef Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic> 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<Scalar>::epsilon();*/ \
m_nonzeroSingularValues = m_diagSize; \
@ -56,14 +56,14 @@ JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, 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<lapack_int>(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<lapack_int>(m_cols) : (m_computeThinV) ? internal::convert_index<lapack_int>(m_diagSize) : 1; \
lapack_int vt_rows = (ShouldComputeFullV) ? internal::convert_index<lapack_int>(m_cols) : (ShouldComputeThinV) ? internal::convert_index<lapack_int>(m_diagSize) : 1; \
if (computeV()) { \
localV.resize(vt_rows, m_cols); \
ldvt = internal::convert_index<lapack_int>(localV.outerStride()); \
@ -78,15 +78,26 @@ JacobiSVD<Matrix<EIGTYPE, Dynamic, Dynamic, EIGCOLROW, Dynamic, Dynamic>, 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

129
Eigen/src/SVD/SVDBase.h
View File

@ -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 <typename MatrixType, int Options>
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<typename Derived> struct traits<SVDBase<Derived> >
: traits<Derived>
{
@ -29,6 +59,32 @@ template<typename Derived> struct traits<SVDBase<Derived> >
typedef int StorageIndex;
enum { Flags = 0 };
};
template <typename MatrixType, int Options>
struct svd_traits : traits<MatrixType> {
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<typename MatrixType::Scalar>::Real RealScalar;
typedef typename Eigen::internal::traits<SVDBase>::StorageIndex StorageIndex;
typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
static constexpr bool ShouldComputeFullU = internal::traits<Derived>::ShouldComputeFullU;
static constexpr bool ShouldComputeThinU = internal::traits<Derived>::ShouldComputeThinU;
static constexpr bool ShouldComputeFullV = internal::traits<Derived>::ShouldComputeFullV;
static constexpr bool ShouldComputeThinV = internal::traits<Derived>::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<Derived>::MatrixUColsAtCompileTime,
MatrixVColsAtCompileTime = internal::traits<Derived>::MatrixVColsAtCompileTime,
MatrixUMaxColsAtCompileTime = internal::traits<Derived>::MatrixUMaxColsAtCompileTime,
MatrixVMaxColsAtCompileTime = internal::traits<Derived>::MatrixVMaxColsAtCompileTime
};
typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, MatrixOptions, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixUType;
typedef Matrix<Scalar, ColsAtCompileTime, ColsAtCompileTime, MatrixOptions, MaxColsAtCompileTime, MaxColsAtCompileTime> 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<Scalar, RowsAtCompileTime, MatrixUColsAtCompileTime, MatrixOptions,
MaxRowsAtCompileTime, MatrixUMaxColsAtCompileTime>::type MatrixUType;
typedef
typename internal::make_proper_matrix_type<Scalar, ColsAtCompileTime, MatrixVColsAtCompileTime, MatrixOptions,
MaxColsAtCompileTime, MatrixVMaxColsAtCompileTime>::type MatrixVType;
typedef typename internal::plain_diag_type<MatrixType, RealScalar>::type SingularValuesType;
Derived& derived() { return *static_cast<Derived*>(this); }
const Derived& derived() const { return *static_cast<const Derived*>(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<Derived>::_solve_impl_transposed(const RhsType &rhs, DstType &dst)
}
#endif
template<typename MatrixType>
bool SVDBase<MatrixType>::allocate(Index rows, Index cols, unsigned int computationOptions)
{
template <typename Derived>
bool SVDBase<Derived>::allocate(Index rows, Index cols, unsigned int computationOptions) {
eigen_assert(rows >= 0 && cols >= 0);
if (m_isAllocated &&
@ -349,14 +421,13 @@ bool SVDBase<MatrixType>::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);

10
bench/dense_solvers.cpp
View File

@ -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<Mat> cpqr(A.rows(),A.cols());
CompleteOrthogonalDecomposition<Mat> cod(A.rows(),A.cols());
FullPivHouseholderQR<Mat> fpqr(A.rows(),A.cols());
JacobiSVD<MatDyn> jsvd(A.rows(),A.cols());
BDCSVD<MatDyn> bdcsvd(A.rows(),A.cols());
JacobiSVD<MatDyn, ComputeThinU|ComputeThinV> jsvd(A.rows(),A.cols());
BDCSVD<MatDyn, ComputeThinU|ComputeThinV> 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();

4
doc/UsingBlasLapackBackends.dox
View File

@ -101,8 +101,8 @@ m1.colPivHouseholderQr();
?geqp3
\endcode</td></tr>
<tr class="alt"><td>Singular value decomposition \n \c EIGEN_USE_LAPACKE </td><td>\code
JacobiSVD<MatrixXd> svd;
svd.compute(m1, ComputeThinV);
JacobiSVD<MatrixXd, ComputeThinV> svd;
svd.compute(m1);
\endcode</td><td>\code
?gesvd
\endcode</td></tr>

12
doc/examples/TutorialLinAlgSVDSolve.cpp
View File

@ -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<ComputeThinU | ComputeThinV>().solve(b) << endl;
}

2
doc/snippets/JacobiSVD_basic.cpp
View File

@ -1,6 +1,6 @@
MatrixXf m = MatrixXf::Random(3,2);
cout << "Here is the matrix m:" << endl << m << endl;
JacobiSVD<MatrixXf> svd(m, ComputeThinU | ComputeThinV);
JacobiSVD<MatrixXf, ComputeThinU | ComputeThinV> 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;

2
lapack/svd.cpp
View File

@ -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;
}
}

111
test/bdcsvd.cpp
View File

@ -16,29 +16,12 @@
#include "main.h"
#include <Eigen/SVD>
#include <iostream>
#include <Eigen/LU>
#define SVD_DEFAULT(M) BDCSVD<M>
#define SVD_FOR_MIN_NORM(M) BDCSVD<M>
#define SVD_STATIC_OPTIONS(M, O) BDCSVD<M, O>
#include "svd_common.h"
// Check all variants of JacobiSVD
template<typename MatrixType>
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<BDCSVD<MatrixType> >(m, false) ));
}
template<typename MatrixType>
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<ComputeFullU | ComputeFullV>().solve(m), m);
VERIFY_IS_APPROX(m.template bdcSvd<ComputeFullU | ComputeFullV>().transpose().solve(m), m);
VERIFY_IS_APPROX(m.template bdcSvd<ComputeFullU | ComputeFullV>().adjoint().solve(m), m);
}
// Compare the Singular values returned with Jacobi and Bdc.
template<typename MatrixType>
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 <typename MatrixType>
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<MatrixType> bdc_svd(m.rows(), m.cols(), computationOptions);
BDCSVD<MatrixType> bdc_svd(m.rows(), m.cols());
bdc_svd.setSwitchSize(algoswap);
bdc_svd.compute(m);
JacobiSVD<MatrixType> 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 <typename MatrixType>
void bdcsvd_all_options(const MatrixType& input = MatrixType()) {
MatrixType m = input;
svd_fill_random(m);
svd_option_checks<MatrixType, 0>(m);
}
template <typename MatrixType>
void bdcsvd_verify_assert(const MatrixType& input = MatrixType()) {
svd_verify_assert<MatrixType>(input);
svd_verify_constructor_options_assert<BDCSVD<MatrixType>>(input);
}
EIGEN_DECLARE_TEST(bdcsvd)
{
CALL_SUBTEST_3(( svd_verify_assert<BDCSVD<Matrix3f> >(Matrix3f()) ));
CALL_SUBTEST_4(( svd_verify_assert<BDCSVD<Matrix4d> >(Matrix4d()) ));
CALL_SUBTEST_7(( svd_verify_assert<BDCSVD<MatrixXf> >(MatrixXf(10,12)) ));
CALL_SUBTEST_8(( svd_verify_assert<BDCSVD<MatrixXcd> >(MatrixXcd(7,5)) ));
CALL_SUBTEST_3((bdcsvd_verify_assert<Matrix3f>()));
CALL_SUBTEST_4((bdcsvd_verify_assert<Matrix4d>()));
CALL_SUBTEST_7((bdcsvd_verify_assert<Matrix<float, 30, 21>>()));
CALL_SUBTEST_7((bdcsvd_verify_assert<Matrix<float, 21, 30>>()));
CALL_SUBTEST_9((bdcsvd_verify_assert<Matrix<std::complex<double>, 20, 27>>()));
CALL_SUBTEST_101(( svd_all_trivial_2x2(bdcsvd<Matrix2cd>) ));
CALL_SUBTEST_102(( svd_all_trivial_2x2(bdcsvd<Matrix2d>) ));
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_3(( bdcsvd<Matrix3f>() ));
CALL_SUBTEST_4(( bdcsvd<Matrix4d>() ));
CALL_SUBTEST_5(( bdcsvd<Matrix<float,3,5> >() ));
CALL_SUBTEST_101((svd_all_trivial_2x2(bdcsvd_all_options<Matrix2cd>)));
CALL_SUBTEST_102((svd_all_trivial_2x2(bdcsvd_all_options<Matrix2d>)));
for (int i = 0; i < g_repeat; i++) {
int r = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2),
c = internal::random<int>(1, EIGEN_TEST_MAX_SIZE/2);
TEST_SET_BUT_UNUSED_VARIABLE(r)
TEST_SET_BUT_UNUSED_VARIABLE(c)
CALL_SUBTEST_6(( bdcsvd(Matrix<double,Dynamic,2>(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>(MatrixXf(r, c))));
CALL_SUBTEST_10((compare_bdc_jacobi<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_8((compare_bdc_jacobi<MatrixXcd>(MatrixXcd(r, c))));
// Test on inf/nan matrix
CALL_SUBTEST_7( (svd_inf_nan<BDCSVD<MatrixXf>, MatrixXf>()) );
CALL_SUBTEST_10( (svd_inf_nan<BDCSVD<MatrixXd>, MatrixXd>()) );
CALL_SUBTEST_7((svd_inf_nan<MatrixXf>()));
CALL_SUBTEST_10((svd_inf_nan<MatrixXd>()));
// Verify some computations using all combinations of the Options template parameter.
CALL_SUBTEST_3((bdcsvd_all_options<Matrix3f>()));
CALL_SUBTEST_3((bdcsvd_all_options<Matrix<float, 2, 3>>()));
CALL_SUBTEST_4((bdcsvd_all_options<Matrix<double, 20, 17>>()));
CALL_SUBTEST_4((bdcsvd_all_options<Matrix<double, 17, 20>>()));
CALL_SUBTEST_5((bdcsvd_all_options<Matrix<double, Dynamic, 30>>(Matrix<double, Dynamic, 30>(r, 30))));
CALL_SUBTEST_5((bdcsvd_all_options<Matrix<double, 20, Dynamic>>(Matrix<double, 20, Dynamic>(20, c))));
CALL_SUBTEST_7((bdcsvd_all_options<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_8((bdcsvd_all_options<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_10((bdcsvd_all_options<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_14((bdcsvd_all_options<Matrix<double, 20, 27, RowMajor>>()));
CALL_SUBTEST_14((bdcsvd_all_options<Matrix<double, 27, 20, RowMajor>>()));
CALL_SUBTEST_15((
svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, ColMajor, 20, 35>, ColPivHouseholderQRPreconditioner>(
r, c)));
CALL_SUBTEST_15(
(svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, ColMajor, 35, 20>, HouseholderQRPreconditioner>(r,
c)));
CALL_SUBTEST_15((
svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, RowMajor, 20, 35>, ColPivHouseholderQRPreconditioner>(
r, c)));
CALL_SUBTEST_15(
(svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, RowMajor, 35, 20>, HouseholderQRPreconditioner>(r,
c)));
}
// test matrixbase method

4
test/boostmultiprec.cpp
View File

@ -200,8 +200,8 @@ EIGEN_DECLARE_TEST(boostmultiprec)
TEST_SET_BUT_UNUSED_VARIABLE(s)
}
CALL_SUBTEST_9(( jacobisvd(Mat(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
CALL_SUBTEST_10(( bdcsvd(Mat(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
CALL_SUBTEST_9(( jacobisvd_all_options(Mat(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
CALL_SUBTEST_10(( bdcsvd_all_options(Mat(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
CALL_SUBTEST_11(( test_simplicial_cholesky_T<Real,int,ColMajor>() ));
}

6
test/diagonalmatrices.cpp
View File

@ -96,9 +96,9 @@ template<typename MatrixType> 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

2
test/geo_alignedbox.cpp
View File

@ -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<MatrixType> svd = q.jacobiSvd(ComputeFullU | ComputeFullV);
const JacobiSVD<MatrixType, ComputeFullU | ComputeFullV> 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;

173
test/jacobisvd.cpp
View File

@ -16,49 +16,9 @@
#define SVD_DEFAULT(M) JacobiSVD<M>
#define SVD_FOR_MIN_NORM(M) JacobiSVD<M,ColPivHouseholderQRPreconditioner>
#define SVD_STATIC_OPTIONS(M, O) JacobiSVD<M, O>
#include "svd_common.h"
// Check all variants of JacobiSVD
template<typename MatrixType>
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<JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> >(m, true) )); // check full only
CALL_SUBTEST(( svd_test_all_computation_options<JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner> >(m, false) ));
CALL_SUBTEST(( svd_test_all_computation_options<JacobiSVD<MatrixType, HouseholderQRPreconditioner> >(m, false) ));
if(m.rows()==m.cols())
CALL_SUBTEST(( svd_test_all_computation_options<JacobiSVD<MatrixType, NoQRPreconditioner> >(m, false) ));
}
template<typename MatrixType> void jacobisvd_verify_assert(const MatrixType& m)
{
svd_verify_assert<JacobiSVD<MatrixType> >(m);
svd_verify_assert<JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner> >(m, true);
svd_verify_assert<JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner> >(m);
svd_verify_assert<JacobiSVD<MatrixType, HouseholderQRPreconditioner> >(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<MatrixType, FullPivHouseholderQRPreconditioner> 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<typename MatrixType>
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<ComputeFullU | ComputeFullV>().solve(m), m);
VERIFY_IS_APPROX(m.template jacobiSvd<ComputeFullU | ComputeFullV>().transpose().solve(m), m);
VERIFY_IS_APPROX(m.template jacobiSvd<ComputeFullU | ComputeFullV>().adjoint().solve(m), m);
}
template <typename MatrixType>
void jacobisvd_all_options(const MatrixType& input = MatrixType()) {
MatrixType m = input;
svd_fill_random(m);
svd_option_checks<MatrixType, 0>(m);
svd_option_checks<MatrixType, ColPivHouseholderQRPreconditioner>(m);
svd_option_checks<MatrixType, HouseholderQRPreconditioner>(m);
svd_option_checks_full_only<MatrixType, FullPivHouseholderQRPreconditioner>(
m); // FullPiv only used when computing full unitaries
}
template <typename MatrixType>
void jacobisvd_verify_assert(const MatrixType& m = MatrixType()) {
svd_verify_assert<MatrixType, 0>(m);
svd_verify_assert<MatrixType, ColPivHouseholderQRPreconditioner>(m);
svd_verify_assert<MatrixType, HouseholderQRPreconditioner>(m);
svd_verify_assert_full_only<MatrixType, FullPivHouseholderQRPreconditioner>(m);
svd_verify_constructor_options_assert<JacobiSVD<MatrixType>>(m);
svd_verify_constructor_options_assert<JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>>(m);
svd_verify_constructor_options_assert<JacobiSVD<MatrixType, HouseholderQRPreconditioner>>(m);
svd_verify_constructor_options_assert<JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>>(m, true);
}
template <typename MatrixType>
void jacobisvd_verify_inputs(const MatrixType& m = MatrixType()) {
// check defaults
typedef JacobiSVD<MatrixType> 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<MatrixType, ComputeThinU>::QRPreconditioner == (int)ColPivHouseholderQRPreconditioner));
VERIFY(((int)JacobiSVD<MatrixType, ComputeThinV>::QRPreconditioner == (int)ColPivHouseholderQRPreconditioner));
VERIFY(((int)JacobiSVD<MatrixType, ComputeThinU | ComputeThinV>::QRPreconditioner ==
(int)ColPivHouseholderQRPreconditioner));
VERIFY(((int)JacobiSVD<MatrixType, ComputeFullU | ComputeFullV>::QRPreconditioner ==
(int)ColPivHouseholderQRPreconditioner));
VERIFY(((int)JacobiSVD<MatrixType, ComputeThinU | ComputeFullV>::QRPreconditioner ==
(int)ColPivHouseholderQRPreconditioner));
VERIFY(((int)JacobiSVD<MatrixType, ComputeFullU | ComputeThinV>::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<Matrix2cd>));
CALL_SUBTEST_12(svd_all_trivial_2x2(jacobisvd<Matrix2d>));
CALL_SUBTEST_4((jacobisvd_verify_inputs<Matrix4d>()));
CALL_SUBTEST_7((jacobisvd_verify_inputs(Matrix<float, 10, Dynamic>(10, 12))));
CALL_SUBTEST_8((jacobisvd_verify_inputs<Matrix<std::complex<double>, 7, 5>>()));
for(int i = 0; i < g_repeat; i++) {
CALL_SUBTEST_3(( jacobisvd<Matrix3f>() ));
CALL_SUBTEST_4(( jacobisvd<Matrix4d>() ));
CALL_SUBTEST_5(( jacobisvd<Matrix<float,3,5> >() ));
CALL_SUBTEST_6(( jacobisvd<Matrix<double,Dynamic,2> >(Matrix<double,Dynamic,2>(10,2)) ));
CALL_SUBTEST_3((jacobisvd_verify_assert<Matrix3f>()));
CALL_SUBTEST_4((jacobisvd_verify_assert<Matrix4d>()));
CALL_SUBTEST_7((jacobisvd_verify_assert<Matrix<float, 10, 12>>()));
CALL_SUBTEST_7((jacobisvd_verify_assert<Matrix<float, 12, 10>>()));
CALL_SUBTEST_7((jacobisvd_verify_assert<MatrixXf>(MatrixXf(10, 12))));
CALL_SUBTEST_8((jacobisvd_verify_assert<MatrixXcd>(MatrixXcd(7, 5))));
CALL_SUBTEST_11(svd_all_trivial_2x2(jacobisvd_all_options<Matrix2cd>));
CALL_SUBTEST_12(svd_all_trivial_2x2(jacobisvd_all_options<Matrix2d>));
for (int i = 0; i < g_repeat; i++) {
int r = internal::random<int>(1, 30),
c = internal::random<int>(1, 30);
TEST_SET_BUT_UNUSED_VARIABLE(r)
TEST_SET_BUT_UNUSED_VARIABLE(c)
CALL_SUBTEST_10(( jacobisvd<MatrixXd>(MatrixXd(r,c)) ));
CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(r,c)) ));
CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(r,c)) ));
(void) r;
(void) c;
CALL_SUBTEST_3((jacobisvd_all_options<Matrix3f>()));
CALL_SUBTEST_3((jacobisvd_all_options<Matrix<float, 2, 3>>()));
CALL_SUBTEST_4((jacobisvd_all_options<Matrix4d>()));
CALL_SUBTEST_4((jacobisvd_all_options<Matrix<double, 10, 16>>()));
CALL_SUBTEST_4((jacobisvd_all_options<Matrix<double, 16, 10>>()));
CALL_SUBTEST_5((jacobisvd_all_options<Matrix<double, Dynamic, 16>>(Matrix<double, Dynamic, 16>(r, 16))));
CALL_SUBTEST_5((jacobisvd_all_options<Matrix<double, 10, Dynamic>>(Matrix<double, 10, Dynamic>(10, c))));
CALL_SUBTEST_7((jacobisvd_all_options<MatrixXf>(MatrixXf(r, c))));
CALL_SUBTEST_8((jacobisvd_all_options<MatrixXcd>(MatrixXcd(r, c))));
CALL_SUBTEST_10((jacobisvd_all_options<MatrixXd>(MatrixXd(r, c))));
CALL_SUBTEST_14((jacobisvd_all_options<Matrix<double, 5, 7, RowMajor>>()));
CALL_SUBTEST_14((jacobisvd_all_options<Matrix<double, 7, 5, RowMajor>>()));
MatrixXcd noQRTest = MatrixXcd(r, r);
svd_fill_random(noQRTest);
CALL_SUBTEST_16((svd_option_checks<MatrixXcd, NoQRPreconditioner>(noQRTest)));
CALL_SUBTEST_15((
svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, ColMajor, 13, 15>, ColPivHouseholderQRPreconditioner>(
r, c)));
CALL_SUBTEST_15(
(svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, ColMajor, 15, 13>, HouseholderQRPreconditioner>(r,
c)));
CALL_SUBTEST_15((
svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, RowMajor, 13, 15>, ColPivHouseholderQRPreconditioner>(
r, c)));
CALL_SUBTEST_15(
(svd_check_max_size_matrix<Matrix<float, Dynamic, Dynamic, RowMajor, 15, 13>, HouseholderQRPreconditioner>(r,
c)));
// Test on inf/nan matrix
CALL_SUBTEST_7( (svd_inf_nan<JacobiSVD<MatrixXf>, MatrixXf>()) );
CALL_SUBTEST_10( (svd_inf_nan<JacobiSVD<MatrixXd>, MatrixXd>()) );
CALL_SUBTEST_7((svd_inf_nan<MatrixXf>()));
CALL_SUBTEST_10((svd_inf_nan<MatrixXd>()));
// bug1395 test compile-time vectors as input
CALL_SUBTEST_13(( jacobisvd_verify_assert(Matrix<double,6,1>()) ));
CALL_SUBTEST_13(( jacobisvd_verify_assert(Matrix<double,1,6>()) ));
CALL_SUBTEST_13(( jacobisvd_verify_assert(Matrix<double,Dynamic,1>(r)) ));
CALL_SUBTEST_13(( jacobisvd_verify_assert(Matrix<double,1,Dynamic>(c)) ));
CALL_SUBTEST_13((jacobisvd_verify_assert<Matrix<double, 6, 1>>()));
CALL_SUBTEST_13((jacobisvd_verify_assert<Matrix<double, 1, 6>>()));
CALL_SUBTEST_13((jacobisvd_verify_assert<Matrix<double, Dynamic, 1>>(Matrix<double, Dynamic, 1>(r))));
CALL_SUBTEST_13((jacobisvd_verify_assert<Matrix<double, 1, Dynamic>>(Matrix<double, 1, Dynamic>(c))));
}
CALL_SUBTEST_7(( jacobisvd<MatrixXf>(MatrixXf(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/2))) ));
CALL_SUBTEST_8(( jacobisvd<MatrixXcd>(MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3), internal::random<int>(EIGEN_TEST_MAX_SIZE/4, EIGEN_TEST_MAX_SIZE/3))) ));
CALL_SUBTEST_7((jacobisvd_all_options<MatrixXd>(
MatrixXd(internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 2),
internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 2)))));
CALL_SUBTEST_8((jacobisvd_all_options<MatrixXcd>(
MatrixXcd(internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 3),
internal::random<int>(EIGEN_TEST_MAX_SIZE / 4, EIGEN_TEST_MAX_SIZE / 3)))));
// test matrixbase method
CALL_SUBTEST_1(( jacobisvd_method<Matrix2cd>() ));

2
test/nomalloc.cpp
View File

@ -152,7 +152,7 @@ void ctms_decompositions()
x = fpQR.solve(b);
// SVD module
Eigen::JacobiSVD<Matrix> jSVD; jSVD.compute(A, ComputeFullU | ComputeFullV);
Eigen::JacobiSVD<Matrix, ComputeFullU | ComputeFullV> jSVD; jSVD.compute(A);
}
void test_zerosized() {

4
test/qr_colpivoting.cpp
View File

@ -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<MatrixType> svd(matrix, ComputeThinU | ComputeThinV);
JacobiSVD<MatrixType, ComputeThinU | ComputeThinV> 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<Scalar, Rows, Cols2> rhs = matrix * exact_solution;
Matrix<Scalar, Cols, Cols2> cod_solution = cod.solve(rhs);
JacobiSVD<MatrixType> svd(matrix, ComputeFullU | ComputeFullV);
JacobiSVD<MatrixType, ComputeFullU | ComputeFullV> svd(matrix);
Matrix<Scalar, Cols, Cols2> svd_solution = svd.solve(rhs);
VERIFY_IS_APPROX(cod_solution, svd_solution);

421
test/svd_common.h
View File

@ -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<typename SvdType, typename MatrixType>
void svd_compare_to_full(const MatrixType& m,
unsigned int computationOptions,
const SvdType& referenceSvd)
{
template <typename MatrixType, typename SvdType, int Options>
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<RealScalar>();
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<typename SvdType, typename MatrixType>
void svd_least_square(const MatrixType& m, unsigned int computationOptions)
{
template <typename SvdType, typename MatrixType>
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<Scalar, ColsAtCompileTime, Dynamic> SolutionType;
RhsType rhs = RhsType::Random(rows, internal::random<Index>(1, cols));
SvdType svd(m, computationOptions);
SvdType svd(m);
if(internal::is_same<RealScalar,double>::value) svd.setThreshold(1e-8);
if (internal::is_same<RealScalar, double>::value)
svd.setThreshold(1e-8);
else if(internal::is_same<RealScalar,float>::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<typename MatrixType>
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 <typename MatrixType, int Options>
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<Scalar,RowsAtCompileTime3,Dynamic> C = Matrix<Scalar,RowsAtCompileTime3,Dynamic>::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<CMatrixType, MatrixType>(m, solver, rows, cols, cols2);
}
// Check full, compare_to_full, least_square, and min_norm for all possible compute-options
template<typename SvdType, typename MatrixType>
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<typename SvdType::StorageIndex,int>::value ));
SvdType fullSvd(m, ComputeFullU|ComputeFullV);
CALL_SUBTEST(( svd_check_full(m, fullSvd) ));
CALL_SUBTEST(( svd_least_square<SvdType>(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<SvdType>(m, ComputeFullU | ComputeThinV) ));
CALL_SUBTEST(( svd_least_square<SvdType>(m, ComputeThinU | ComputeFullV) ));
CALL_SUBTEST(( svd_least_square<SvdType>(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<typename Scalar>
@ -300,29 +247,28 @@ template<typename T> 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<typename SvdType, typename MatrixType>
void svd_inf_nan()
{
SvdType svd;
template <typename MatrixType>
void svd_inf_nan() {
SVD_STATIC_OPTIONS(MatrixType, ComputeFullU | ComputeFullV) svd;
typedef typename MatrixType::Scalar Scalar;
Scalar some_inf = Scalar(1) / zero<Scalar>();
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<Scalar>::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<int>(0,9), internal::random<int>(0,9)) = some_inf;
svd.compute(m, ComputeFullU | ComputeFullV);
svd.compute(m);
VERIFY(svd.info() == InvalidInput);
m = MatrixType::Zero(10,10);
m(internal::random<int>(0,9), internal::random<int>(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<Scalar>::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<typename MatrixType>
void svd_all_trivial_2x2( void (*cb)(const MatrixType&,bool) )
{
template <typename MatrixType>
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<typename SvdType,typename MatrixType>
void svd_verify_assert(const MatrixType& m, bool fullOnly = false)
{
template <typename MatrixType, int QRPreconditioner = 0>
void svd_verify_assert_full_only(const MatrixType& m = MatrixType()) {
enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime };
typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, 1> 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 <typename MatrixType, int QRPreconditioner = 0>
void svd_verify_assert(const MatrixType& m = MatrixType()) {
enum { RowsAtCompileTime = MatrixType::RowsAtCompileTime };
typedef Matrix<typename MatrixType::Scalar, RowsAtCompileTime, 1> 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<MatrixType, QRPreconditioner>(m);
}
template <typename MatrixType, int Options>
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<SVDType, MatrixType>(m);
// svd_min_norm generates non-square matrices so it can't be used with NoQRPreconditioner
if ((Options & internal::QRPreconditionerBits) != NoQRPreconditioner) svd_min_norm<MatrixType, Options>(m);
}
}
template <typename SvdType, typename MatrixType>
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 <typename MatrixType, int QRPreconditioner = 0>
void svd_option_checks(const MatrixType& m) {
svd_compute_checks<MatrixType, QRPreconditioner>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeThinU>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeThinV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeThinU | ComputeThinV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullU>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullU | ComputeFullV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeThinU | ComputeFullV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullU | ComputeThinV>(m);
typedef SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeFullU | ComputeFullV) FullSvdType;
FullSvdType fullSvd(m);
svd_check_full(m, fullSvd);
svd_compare_to_full<MatrixType, FullSvdType, QRPreconditioner | ComputeFullU | ComputeFullV>(m, fullSvd);
typedef SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner) DynamicSvd;
svd_check_constructor_options<DynamicSvd>(m, 0);
svd_check_constructor_options<DynamicSvd>(m, ComputeThinU);
svd_check_constructor_options<DynamicSvd>(m, ComputeThinV);
svd_check_constructor_options<DynamicSvd>(m, ComputeThinU | ComputeThinV);
svd_check_constructor_options<DynamicSvd>(m, ComputeFullU);
svd_check_constructor_options<DynamicSvd>(m, ComputeFullV);
svd_check_constructor_options<DynamicSvd>(m, ComputeFullU | ComputeFullV);
svd_check_constructor_options<DynamicSvd>(m, ComputeThinU | ComputeFullV);
svd_check_constructor_options<DynamicSvd>(m, ComputeFullU | ComputeThinV);
}
template <typename MatrixType, int QRPreconditioner = 0>
void svd_option_checks_full_only(const MatrixType& m) {
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullU>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullV>(m);
svd_compute_checks<MatrixType, QRPreconditioner | ComputeFullU | ComputeFullV>(m);
SVD_STATIC_OPTIONS(MatrixType, QRPreconditioner | ComputeFullU | ComputeFullV) fullSvd(m);
svd_check_full(m, fullSvd);
}
template <typename MatrixType, int QRPreconditioner = 0>
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<typename MatrixType::Scalar> 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 <typename SvdType, typename MatrixType>
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