diff --git a/Eigen/src/Sparse/UmfPackSupport.h b/Eigen/src/Sparse/UmfPackSupport.h index d497e27d6..751e390f6 100644 --- a/Eigen/src/Sparse/UmfPackSupport.h +++ b/Eigen/src/Sparse/UmfPackSupport.h @@ -25,6 +25,51 @@ #ifndef EIGEN_UMFPACKSUPPORT_H #define EIGEN_UMFPACKSUPPORT_H +/* TODO extract L, extrac U, compute det, etc... */ + +inline int umfpack_symbolic(int n_row,int n_col, + const int Ap[], const int Ai[], const double Ax[], void **Symbolic, + const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) +{ + return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info); +} + +inline int umfpack_symbolic(int n_row,int n_col, + const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic, + const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO]) +{ + return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&Ax[0].real(),0,Symbolic,Control,Info); +} + +inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[], + void *Symbolic, void **Numeric, + const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) +{ + return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info); +} + +inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[], + void *Symbolic, void **Numeric, + const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO]) +{ + return umfpack_zi_numeric(Ap,Ai,&Ax[0].real(),0,Symbolic,Numeric,Control,Info); +} + +inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[], + double X[], const double B[], void *Numeric, + const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) +{ + return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info); +} + +inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[], + std::complex<double> X[], const std::complex<double> B[], void *Numeric, + const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) +{ + return umfpack_zi_solve(sys,Ap,Ai,&Ax[0].real(),0,&X[0].real(),0,&B[0].real(),0,Numeric,Control,Info); +} + + template<typename MatrixType> class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType> { @@ -69,8 +114,9 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType> template<typename MatrixType> void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a) { - const int size = a.rows(); - ei_assert((MatrixType::Flags&RowMajorBit)==0); + const int rows = a.rows(); + const int cols = a.cols(); + ei_assert((MatrixType::Flags&RowMajorBit)==0 && "Row major matrices are not supported yet"); m_matrixRef = &a; @@ -79,10 +125,10 @@ void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a) void* symbolic; int errorCode = 0; - errorCode = umfpack_di_symbolic(size, size, a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(), + errorCode = umfpack_symbolic(rows, cols, a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(), &symbolic, 0, 0); if (errorCode==0) - errorCode = umfpack_di_numeric(a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(), + errorCode = umfpack_numeric(a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(), symbolic, &m_numeric, 0, 0); umfpack_di_free_symbolic(&symbolic); @@ -119,7 +165,7 @@ bool SparseLU<MatrixType,UmfPack>::solve(const MatrixBase<BDerived> &b, MatrixBa int errorCode; for (int j=0; j<rhsCols; ++j) { - errorCode = umfpack_di_solve(UMFPACK_A, + errorCode = umfpack_solve(UMFPACK_A, m_matrixRef->_outerIndexPtr(), m_matrixRef->_innerIndexPtr(), m_matrixRef->_valuePtr(), &x->col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0); if (errorCode!=0)