UmfPack support: add support for complex<double>

1 files changed, 51 insertions, 5 deletions

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)