diff options
Diffstat (limited to 'Eigen/src/UmfPackSupport/UmfPackSupport.h')
| -rw-r--r-- | Eigen/src/UmfPackSupport/UmfPackSupport.h | 251 |
1 files changed, 117 insertions, 134 deletions
diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h index 29c60c378..dc74de935 100644 --- a/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -16,6 +16,13 @@ namespace Eigen { // 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_free_numeric(void **Numeric, double) { umfpack_di_free_numeric(Numeric); *Numeric = 0; } @@ -107,15 +114,6 @@ inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *N 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 +126,46 @@ namespace internal { * 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; - 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 +177,8 @@ class UmfPackLU : internal::noncopyable 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. * @@ -207,37 +224,11 @@ class UmfPackLU : internal::noncopyable { 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. @@ -250,11 +241,44 @@ class UmfPackLU : internal::noncopyable 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 +292,14 @@ class UmfPackLU : internal::noncopyable if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); - grapInput(matrix.derived()); + grab(matrix.derived()); factorize_impl(); } - #ifndef EIGEN_PARSED_BY_DOXYGEN /** \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,52 +313,17 @@ class UmfPackLU : internal::noncopyable m_isInitialized = false; m_numeric = 0; m_symbolic = 0; - m_outerIndexPtr = 0; - m_innerIndexPtr = 0; - m_valuePtr = 0; m_extractedDataAreDirty = true; } - 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() { + umfpack_defaults(m_control.data(), Scalar()); int errorCode = 0; - errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, - &m_symbolic, 0, 0); + 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(), 0); m_isInitialized = true; m_info = errorCode ? InvalidInput : Success; @@ -347,36 +334,52 @@ class UmfPackLU : internal::noncopyable void factorize_impl() { - int errorCode; - errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, - m_symbolic, &m_numeric, 0, 0); + m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), + m_symbolic, &m_numeric, m_control.data(), 0); - m_info = errorCode ? NumericalIssue : Success; + 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 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 +421,30 @@ typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() cons 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(), 0); + if(x.innerStride()!=1) + x.col(j) = x_tmp; if (errorCode!=0) return false; } @@ -438,37 +452,6 @@ bool UmfPackLU<MatrixType>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDe 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 |