From b509cf074224c470f62f356fe7cbdfc572ad8978 Mon Sep 17 00:00:00 2001 From: Gael Guennebaud Date: Fri, 1 Jun 2012 16:31:36 +0200 Subject: [PATCH] Fix bug #468: generalize UmfPack support to accept any input at the cost of an implicit copy. --- Eigen/src/UmfPackSupport/UmfPackSupport.h | 63 ++++++++++++++++------- 1 file changed, 45 insertions(+), 18 deletions(-) diff --git a/Eigen/src/UmfPackSupport/UmfPackSupport.h b/Eigen/src/UmfPackSupport/UmfPackSupport.h index 88d9d0201..9b09241b9 100644 --- a/Eigen/src/UmfPackSupport/UmfPackSupport.h +++ b/Eigen/src/UmfPackSupport/UmfPackSupport.h @@ -126,11 +126,11 @@ inline int umfpack_get_determinant(std::complex *Mx, double *Ex, void *N * \brief A sparse LU factorization and solver based on UmfPack * * This class allows to solve for A.X = B sparse linear problems via a LU factorization - * using the UmfPack library. The sparse matrix A must be in a compressed column-major form, squared and full rank. + * using the UmfPack library. The sparse matrix A must be squared and full rank. * The vectors or matrices X and B can be either dense or sparse. * - * WARNING The Eigen column-major SparseMatrix is not always in compressed form. - * The user should call makeCompressed() to get a matrix in CSC suitable for UMFPACK + * \WARNING The input matrix A should be in a \b compressed and \b column-major form. + * 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 @@ -147,6 +147,7 @@ class UmfPackLU typedef Matrix IntRowVectorType; typedef Matrix IntColVectorType; typedef SparseMatrix LUMatrixType; + typedef SparseMatrix UmfpackMatrixType; public: @@ -164,8 +165,8 @@ class UmfPackLU if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); } - inline Index rows() const { return m_matrixRef->rows(); } - inline Index cols() const { return m_matrixRef->cols(); } + inline Index rows() const { return m_copyMatrix.rows(); } + inline Index cols() const { return m_copyMatrix.cols(); } /** \brief Reports whether previous computation was successful. * @@ -203,7 +204,8 @@ class UmfPackLU } /** Computes the sparse Cholesky decomposition of \a matrix - * Note that the matrix should be in compressed format. Please, use makeCompressed() to get it !! + * Note that the matrix should be column-major, and in compressed format for best performance. + * \sa SparseMatrix::makeCompressed(). */ void compute(const MatrixType& matrix) { @@ -218,9 +220,9 @@ class UmfPackLU template inline const internal::solve_retval solve(const MatrixBase& b) const { - eigen_assert(m_isInitialized && "UmfPAckLU is not initialized."); + 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"); + && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b"); return internal::solve_retval(*this, b.derived()); } @@ -241,19 +243,19 @@ class UmfPackLU * * This function is particularly useful when solving for several problems having the same structure. * - * \sa factorize() + * \sa factorize(), compute() */ void analyzePattern(const MatrixType& matrix) { - eigen_assert((MatrixType::Flags&RowMajorBit)==0 && "UmfPackLU: Row major matrices are not supported yet"); - if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar()); if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); + + grapInput(matrix); int errorCode = 0; - errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), matrix.outerIndexPtr(), matrix.innerIndexPtr(), matrix.valuePtr(), + errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, &m_symbolic, 0, 0); m_isInitialized = true; @@ -264,9 +266,9 @@ class UmfPackLU /** Performs a numeric decomposition of \a matrix * - * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. + * The given matrix must has the same sparcity than the matrix on which the pattern anylysis has been performed. * - * \sa analyzePattern() + * \sa analyzePattern(), compute() */ void factorize(const MatrixType& matrix) { @@ -274,10 +276,10 @@ class UmfPackLU if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar()); - m_matrixRef = &matrix; + grapInput(matrix); int errorCode; - errorCode = umfpack_numeric(matrix.outerIndexPtr(), matrix.innerIndexPtr(), matrix.valuePtr(), + errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, m_symbolic, &m_numeric, 0, 0); m_info = errorCode ? NumericalIssue : Success; @@ -303,6 +305,28 @@ class UmfPackLU m_isInitialized = false; m_numeric = 0; m_symbolic = 0; + m_outerIndexPtr = 0; + m_innerIndexPtr = 0; + m_valuePtr = 0; + } + + void grapInput(const MatrixType& mat) + { + 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(); + } } // cached data to reduce reallocation, etc. @@ -311,7 +335,10 @@ class UmfPackLU mutable IntColVectorType m_p; mutable IntRowVectorType m_q; - const MatrixType* m_matrixRef; + UmfpackMatrixType m_copyMatrix; + const Scalar* m_valuePtr; + const int* m_outerIndexPtr; + const int* m_innerIndexPtr; void* m_numeric; void* m_symbolic; @@ -374,7 +401,7 @@ bool UmfPackLU::_solve(const MatrixBase &b, MatrixBaseouterIndexPtr(), m_matrixRef->innerIndexPtr(), m_matrixRef->valuePtr(), + m_outerIndexPtr, m_innerIndexPtr, m_valuePtr, &x.col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0); if (errorCode!=0) return false;