Squashed 'third_party/eigen/' changes from 61d72f6..cf794d3
Change-Id: I9b814151b01f49af6337a8605d0c42a3a1ed4c72
git-subtree-dir: third_party/eigen
git-subtree-split: cf794d3b741a6278df169e58461f8529f43bce5d
diff --git a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
index de00877..953d57c 100644
--- a/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
+++ b/Eigen/src/SPQRSupport/SuiteSparseQRSupport.h
@@ -2,6 +2,7 @@
// for linear algebra.
//
// Copyright (C) 2012 Desire Nuentsa <desire.nuentsa_wakam@inria.fr>
+// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
@@ -32,45 +33,54 @@
} // End namespace internal
/**
- * \ingroup SPQRSupport_Module
- * \class SPQR
- * \brief Sparse QR factorization based on SuiteSparseQR library
- *
- * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition
- * of sparse matrices. The result is then used to solve linear leasts_square systems.
- * Clearly, a QR factorization is returned such that A*P = Q*R where :
- *
- * P is the column permutation. Use colsPermutation() to get it.
- *
- * Q is the orthogonal matrix represented as Householder reflectors.
- * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
- * You can then apply it to a vector.
- *
- * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix.
- * NOTE : The Index type of R is always UF_long. You can get it with SPQR::Index
- *
- * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
- * NOTE
- *
- */
+ * \ingroup SPQRSupport_Module
+ * \class SPQR
+ * \brief Sparse QR factorization based on SuiteSparseQR library
+ *
+ * This class is used to perform a multithreaded and multifrontal rank-revealing QR decomposition
+ * of sparse matrices. The result is then used to solve linear leasts_square systems.
+ * Clearly, a QR factorization is returned such that A*P = Q*R where :
+ *
+ * P is the column permutation. Use colsPermutation() to get it.
+ *
+ * Q is the orthogonal matrix represented as Householder reflectors.
+ * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose.
+ * You can then apply it to a vector.
+ *
+ * R is the sparse triangular factor. Use matrixQR() to get it as SparseMatrix.
+ * NOTE : The Index type of R is always SuiteSparse_long. You can get it with SPQR::Index
+ *
+ * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<>
+ *
+ * \implsparsesolverconcept
+ *
+ *
+ */
template<typename _MatrixType>
-class SPQR
+class SPQR : public SparseSolverBase<SPQR<_MatrixType> >
{
+ protected:
+ typedef SparseSolverBase<SPQR<_MatrixType> > Base;
+ using Base::m_isInitialized;
public:
typedef typename _MatrixType::Scalar Scalar;
typedef typename _MatrixType::RealScalar RealScalar;
- typedef UF_long Index ;
- typedef SparseMatrix<Scalar, ColMajor, Index> MatrixType;
- typedef PermutationMatrix<Dynamic, Dynamic> PermutationType;
+ typedef SuiteSparse_long StorageIndex ;
+ typedef SparseMatrix<Scalar, ColMajor, StorageIndex> MatrixType;
+ typedef Map<PermutationMatrix<Dynamic, Dynamic, StorageIndex> > PermutationType;
+ enum {
+ ColsAtCompileTime = Dynamic,
+ MaxColsAtCompileTime = Dynamic
+ };
public:
SPQR()
- : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
+ : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
{
cholmod_l_start(&m_cc);
}
- SPQR(const _MatrixType& matrix)
- : m_isInitialized(false), m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
+ explicit SPQR(const _MatrixType& matrix)
+ : m_ordering(SPQR_ORDERING_DEFAULT), m_allow_tol(SPQR_DEFAULT_TOL), m_tolerance (NumTraits<Scalar>::epsilon()), m_useDefaultThreshold(true)
{
cholmod_l_start(&m_cc);
compute(matrix);
@@ -103,23 +113,22 @@
RealScalar pivotThreshold = m_tolerance;
if(m_useDefaultThreshold)
{
- using std::max;
RealScalar max2Norm = 0.0;
- for (int j = 0; j < mat.cols(); j++) max2Norm = (max)(max2Norm, mat.col(j).norm());
+ for (int j = 0; j < mat.cols(); j++) max2Norm = numext::maxi(max2Norm, mat.col(j).norm());
if(max2Norm==RealScalar(0))
max2Norm = RealScalar(1);
pivotThreshold = 20 * (mat.rows() + mat.cols()) * max2Norm * NumTraits<RealScalar>::epsilon();
}
-
cholmod_sparse A;
A = viewAsCholmod(mat);
+ m_rows = matrix.rows();
Index col = matrix.cols();
m_rank = SuiteSparseQR<Scalar>(m_ordering, pivotThreshold, col, &A,
&m_cR, &m_E, &m_H, &m_HPinv, &m_HTau, &m_cc);
if (!m_cR)
{
- m_info = NumericalIssue;
+ m_info = NumericalIssue;
m_isInitialized = false;
return;
}
@@ -130,28 +139,15 @@
/**
* Get the number of rows of the input matrix and the Q matrix
*/
- inline Index rows() const {return m_cR->nrow; }
+ inline Index rows() const {return m_rows; }
/**
* Get the number of columns of the input matrix.
*/
inline Index cols() const { return m_cR->ncol; }
-
- /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
- *
- * \sa compute()
- */
- template<typename Rhs>
- inline const internal::solve_retval<SPQR, Rhs> solve(const MatrixBase<Rhs>& B) const
- {
- eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
- eigen_assert(this->rows()==B.rows()
- && "SPQR::solve(): invalid number of rows of the right hand side matrix B");
- return internal::solve_retval<SPQR, Rhs>(*this, B.derived());
- }
template<typename Rhs, typename Dest>
- void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
+ void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
{
eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
eigen_assert(b.cols()==1 && "This method is for vectors only");
@@ -184,7 +180,7 @@
{
eigen_assert(m_isInitialized && " The QR factorization should be computed first, call compute()");
if(!m_isRUpToDate) {
- m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::Index>(*m_cR);
+ m_R = viewAsEigen<Scalar,ColMajor, typename MatrixType::StorageIndex>(*m_cR);
m_isRUpToDate = true;
}
return m_R;
@@ -198,11 +194,7 @@
PermutationType colsPermutation() const
{
eigen_assert(m_isInitialized && "Decomposition is not initialized.");
- Index n = m_cR->ncol;
- PermutationType colsPerm(n);
- for(Index j = 0; j <n; j++) colsPerm.indices()(j) = m_E[j];
- return colsPerm;
-
+ return PermutationType(m_E, m_cR->ncol);
}
/**
* Gets the rank of the matrix.
@@ -237,7 +229,6 @@
return m_info;
}
protected:
- bool m_isInitialized;
bool m_analysisIsOk;
bool m_factorizationIsOk;
mutable bool m_isRUpToDate;
@@ -247,13 +238,14 @@
RealScalar m_tolerance; // treat columns with 2-norm below this tolerance as zero
mutable cholmod_sparse *m_cR; // The sparse R factor in cholmod format
mutable MatrixType m_R; // The sparse matrix R in Eigen format
- mutable Index *m_E; // The permutation applied to columns
+ mutable StorageIndex *m_E; // The permutation applied to columns
mutable cholmod_sparse *m_H; //The householder vectors
- mutable Index *m_HPinv; // The row permutation of H
+ mutable StorageIndex *m_HPinv; // The row permutation of H
mutable cholmod_dense *m_HTau; // The Householder coefficients
mutable Index m_rank; // The rank of the matrix
mutable cholmod_common m_cc; // Workspace and parameters
bool m_useDefaultThreshold; // Use default threshold
+ Index m_rows;
template<typename ,typename > friend struct SPQR_QProduct;
};
@@ -261,7 +253,7 @@
struct SPQR_QProduct : ReturnByValue<SPQR_QProduct<SPQRType,Derived> >
{
typedef typename SPQRType::Scalar Scalar;
- typedef typename SPQRType::Index Index;
+ typedef typename SPQRType::StorageIndex StorageIndex;
//Define the constructor to get reference to argument types
SPQR_QProduct(const SPQRType& spqr, const Derived& other, bool transpose) : m_spqr(spqr),m_other(other),m_transpose(transpose) {}
@@ -317,22 +309,5 @@
const SPQRType& m_spqr;
};
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<SPQR<_MatrixType>, Rhs>
- : solve_retval_base<SPQR<_MatrixType>, Rhs>
-{
- typedef SPQR<_MatrixType> Dec;
- EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- dec()._solve(rhs(),dst);
- }
-};
-
-} // end namespace internal
-
}// End namespace Eigen
#endif