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/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h
index 29c60c3..91c09ab 100644
--- a/Eigen/src/UmfPackSupport/UmfPackSupport.h
+++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h
@@ -10,12 +10,37 @@
#ifndef EIGEN_UMFPACKSUPPORT_H
#define EIGEN_UMFPACKSUPPORT_H
-namespace Eigen {
+namespace Eigen {
/* TODO extract L, extract U, compute det, etc... */
// generic double/complex<double> wrapper functions:
+
+inline void umfpack_defaults(double control[UMFPACK_CONTROL], double)
+{ umfpack_di_defaults(control); }
+
+inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>)
+{ umfpack_zi_defaults(control); }
+
+inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double)
+{ umfpack_di_report_info(control, info);}
+
+inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>)
+{ umfpack_zi_report_info(control, info);}
+
+inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double)
+{ umfpack_di_report_status(control, status);}
+
+inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>)
+{ umfpack_zi_report_status(control, status);}
+
+inline void umfpack_report_control(double control[UMFPACK_CONTROL], double)
+{ umfpack_di_report_control(control);}
+
+inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>)
+{ umfpack_zi_report_control(control);}
+
inline void umfpack_free_numeric(void **Numeric, double)
{ umfpack_di_free_numeric(Numeric); *Numeric = 0; }
@@ -107,15 +132,6 @@
return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
}
-namespace internal {
- template<typename T> struct umfpack_helper_is_sparse_plain : false_type {};
- template<typename Scalar, int Options, typename StorageIndex>
- struct umfpack_helper_is_sparse_plain<SparseMatrix<Scalar,Options,StorageIndex> >
- : true_type {};
- template<typename Scalar, int Options, typename StorageIndex>
- struct umfpack_helper_is_sparse_plain<MappedSparseMatrix<Scalar,Options,StorageIndex> >
- : true_type {};
-}
/** \ingroup UmfPackSupport_Module
* \brief A sparse LU factorization and solver based on UmfPack
@@ -128,27 +144,47 @@
* Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix.
* \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
*
- * \sa \ref TutorialSparseDirectSolvers
+ * \implsparsesolverconcept
+ *
+ * \sa \ref TutorialSparseSolverConcept, class SparseLU
*/
template<typename _MatrixType>
-class UmfPackLU : internal::noncopyable
+class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
{
+ protected:
+ typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base;
+ using Base::m_isInitialized;
public:
+ using Base::_solve_impl;
typedef _MatrixType MatrixType;
typedef typename MatrixType::Scalar Scalar;
typedef typename MatrixType::RealScalar RealScalar;
- typedef typename MatrixType::Index Index;
+ typedef typename MatrixType::StorageIndex StorageIndex;
typedef Matrix<Scalar,Dynamic,1> Vector;
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
typedef SparseMatrix<Scalar> LUMatrixType;
typedef SparseMatrix<Scalar,ColMajor,int> UmfpackMatrixType;
+ typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef;
+ enum {
+ ColsAtCompileTime = MatrixType::ColsAtCompileTime,
+ MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
+ };
public:
- UmfPackLU() { init(); }
+ typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl;
+ typedef Array<double, UMFPACK_INFO, 1> UmfpackInfo;
- UmfPackLU(const MatrixType& matrix)
+ UmfPackLU()
+ : m_dummy(0,0), mp_matrix(m_dummy)
+ {
+ init();
+ }
+
+ template<typename InputMatrixType>
+ explicit UmfPackLU(const InputMatrixType& matrix)
+ : mp_matrix(matrix)
{
init();
compute(matrix);
@@ -160,8 +196,8 @@
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
}
- inline Index rows() const { return m_copyMatrix.rows(); }
- inline Index cols() const { return m_copyMatrix.cols(); }
+ inline Index rows() const { return mp_matrix.rows(); }
+ inline Index cols() const { return mp_matrix.cols(); }
/** \brief Reports whether previous computation was successful.
*
@@ -198,7 +234,7 @@
return m_q;
}
- /** Computes the sparse Cholesky decomposition of \a matrix
+ /** Computes the sparse Cholesky decomposition of \a matrix
* Note that the matrix should be column-major, and in compressed format for best performance.
* \sa SparseMatrix::makeCompressed().
*/
@@ -207,37 +243,11 @@
{
if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
- grapInput(matrix.derived());
+ grab(matrix.derived());
analyzePattern_impl();
factorize_impl();
}
- /** \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<UmfPackLU, Rhs> solve(const MatrixBase<Rhs>& b) const
- {
- eigen_assert(m_isInitialized && "UmfPackLU is not initialized.");
- eigen_assert(rows()==b.rows()
- && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
- return internal::solve_retval<UmfPackLU, 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<UmfPackLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
- {
- eigen_assert(m_isInitialized && "UmfPackLU is not initialized.");
- eigen_assert(rows()==b.rows()
- && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
- return internal::sparse_solve_retval<UmfPackLU, Rhs>(*this, b.derived());
- }
-
/** Performs a symbolic decomposition on the sparcity of \a matrix.
*
* This function is particularly useful when solving for several problems having the same structure.
@@ -249,12 +259,45 @@
{
if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
-
- grapInput(matrix.derived());
+
+ grab(matrix.derived());
analyzePattern_impl();
}
+ /** Provides the return status code returned by UmfPack during the numeric
+ * factorization.
+ *
+ * \sa factorize(), compute()
+ */
+ inline int umfpackFactorizeReturncode() const
+ {
+ eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()");
+ return m_fact_errorCode;
+ }
+
+ /** Provides access to the control settings array used by UmfPack.
+ *
+ * If this array contains NaN's, the default values are used.
+ *
+ * See UMFPACK documentation for details.
+ */
+ inline const UmfpackControl& umfpackControl() const
+ {
+ return m_control;
+ }
+
+ /** Provides access to the control settings array used by UmfPack.
+ *
+ * If this array contains NaN's, the default values are used.
+ *
+ * See UMFPACK documentation for details.
+ */
+ inline UmfpackControl& umfpackControl()
+ {
+ return m_control;
+ }
+
/** Performs a numeric decomposition of \a matrix
*
* The given matrix must has the same sparcity than the matrix on which the pattern anylysis has been performed.
@@ -268,16 +311,42 @@
if(m_numeric)
umfpack_free_numeric(&m_numeric,Scalar());
- grapInput(matrix.derived());
-
+ grab(matrix.derived());
+
factorize_impl();
}
- #ifndef EIGEN_PARSED_BY_DOXYGEN
+ /** Prints the current UmfPack control settings.
+ *
+ * \sa umfpackControl()
+ */
+ void umfpackReportControl()
+ {
+ umfpack_report_control(m_control.data(), Scalar());
+ }
+
+ /** Prints statistics collected by UmfPack.
+ *
+ * \sa analyzePattern(), compute()
+ */
+ void umfpackReportInfo()
+ {
+ eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
+ umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar());
+ }
+
+ /** Prints the status of the previous factorization operation performed by UmfPack (symbolic or numerical factorization).
+ *
+ * \sa analyzePattern(), compute()
+ */
+ void umfpackReportStatus() {
+ eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
+ umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar());
+ }
+
/** \internal */
template<typename BDerived,typename XDerived>
- bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
- #endif
+ bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
Scalar determinant() const;
@@ -291,92 +360,75 @@
m_isInitialized = false;
m_numeric = 0;
m_symbolic = 0;
- m_outerIndexPtr = 0;
- m_innerIndexPtr = 0;
- m_valuePtr = 0;
m_extractedDataAreDirty = true;
+
+ umfpack_defaults(m_control.data(), Scalar());
}
-
- template<typename InputMatrixType>
- void grapInput_impl(const InputMatrixType& mat, internal::true_type)
- {
- m_copyMatrix.resize(mat.rows(), mat.cols());
- if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() )
- {
- // non supported input -> copy
- m_copyMatrix = mat;
- m_outerIndexPtr = m_copyMatrix.outerIndexPtr();
- m_innerIndexPtr = m_copyMatrix.innerIndexPtr();
- m_valuePtr = m_copyMatrix.valuePtr();
- }
- else
- {
- m_outerIndexPtr = mat.outerIndexPtr();
- m_innerIndexPtr = mat.innerIndexPtr();
- m_valuePtr = mat.valuePtr();
- }
- }
-
- template<typename InputMatrixType>
- void grapInput_impl(const InputMatrixType& mat, internal::false_type)
- {
- m_copyMatrix = mat;
- m_outerIndexPtr = m_copyMatrix.outerIndexPtr();
- m_innerIndexPtr = m_copyMatrix.innerIndexPtr();
- m_valuePtr = m_copyMatrix.valuePtr();
- }
-
- template<typename InputMatrixType>
- void grapInput(const InputMatrixType& mat)
- {
- grapInput_impl(mat, internal::umfpack_helper_is_sparse_plain<InputMatrixType>());
- }
-
+
void analyzePattern_impl()
{
- int errorCode = 0;
- errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
- &m_symbolic, 0, 0);
+ m_fact_errorCode = umfpack_symbolic(internal::convert_index<int>(mp_matrix.rows()),
+ internal::convert_index<int>(mp_matrix.cols()),
+ mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
+ &m_symbolic, m_control.data(), m_umfpackInfo.data());
m_isInitialized = true;
- m_info = errorCode ? InvalidInput : Success;
+ m_info = m_fact_errorCode ? InvalidInput : Success;
m_analysisIsOk = true;
m_factorizationIsOk = false;
m_extractedDataAreDirty = true;
}
-
+
void factorize_impl()
{
- int errorCode;
- errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
- m_symbolic, &m_numeric, 0, 0);
- m_info = errorCode ? NumericalIssue : Success;
+ m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
+ m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data());
+
+ m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
m_factorizationIsOk = true;
m_extractedDataAreDirty = true;
}
+ template<typename MatrixDerived>
+ void grab(const EigenBase<MatrixDerived> &A)
+ {
+ mp_matrix.~UmfpackMatrixRef();
+ ::new (&mp_matrix) UmfpackMatrixRef(A.derived());
+ }
+
+ void grab(const UmfpackMatrixRef &A)
+ {
+ if(&(A.derived()) != &mp_matrix)
+ {
+ mp_matrix.~UmfpackMatrixRef();
+ ::new (&mp_matrix) UmfpackMatrixRef(A);
+ }
+ }
+
// cached data to reduce reallocation, etc.
mutable LUMatrixType m_l;
+ int m_fact_errorCode;
+ UmfpackControl m_control;
+ mutable UmfpackInfo m_umfpackInfo;
+
mutable LUMatrixType m_u;
mutable IntColVectorType m_p;
mutable IntRowVectorType m_q;
- UmfpackMatrixType m_copyMatrix;
- const Scalar* m_valuePtr;
- const int* m_outerIndexPtr;
- const int* m_innerIndexPtr;
+ UmfpackMatrixType m_dummy;
+ UmfpackMatrixRef mp_matrix;
+
void* m_numeric;
void* m_symbolic;
mutable ComputationInfo m_info;
- bool m_isInitialized;
int m_factorizationIsOk;
int m_analysisIsOk;
mutable bool m_extractedDataAreDirty;
-
+
private:
- UmfPackLU(UmfPackLU& ) { }
+ UmfPackLU(const UmfPackLU& ) { }
};
@@ -418,19 +470,30 @@
template<typename MatrixType>
template<typename BDerived,typename XDerived>
-bool UmfPackLU<MatrixType>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
+bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
{
- const int rhsCols = b.cols();
+ Index rhsCols = b.cols();
eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
-
+
int errorCode;
+ Scalar* x_ptr = 0;
+ Matrix<Scalar,Dynamic,1> x_tmp;
+ if(x.innerStride()!=1)
+ {
+ x_tmp.resize(x.rows());
+ x_ptr = x_tmp.data();
+ }
for (int j=0; j<rhsCols; ++j)
{
+ if(x.innerStride()==1)
+ x_ptr = &x.col(j).coeffRef(0);
errorCode = umfpack_solve(UMFPACK_A,
- m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
- &x.col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0);
+ mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
+ x_ptr, &b.const_cast_derived().col(j).coeffRef(0), m_numeric, m_control.data(), m_umfpackInfo.data());
+ if(x.innerStride()!=1)
+ x.col(j) = x_tmp;
if (errorCode!=0)
return false;
}
@@ -438,37 +501,6 @@
return true;
}
-
-namespace internal {
-
-template<typename _MatrixType, typename Rhs>
-struct solve_retval<UmfPackLU<_MatrixType>, Rhs>
- : solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
-{
- typedef UmfPackLU<_MatrixType> Dec;
- EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- dec()._solve(rhs(),dst);
- }
-};
-
-template<typename _MatrixType, typename Rhs>
-struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs>
- : sparse_solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
-{
- typedef UmfPackLU<_MatrixType> Dec;
- EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
-
- template<typename Dest> void evalTo(Dest& dst) const
- {
- this->defaultEvalTo(dst);
- }
-};
-
-} // end namespace internal
-
} // end namespace Eigen
#endif // EIGEN_UMFPACKSUPPORT_H