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/CholmodSupport/CMakeLists.txt b/Eigen/src/CholmodSupport/CMakeLists.txt
deleted file mode 100644
index 814dfa6..0000000
--- a/Eigen/src/CholmodSupport/CMakeLists.txt
+++ /dev/null
@@ -1,6 +0,0 @@
-FILE(GLOB Eigen_CholmodSupport_SRCS "*.h")
-
-INSTALL(FILES
- ${Eigen_CholmodSupport_SRCS}
- DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/CholmodSupport COMPONENT Devel
- )
diff --git a/Eigen/src/CholmodSupport/CholmodSupport.h b/Eigen/src/CholmodSupport/CholmodSupport.h
index c449960..5719720 100644
--- a/Eigen/src/CholmodSupport/CholmodSupport.h
+++ b/Eigen/src/CholmodSupport/CholmodSupport.h
@@ -14,46 +14,52 @@
namespace internal {
-template<typename Scalar, typename CholmodType>
-void cholmod_configure_matrix(CholmodType& mat)
-{
- if (internal::is_same<Scalar,float>::value)
- {
- mat.xtype = CHOLMOD_REAL;
- mat.dtype = CHOLMOD_SINGLE;
- }
- else if (internal::is_same<Scalar,double>::value)
- {
+template<typename Scalar> struct cholmod_configure_matrix;
+
+template<> struct cholmod_configure_matrix<double> {
+ template<typename CholmodType>
+ static void run(CholmodType& mat) {
mat.xtype = CHOLMOD_REAL;
mat.dtype = CHOLMOD_DOUBLE;
}
- else if (internal::is_same<Scalar,std::complex<float> >::value)
- {
- mat.xtype = CHOLMOD_COMPLEX;
- mat.dtype = CHOLMOD_SINGLE;
- }
- else if (internal::is_same<Scalar,std::complex<double> >::value)
- {
+};
+
+template<> struct cholmod_configure_matrix<std::complex<double> > {
+ template<typename CholmodType>
+ static void run(CholmodType& mat) {
mat.xtype = CHOLMOD_COMPLEX;
mat.dtype = CHOLMOD_DOUBLE;
}
- else
- {
- eigen_assert(false && "Scalar type not supported by CHOLMOD");
- }
-}
+};
+
+// Other scalar types are not yet suppotred by Cholmod
+// template<> struct cholmod_configure_matrix<float> {
+// template<typename CholmodType>
+// static void run(CholmodType& mat) {
+// mat.xtype = CHOLMOD_REAL;
+// mat.dtype = CHOLMOD_SINGLE;
+// }
+// };
+//
+// template<> struct cholmod_configure_matrix<std::complex<float> > {
+// template<typename CholmodType>
+// static void run(CholmodType& mat) {
+// mat.xtype = CHOLMOD_COMPLEX;
+// mat.dtype = CHOLMOD_SINGLE;
+// }
+// };
} // namespace internal
/** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object.
* Note that the data are shared.
*/
-template<typename _Scalar, int _Options, typename _Index>
-cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
+template<typename _Scalar, int _Options, typename _StorageIndex>
+cholmod_sparse viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_StorageIndex> > mat)
{
cholmod_sparse res;
res.nzmax = mat.nonZeros();
- res.nrow = mat.rows();;
+ res.nrow = mat.rows();
res.ncol = mat.cols();
res.p = mat.outerIndexPtr();
res.i = mat.innerIndexPtr();
@@ -74,11 +80,11 @@
res.dtype = 0;
res.stype = -1;
- if (internal::is_same<_Index,int>::value)
+ if (internal::is_same<_StorageIndex,int>::value)
{
res.itype = CHOLMOD_INT;
}
- else if (internal::is_same<_Index,UF_long>::value)
+ else if (internal::is_same<_StorageIndex,long>::value)
{
res.itype = CHOLMOD_LONG;
}
@@ -88,7 +94,7 @@
}
// setup res.xtype
- internal::cholmod_configure_matrix<_Scalar>(res);
+ internal::cholmod_configure_matrix<_Scalar>::run(res);
res.stype = 0;
@@ -98,16 +104,23 @@
template<typename _Scalar, int _Options, typename _Index>
const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
{
- cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
+ cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived()));
+ return res;
+}
+
+template<typename _Scalar, int _Options, typename _Index>
+const cholmod_sparse viewAsCholmod(const SparseVector<_Scalar,_Options,_Index>& mat)
+{
+ cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived()));
return res;
}
/** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix.
* The data are not copied but shared. */
template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
-cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
+cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<const SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
{
- cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
+ cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.matrix().const_cast_derived()));
if(UpLo==Upper) res.stype = 1;
if(UpLo==Lower) res.stype = -1;
@@ -131,19 +144,19 @@
res.x = (void*)(mat.derived().data());
res.z = 0;
- internal::cholmod_configure_matrix<Scalar>(res);
+ internal::cholmod_configure_matrix<Scalar>::run(res);
return res;
}
/** Returns a view of the Cholmod sparse matrix \a cm as an Eigen sparse matrix.
* The data are not copied but shared. */
-template<typename Scalar, int Flags, typename Index>
-MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
+template<typename Scalar, int Flags, typename StorageIndex>
+MappedSparseMatrix<Scalar,Flags,StorageIndex> viewAsEigen(cholmod_sparse& cm)
{
- return MappedSparseMatrix<Scalar,Flags,Index>
- (cm.nrow, cm.ncol, static_cast<Index*>(cm.p)[cm.ncol],
- static_cast<Index*>(cm.p), static_cast<Index*>(cm.i),static_cast<Scalar*>(cm.x) );
+ return MappedSparseMatrix<Scalar,Flags,StorageIndex>
+ (cm.nrow, cm.ncol, static_cast<StorageIndex*>(cm.p)[cm.ncol],
+ static_cast<StorageIndex*>(cm.p), static_cast<StorageIndex*>(cm.i),static_cast<Scalar*>(cm.x) );
}
enum CholmodMode {
@@ -157,29 +170,39 @@
* \sa class CholmodSupernodalLLT, class CholmodSimplicialLDLT, class CholmodSimplicialLLT
*/
template<typename _MatrixType, int _UpLo, typename Derived>
-class CholmodBase : internal::noncopyable
+class CholmodBase : public SparseSolverBase<Derived>
{
+ protected:
+ typedef SparseSolverBase<Derived> Base;
+ using Base::derived;
+ using Base::m_isInitialized;
public:
typedef _MatrixType MatrixType;
enum { UpLo = _UpLo };
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
typedef MatrixType CholMatrixType;
- typedef typename MatrixType::Index Index;
+ typedef typename MatrixType::StorageIndex StorageIndex;
+ enum {
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+ MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+ };
public:
CholmodBase()
- : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
+ : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false)
{
- m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
+ EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
+ m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
cholmod_start(&m_cholmod);
}
- CholmodBase(const MatrixType& matrix)
- : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
+ explicit CholmodBase(const MatrixType& matrix)
+ : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false)
{
- m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
+ EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
+ m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
cholmod_start(&m_cholmod);
compute(matrix);
}
@@ -191,11 +214,8 @@
cholmod_finish(&m_cholmod);
}
- inline Index cols() const { return m_cholmodFactor->n; }
- inline Index rows() const { return m_cholmodFactor->n; }
-
- Derived& derived() { return *static_cast<Derived*>(this); }
- const Derived& derived() const { return *static_cast<const Derived*>(this); }
+ inline StorageIndex cols() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
+ inline StorageIndex rows() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
/** \brief Reports whether previous computation was successful.
*
@@ -216,34 +236,6 @@
return derived();
}
- /** \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<CholmodBase, Rhs>
- solve(const MatrixBase<Rhs>& b) const
- {
- eigen_assert(m_isInitialized && "LLT is not initialized.");
- eigen_assert(rows()==b.rows()
- && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
- return internal::solve_retval<CholmodBase, Rhs>(*this, b.derived());
- }
-
- /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
- *
- * \sa compute()
- */
- template<typename Rhs>
- inline const internal::sparse_solve_retval<CholmodBase, Rhs>
- solve(const SparseMatrixBase<Rhs>& b) const
- {
- eigen_assert(m_isInitialized && "LLT is not initialized.");
- eigen_assert(rows()==b.rows()
- && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
- return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived());
- }
-
/** Performs a symbolic decomposition on the sparsity pattern of \a matrix.
*
* This function is particularly useful when solving for several problems having the same structure.
@@ -277,7 +269,7 @@
eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
-
+
// If the factorization failed, minor is the column at which it did. On success minor == n.
this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
m_factorizationIsOk = true;
@@ -290,20 +282,22 @@
#ifndef EIGEN_PARSED_BY_DOXYGEN
/** \internal */
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_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
const Index size = m_cholmodFactor->n;
EIGEN_UNUSED_VARIABLE(size);
eigen_assert(size==b.rows());
+
+ // Cholmod needs column-major stoarge without inner-stride, which corresponds to the default behavior of Ref.
+ Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b.derived());
- // note: cd stands for Cholmod Dense
- Rhs& b_ref(b.const_cast_derived());
cholmod_dense b_cd = viewAsCholmod(b_ref);
cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
if(!x_cd)
{
this->m_info = NumericalIssue;
+ return;
}
// TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
@@ -311,8 +305,8 @@
}
/** \internal */
- template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
- void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
+ template<typename RhsDerived, typename DestDerived>
+ void _solve_impl(const SparseMatrixBase<RhsDerived> &b, SparseMatrixBase<DestDerived> &dest) const
{
eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
const Index size = m_cholmodFactor->n;
@@ -320,14 +314,16 @@
eigen_assert(size==b.rows());
// note: cs stands for Cholmod Sparse
- cholmod_sparse b_cs = viewAsCholmod(b);
+ Ref<SparseMatrix<typename RhsDerived::Scalar,ColMajor,typename RhsDerived::StorageIndex> > b_ref(b.const_cast_derived());
+ cholmod_sparse b_cs = viewAsCholmod(b_ref);
cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
if(!x_cs)
{
this->m_info = NumericalIssue;
+ return;
}
// TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
- dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
+ dest.derived() = viewAsEigen<typename DestDerived::Scalar,ColMajor,typename DestDerived::StorageIndex>(*x_cs);
cholmod_free_sparse(&x_cs, &m_cholmod);
}
#endif // EIGEN_PARSED_BY_DOXYGEN
@@ -344,10 +340,61 @@
*/
Derived& setShift(const RealScalar& offset)
{
- m_shiftOffset[0] = offset;
+ m_shiftOffset[0] = double(offset);
return derived();
}
+ /** \returns the determinant of the underlying matrix from the current factorization */
+ Scalar determinant() const
+ {
+ using std::exp;
+ return exp(logDeterminant());
+ }
+
+ /** \returns the log determinant of the underlying matrix from the current factorization */
+ Scalar logDeterminant() const
+ {
+ using std::log;
+ using numext::real;
+ eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
+
+ RealScalar logDet = 0;
+ Scalar *x = static_cast<Scalar*>(m_cholmodFactor->x);
+ if (m_cholmodFactor->is_super)
+ {
+ // Supernodal factorization stored as a packed list of dense column-major blocs,
+ // as described by the following structure:
+
+ // super[k] == index of the first column of the j-th super node
+ StorageIndex *super = static_cast<StorageIndex*>(m_cholmodFactor->super);
+ // pi[k] == offset to the description of row indices
+ StorageIndex *pi = static_cast<StorageIndex*>(m_cholmodFactor->pi);
+ // px[k] == offset to the respective dense block
+ StorageIndex *px = static_cast<StorageIndex*>(m_cholmodFactor->px);
+
+ Index nb_super_nodes = m_cholmodFactor->nsuper;
+ for (Index k=0; k < nb_super_nodes; ++k)
+ {
+ StorageIndex ncols = super[k + 1] - super[k];
+ StorageIndex nrows = pi[k + 1] - pi[k];
+
+ Map<const Array<Scalar,1,Dynamic>, 0, InnerStride<> > sk(x + px[k], ncols, InnerStride<>(nrows+1));
+ logDet += sk.real().log().sum();
+ }
+ }
+ else
+ {
+ // Simplicial factorization stored as standard CSC matrix.
+ StorageIndex *p = static_cast<StorageIndex*>(m_cholmodFactor->p);
+ Index size = m_cholmodFactor->n;
+ for (Index k=0; k<size; ++k)
+ logDet += log(real( x[p[k]] ));
+ }
+ if (m_cholmodFactor->is_ll)
+ logDet *= 2.0;
+ return logDet;
+ };
+
template<typename Stream>
void dumpMemory(Stream& /*s*/)
{}
@@ -355,9 +402,8 @@
protected:
mutable cholmod_common m_cholmod;
cholmod_factor* m_cholmodFactor;
- RealScalar m_shiftOffset[2];
+ double m_shiftOffset[2];
mutable ComputationInfo m_info;
- bool m_isInitialized;
int m_factorizationIsOk;
int m_analysisIsOk;
};
@@ -376,9 +422,13 @@
* \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower.
*
+ * \implsparsesolverconcept
+ *
* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
*
- * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLLT
+ * \warning Only double precision real and complex scalar types are supported by Cholmod.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class CholmodSupernodalLLT, class SimplicialLLT
*/
template<typename _MatrixType, int _UpLo = Lower>
class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
@@ -395,7 +445,7 @@
CholmodSimplicialLLT(const MatrixType& matrix) : Base()
{
init();
- compute(matrix);
+ this->compute(matrix);
}
~CholmodSimplicialLLT() {}
@@ -423,9 +473,13 @@
* \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower.
*
+ * \implsparsesolverconcept
+ *
* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
*
- * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLDLT
+ * \warning Only double precision real and complex scalar types are supported by Cholmod.
+ *
+ * \sa \ref TutorialSparseSolverConcept, class CholmodSupernodalLLT, class SimplicialLDLT
*/
template<typename _MatrixType, int _UpLo = Lower>
class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
@@ -442,7 +496,7 @@
CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
{
init();
- compute(matrix);
+ this->compute(matrix);
}
~CholmodSimplicialLDLT() {}
@@ -468,9 +522,13 @@
* \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower.
*
+ * \implsparsesolverconcept
+ *
* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
*
- * \sa \ref TutorialSparseDirectSolvers
+ * \warning Only double precision real and complex scalar types are supported by Cholmod.
+ *
+ * \sa \ref TutorialSparseSolverConcept
*/
template<typename _MatrixType, int _UpLo = Lower>
class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
@@ -487,7 +545,7 @@
CholmodSupernodalLLT(const MatrixType& matrix) : Base()
{
init();
- compute(matrix);
+ this->compute(matrix);
}
~CholmodSupernodalLLT() {}
@@ -515,9 +573,13 @@
* \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
* or Upper. Default is Lower.
*
+ * \implsparsesolverconcept
+ *
* This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
*
- * \sa \ref TutorialSparseDirectSolvers
+ * \warning Only double precision real and complex scalar types are supported by Cholmod.
+ *
+ * \sa \ref TutorialSparseSolverConcept
*/
template<typename _MatrixType, int _UpLo = Lower>
class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
@@ -534,7 +596,7 @@
CholmodDecomposition(const MatrixType& matrix) : Base()
{
init();
- compute(matrix);
+ this->compute(matrix);
}
~CholmodDecomposition() {}
@@ -572,36 +634,6 @@
}
};
-namespace internal {
-
-template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
-struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
- : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
-{
- typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
- EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- dec()._solve(rhs(),dst);
- }
-};
-
-template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
-struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
- : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
-{
- typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
- EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- dec()._solve(rhs(),dst);
- }
-};
-
-} // end namespace internal
-
} // end namespace Eigen
#endif // EIGEN_CHOLMODSUPPORT_H