add new interface to SuperLU

unsupported/Eigen/SuperLUSupport

@@ -24,10 +24,11 @@ namespace Eigen {
   * \endcode
   */
 
-struct SuperLU {};
-
 #include "src/SparseExtra/SuperLUSupport.h"
 
+struct SuperLULegacy {};
+#include "src/SparseExtra/SuperLUSupportLegacy.h"
+
 } // namespace Eigen
 
 #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"

unsupported/Eigen/src/SparseExtra/SuperLUSupport.h

@@ -194,7 +194,7 @@ struct SluMatrix : SuperMatrix
       res.ncol      = mat.cols();
     }
 
-    res.Mtype     = SLU_GE;
+    res.Mtype       = SLU_GE;
 
     res.storage.nnz       = mat.nonZeros();
     res.storage.values    = mat.derived()._valuePtr();
@@ -252,7 +252,7 @@ struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
       res.ncol      = mat.cols();
     }
 
-    res.Mtype     = SLU_GE;
+    res.Mtype       = SLU_GE;
 
     res.storage.nnz       = mat.nonZeros();
     res.storage.values    = mat._valuePtr();
@@ -295,83 +295,146 @@ MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
 
 } // end namespace internal
 
-template<typename MatrixType>
+
-class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
+template<typename _MatrixType, typename Derived>
+class SuperLUBase
 {
-  protected:
+  public:
-    typedef SparseLU<MatrixType> Base;
+    typedef _MatrixType MatrixType;
-    typedef typename Base::Scalar Scalar;
+    typedef typename MatrixType::Scalar Scalar;
-    typedef typename Base::RealScalar RealScalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef typename MatrixType::Index Index;
     typedef Matrix<Scalar,Dynamic,1> Vector;
     typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
-    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
+    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;    
-    typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType;
+    typedef SparseMatrix<Scalar> LUMatrixType;
-    typedef SparseMatrix<Scalar,Upper> UMatrixType;
-    using Base::m_flags;
-    using Base::m_status;
 
   public:
 
-    SparseLU(int flags = NaturalOrdering)
+    SuperLUBase() {}
-      : Base(flags)
-    {
-    }
 
-    SparseLU(const MatrixType& matrix, int flags = NaturalOrdering)
+    ~SuperLUBase()
-      : Base(flags)
-    {
-      compute(matrix);
-    }
-
-    ~SparseLU()
     {
       Destroy_SuperNode_Matrix(&m_sluL);
       Destroy_CompCol_Matrix(&m_sluU);
     }
-
+    
-    inline const LMatrixType& matrixL() const
+    Derived& derived() { return *static_cast<Derived*>(this); }
+    const Derived& derived() const { return *static_cast<const Derived*>(this); }
+    
+    inline Index rows() const { return m_matrix.rows(); }
+    inline Index cols() const { return m_matrix.cols(); }
+    
+    /** \returns a reference to the Super LU option object to configure the  Super LU algorithms. */
+    inline superlu_options_t& options() { return m_sluOptions; }
+    
+    /** \brief Reports whether previous computation was successful.
+      *
+      * \returns \c Success if computation was succesful,
+      *          \c NumericalIssue if the matrix.appears to be negative.
+      */
+    ComputationInfo info() const
     {
-      if (m_extractedDataAreDirty) extractData();
+      eigen_assert(m_isInitialized && "Decomposition is not initialized.");
-      return m_l;
+      return m_info;
     }
 
-    inline const UMatrixType& matrixU() const
+    /** Computes the sparse Cholesky decomposition of \a matrix */
+    void compute(const MatrixType& matrix)
     {
-      if (m_extractedDataAreDirty) extractData();
+      derived().analyzePattern(matrix);
-      return m_u;
+      derived().factorize(matrix);
     }
-
+    
-    inline const IntColVectorType& permutationP() const
+    /** \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<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& b) const
     {
-      if (m_extractedDataAreDirty) extractData();
+      eigen_assert(m_isInitialized && "SuperLU is not initialized.");
-      return m_p;
+      eigen_assert(rows()==b.rows()
+                && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
+      return internal::solve_retval<SuperLUBase, Rhs>(*this, b.derived());
     }
-
+    
-    inline const IntRowVectorType& permutationQ() const
+    /** \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<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
+//     {
+//       eigen_assert(m_isInitialized && "SuperLU is not initialized.");
+//       eigen_assert(rows()==b.rows()
+//                 && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
+//       return internal::sparse_solve_retval<SuperLU, 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.
+      * 
+      * \sa factorize()
+      */
+    void analyzePattern(const MatrixType& /*matrix*/)
     {
-      if (m_extractedDataAreDirty) extractData();
+      m_isInitialized = true;
-      return m_q;
+      m_info = Success;
+      m_analysisIsOk = true;
+      m_factorizationIsOk = false;
     }
-
+    
-    Scalar determinant() const;
+    template<typename Stream>
-
+    void dumpMemory(Stream& s)
-    template<typename BDerived, typename XDerived>
+    {}
-    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed = SvNoTrans) const;
+    
-
-    void compute(const MatrixType& matrix);
-
   protected:
+    
+    void initFactorization(const MatrixType& a)
+    {
+      const int size = a.rows();
+      m_matrix = a;
 
+      m_sluA = internal::asSluMatrix(m_matrix);
+      memset(&m_sluL,0,sizeof m_sluL);
+      memset(&m_sluU,0,sizeof m_sluU);
+
+      m_p.resize(size);
+      m_q.resize(size);
+      m_sluRscale.resize(size);
+      m_sluCscale.resize(size);
+      m_sluEtree.resize(size);
+
+      // set empty B and X
+      m_sluB.setStorageType(SLU_DN);
+      m_sluB.setScalarType<Scalar>();
+      m_sluB.Mtype          = SLU_GE;
+      m_sluB.storage.values = 0;
+      m_sluB.nrow           = 0;
+      m_sluB.ncol           = 0;
+      m_sluB.storage.lda    = size;
+      m_sluX                = m_sluB;
+      
+      m_extractedDataAreDirty = true;
+    }
+    
+    void init()
+    {
+      m_info = InvalidInput;
+      m_isInitialized = false;
+    }
+    
     void extractData() const;
 
-  protected:
     // cached data to reduce reallocation, etc.
-    mutable LMatrixType m_l;
+    mutable LUMatrixType m_l;
-    mutable UMatrixType m_u;
+    mutable LUMatrixType m_u;
     mutable IntColVectorType m_p;
     mutable IntRowVectorType m_q;
 
-    mutable SparseMatrix<Scalar> m_matrix;
+    mutable LUMatrixType m_matrix;  // copy of the factorized matrix
     mutable SluMatrix m_sluA;
     mutable SuperMatrix m_sluL, m_sluU;
     mutable SluMatrix m_sluB, m_sluX;
@@ -381,171 +444,212 @@ class SparseLU<MatrixType,SuperLU> : public SparseLU<MatrixType>
     mutable std::vector<RealScalar> m_sluRscale, m_sluCscale;
     mutable std::vector<RealScalar> m_sluFerr, m_sluBerr;
     mutable char m_sluEqued;
+
+    mutable ComputationInfo m_info;
+    bool m_isInitialized;
+    int m_factorizationIsOk;
+    int m_analysisIsOk;
     mutable bool m_extractedDataAreDirty;
 };
 
-template<typename MatrixType>
+
-void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
+template<typename _MatrixType>
+class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
 {
-  const int size = a.rows();
+  public:
-  m_matrix = a;
+    typedef SuperLUBase<_MatrixType,SuperLU> Base;
+    typedef _MatrixType MatrixType;
+    typedef typename Base::Scalar Scalar;
+    typedef typename Base::RealScalar RealScalar;
+    typedef typename Base::Index Index;
+    typedef typename Base::IntRowVectorType IntRowVectorType;
+    typedef typename Base::IntColVectorType IntColVectorType;    
+    typedef typename Base::LUMatrixType LUMatrixType;
+    typedef TriangularView<LUMatrixType, Lower|UnitDiag>  LMatrixType;
+    typedef TriangularView<LUMatrixType,  Upper>           UMatrixType;
 
-  set_default_options(&m_sluOptions);
+  public:
-  m_sluOptions.ColPerm = NATURAL;
-  m_sluOptions.PrintStat = NO;
-  m_sluOptions.ConditionNumber = NO;
-  m_sluOptions.Trans = NOTRANS;
-  // m_sluOptions.Equil = NO;
 
-  switch (Base::orderingMethod())
+    SuperLU() : Base() { init(); }
+
+    SuperLU(const MatrixType& matrix) : Base()
+    {
+      init();
+      compute(matrix);
+    }
+
+    ~SuperLU()
+    {
+    }
+    
+    /** Performs a symbolic decomposition on the sparcity of \a matrix.
+      *
+      * This function is particularly useful when solving for several problems having the same structure.
+      * 
+      * \sa factorize()
+      */
+    void analyzePattern(const MatrixType& matrix)
+    {
+      Base::analyzePattern(matrix);
+    }
+    
+    /** 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.
+      *
+      * \sa analyzePattern()
+      */
+    void factorize(const MatrixType& matrix);
+    
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
+    /** \internal */
+    template<typename Rhs,typename Dest>
+    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
+    #endif // EIGEN_PARSED_BY_DOXYGEN
+    
+    inline const LMatrixType& matrixL() const
+    {
+      if (m_extractedDataAreDirty) this->extractData();
+      return m_l;
+    }
+
+    inline const UMatrixType& matrixU() const
+    {
+      if (m_extractedDataAreDirty) this->extractData();
+      return m_u;
+    }
+
+    inline const IntColVectorType& permutationP() const
+    {
+      if (m_extractedDataAreDirty) this->extractData();
+      return m_p;
+    }
+
+    inline const IntRowVectorType& permutationQ() const
+    {
+      if (m_extractedDataAreDirty) this->extractData();
+      return m_q;
+    }
+    
+    Scalar determinant() const;
+    
+  protected:
+    
+    using Base::m_matrix;
+    using Base::m_sluOptions;
+    using Base::m_sluA;
+    using Base::m_sluB;
+    using Base::m_sluX;
+    using Base::m_p;
+    using Base::m_q;
+    using Base::m_sluEtree;
+    using Base::m_sluEqued;
+    using Base::m_sluRscale;
+    using Base::m_sluCscale;
+    using Base::m_sluL;
+    using Base::m_sluU;
+    using Base::m_sluStat;
+    using Base::m_sluFerr;
+    using Base::m_sluBerr;
+    using Base::m_l;
+    using Base::m_u;
+    
+    using Base::m_analysisIsOk;
+    using Base::m_factorizationIsOk;
+    using Base::m_extractedDataAreDirty;
+    using Base::m_isInitialized;
+    using Base::m_info;
+    
+    void init()
+    {
+      Base::init();
+      
+      set_default_options(&this->m_sluOptions);
+      m_sluOptions.PrintStat        = NO;
+      m_sluOptions.ConditionNumber  = NO;
+      m_sluOptions.Trans            = NOTRANS;
+      m_sluOptions.ColPerm          = COLAMD;
+    }
+};
+
+template<typename MatrixType>
+void SuperLU<MatrixType>::factorize(const MatrixType& a)
+{
+  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
+  if(!m_analysisIsOk)
   {
-      case NaturalOrdering          : m_sluOptions.ColPerm = NATURAL; break;
+    m_info = InvalidInput;
-      case MinimumDegree_AT_PLUS_A  : m_sluOptions.ColPerm = MMD_AT_PLUS_A; break;
+    return;
-      case MinimumDegree_ATA        : m_sluOptions.ColPerm = MMD_ATA; break;
+  }
-      case ColApproxMinimumDegree   : m_sluOptions.ColPerm = COLAMD; break;
+  
-      default:
+  initFactorization(a);
-        //std::cerr << "Eigen: ordering method \"" << Base::orderingMethod() << "\" not supported by the SuperLU backend\n";
+  
-        m_sluOptions.ColPerm = NATURAL;
-  };
-
-  m_sluA = internal::asSluMatrix(m_matrix);
-  memset(&m_sluL,0,sizeof m_sluL);
-  memset(&m_sluU,0,sizeof m_sluU);
-  //m_sluEqued = 'B';
   int info = 0;
-
+  RealScalar recip_pivot_growth, rcond;
-  m_p.resize(size);
-  m_q.resize(size);
-  m_sluRscale.resize(size);
-  m_sluCscale.resize(size);
-  m_sluEtree.resize(size);
-
-  RealScalar recip_pivot_gross, rcond;
   RealScalar ferr, berr;
 
-  // set empty B and X
-  m_sluB.setStorageType(SLU_DN);
-  m_sluB.setScalarType<Scalar>();
-  m_sluB.Mtype = SLU_GE;
-  m_sluB.storage.values = 0;
-  m_sluB.nrow = m_sluB.ncol = 0;
-  m_sluB.storage.lda = size;
-  m_sluX = m_sluB;
-
   StatInit(&m_sluStat);
-  if (m_flags&IncompleteFactorization)
+  SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
-  {
+                &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
-    #ifdef EIGEN_SUPERLU_HAS_ILU
+                &m_sluL, &m_sluU,
-    ilu_set_default_options(&m_sluOptions);
+                NULL, 0,
-
+                &m_sluB, &m_sluX,
-    // no attempt to preserve column sum
+                &recip_pivot_growth, &rcond,
-    m_sluOptions.ILU_MILU = SILU;
+                &ferr, &berr,
-
+                &m_sluStat, &info, Scalar());
-    // only basic ILU(k) support -- no direct control over memory consumption
-    // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
-    // and set ILU_FillFactor to max memory growth
-    m_sluOptions.ILU_DropRule = DROP_BASIC;
-    m_sluOptions.ILU_DropTol = Base::m_precision;
-
-    SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
-      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
-      &m_sluL, &m_sluU,
-      NULL, 0,
-      &m_sluB, &m_sluX,
-      &recip_pivot_gross, &rcond,
-      &m_sluStat, &info, Scalar());
-    #else
-    //std::cerr << "Incomplete factorization is only available in SuperLU v4\n";
-    Base::m_succeeded = false;
-    return;
-    #endif
-  }
-  else
-  {
-    SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
-      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
-      &m_sluL, &m_sluU,
-      NULL, 0,
-      &m_sluB, &m_sluX,
-      &recip_pivot_gross, &rcond,
-      &ferr, &berr,
-      &m_sluStat, &info, Scalar());
-  }
   StatFree(&m_sluStat);
 
   m_extractedDataAreDirty = true;
 
   // FIXME how to better check for errors ???
-  Base::m_succeeded = (info == 0);
+  m_info = info == 0 ? Success : NumericalIssue;
+  m_factorizationIsOk = true;
 }
 
 template<typename MatrixType>
-template<typename BDerived,typename XDerived>
+template<typename Rhs,typename Dest>
-bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b,
+void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
-                        MatrixBase<XDerived> *x, const int transposed) const
 {
+  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
+
   const int size = m_matrix.rows();
   const int rhsCols = b.cols();
   eigen_assert(size==b.rows());
 
-  switch (transposed) {
+  m_sluOptions.Trans = NOTRANS;
-      case SvNoTrans    :  m_sluOptions.Trans = NOTRANS; break;
-      case SvTranspose  :  m_sluOptions.Trans = TRANS;   break;
-      case SvAdjoint    :  m_sluOptions.Trans = CONJ;    break;
-      default:
-        //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the SuperLU backend\n";
-        m_sluOptions.Trans = NOTRANS;
-  }
-
   m_sluOptions.Fact = FACTORED;
   m_sluOptions.IterRefine = NOREFINE;
+  
 
   m_sluFerr.resize(rhsCols);
   m_sluBerr.resize(rhsCols);
   m_sluB = SluMatrix::Map(b.const_cast_derived());
-  m_sluX = SluMatrix::Map(x->derived());
+  m_sluX = SluMatrix::Map(x.derived());
+  
+  typename Rhs::PlainObject b_cpy;
+  if(m_sluEqued!='N')
+  {
+    b_cpy = b;
+    m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
+  }
 
   StatInit(&m_sluStat);
   int info = 0;
-  RealScalar recip_pivot_gross, rcond;
+  RealScalar recip_pivot_growth, rcond;
-
+  SuperLU_gssvx(&m_sluOptions, &m_sluA,
-  if (m_flags&IncompleteFactorization)
+                m_q.data(), m_p.data(),
-  {
+                &m_sluEtree[0], &m_sluEqued,
-    #ifdef EIGEN_SUPERLU_HAS_ILU
+                &m_sluRscale[0], &m_sluCscale[0],
-    SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
+                &m_sluL, &m_sluU,
-      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
+                NULL, 0,
-      &m_sluL, &m_sluU,
+                &m_sluB, &m_sluX,
-      NULL, 0,
+                &recip_pivot_growth, &rcond,
-      &m_sluB, &m_sluX,
+                &m_sluFerr[0], &m_sluBerr[0],
-      &recip_pivot_gross, &rcond,
+                &m_sluStat, &info, Scalar());
-      &m_sluStat, &info, Scalar());
-    #else
-    //std::cerr << "Incomplete factorization is only available in SuperLU v4\n";
-    return false;
-    #endif
-  }
-  else
-  {
-    SuperLU_gssvx(
-      &m_sluOptions, &m_sluA,
-      m_q.data(), m_p.data(),
-      &m_sluEtree[0], &m_sluEqued,
-      &m_sluRscale[0], &m_sluCscale[0],
-      &m_sluL, &m_sluU,
-      NULL, 0,
-      &m_sluB, &m_sluX,
-      &recip_pivot_gross, &rcond,
-      &m_sluFerr[0], &m_sluBerr[0],
-      &m_sluStat, &info, Scalar());
-  }
   StatFree(&m_sluStat);
-
+  m_info = info==0 ? Success : NumericalIssue;
-  // reset to previous state
-  m_sluOptions.Trans = NOTRANS;
-  return info==0;
 }
 
-//
 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
 //
 //  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
@@ -553,9 +657,10 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b,
 //  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
 //  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
 //
-template<typename MatrixType>
+template<typename MatrixType, typename Derived>
-void SparseLU<MatrixType,SuperLU>::extractData() const
+void SuperLUBase<MatrixType,Derived>::extractData() const
 {
+  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
   if (m_extractedDataAreDirty)
   {
     int         upper;
@@ -639,13 +744,14 @@ void SparseLU<MatrixType,SuperLU>::extractData() const
 }
 
 template<typename MatrixType>
-typename SparseLU<MatrixType,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::determinant() const
+typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
 {
-  assert((!NumTraits<Scalar>::IsComplex) && "This function is not implemented for complex yet");
+  eigen_assert((!NumTraits<Scalar>::IsComplex) && "This function is not implemented for complex yet");
+  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for computing the determinant, you must first call either compute() or analyzePattern()/factorize()");
+  
   if (m_extractedDataAreDirty)
-    extractData();
+    this->extractData();
 
-  // TODO this code could be moved to the default/base backend
   // FIXME perhaps we have to take into account the scale factors m_sluRscale and m_sluCscale ???
   Scalar det = Scalar(1);
   for (int j=0; j<m_u.cols(); ++j)
@@ -659,9 +765,210 @@ typename SparseLU<MatrixType,SuperLU>::Scalar SparseLU<MatrixType,SuperLU>::dete
         det *= m_u._valuePtr()[lastId];
       }
     }
-//       std::cout << m_sluRscale[j] << " " << m_sluCscale[j] << "   \n";
   }
   return det;
 }
 
+#ifdef EIGEN_SUPERLU_HAS_ILU
+template<typename _MatrixType>
+class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
+{
+  public:
+    typedef SuperLUBase<_MatrixType,SuperILU> Base;
+    typedef _MatrixType MatrixType;
+    typedef typename Base::Scalar Scalar;
+    typedef typename Base::RealScalar RealScalar;
+    typedef typename Base::Index Index;
+
+  public:
+
+    SuperILU() : Base() { init(); }
+
+    SuperILU(const MatrixType& matrix) : Base()
+    {
+      init();
+      compute(matrix);
+    }
+
+    ~SuperILU()
+    {
+    }
+    
+    /** Performs a symbolic decomposition on the sparcity of \a matrix.
+      *
+      * This function is particularly useful when solving for several problems having the same structure.
+      * 
+      * \sa factorize()
+      */
+    void analyzePattern(const MatrixType& matrix)
+    {
+      Base::analyzePattern(matrix);
+    }
+    
+    /** 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.
+      *
+      * \sa analyzePattern()
+      */
+    void factorize(const MatrixType& matrix);
+    
+    #ifndef EIGEN_PARSED_BY_DOXYGEN
+    /** \internal */
+    template<typename Rhs,typename Dest>
+    void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
+    #endif // EIGEN_PARSED_BY_DOXYGEN
+    
+  protected:
+    
+    using Base::m_matrix;
+    using Base::m_sluOptions;
+    using Base::m_sluA;
+    using Base::m_sluB;
+    using Base::m_sluX;
+    using Base::m_p;
+    using Base::m_q;
+    using Base::m_sluEtree;
+    using Base::m_sluEqued;
+    using Base::m_sluRscale;
+    using Base::m_sluCscale;
+    using Base::m_sluL;
+    using Base::m_sluU;
+    using Base::m_sluStat;
+    using Base::m_sluFerr;
+    using Base::m_sluBerr;
+    using Base::m_l;
+    using Base::m_u;
+    
+    using Base::m_analysisIsOk;
+    using Base::m_factorizationIsOk;
+    using Base::m_extractedDataAreDirty;
+    using Base::m_isInitialized;
+    using Base::m_info;
+
+    void init()
+    {
+      Base::init();
+      
+      ilu_set_default_options(&m_sluOptions);
+      m_sluOptions.PrintStat        = NO;
+      m_sluOptions.ConditionNumber  = NO;
+      m_sluOptions.Trans            = NOTRANS;
+      m_sluOptions.ColPerm          = MMD_AT_PLUS_A;
+      
+      // no attempt to preserve column sum
+      m_sluOptions.ILU_MILU = SILU;
+      // only basic ILU(k) support -- no direct control over memory consumption
+      // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
+      // and set ILU_FillFactor to max memory growth
+      m_sluOptions.ILU_DropRule = DROP_BASIC;
+      m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
+    }
+};
+
+template<typename MatrixType>
+void SuperILU<MatrixType>::factorize(const MatrixType& a)
+{
+  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
+  if(!m_analysisIsOk)
+  {
+    m_info = InvalidInput;
+    return;
+  }
+  
+  this->initFactorization(a);
+
+  int info = 0;
+  RealScalar recip_pivot_growth, rcond;
+
+  StatInit(&m_sluStat);
+  SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
+                &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
+                &m_sluL, &m_sluU,
+                NULL, 0,
+                &m_sluB, &m_sluX,
+                &recip_pivot_growth, &rcond,
+                &m_sluStat, &info, Scalar());
+  StatFree(&m_sluStat);
+
+  // FIXME how to better check for errors ???
+  m_info = info == 0 ? Success : NumericalIssue;
+  m_factorizationIsOk = true;
+}
+
+template<typename MatrixType>
+template<typename Rhs,typename Dest>
+void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
+{
+  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
+
+  const int size = m_matrix.rows();
+  const int rhsCols = b.cols();
+  eigen_assert(size==b.rows());
+
+  m_sluOptions.Trans = NOTRANS;
+  m_sluOptions.Fact = FACTORED;
+  m_sluOptions.IterRefine = NOREFINE;
+
+  m_sluFerr.resize(rhsCols);
+  m_sluBerr.resize(rhsCols);
+  m_sluB = SluMatrix::Map(b.const_cast_derived());
+  m_sluX = SluMatrix::Map(x.derived());
+
+  typename Rhs::PlainObject b_cpy;
+  if(m_sluEqued!='N')
+  {
+    b_cpy = b;
+    m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
+  }
+  
+  int info = 0;
+  RealScalar recip_pivot_growth, rcond;
+
+  StatInit(&m_sluStat);
+  SuperLU_gsisx(&m_sluOptions, &m_sluA,
+                m_q.data(), m_p.data(),
+                &m_sluEtree[0], &m_sluEqued,
+                &m_sluRscale[0], &m_sluCscale[0],
+                &m_sluL, &m_sluU,
+                NULL, 0,
+                &m_sluB, &m_sluX,
+                &recip_pivot_growth, &rcond,
+                &m_sluStat, &info, Scalar());
+  StatFree(&m_sluStat);
+
+  m_info = info==0 ? Success : NumericalIssue;
+}
+#endif
+
+namespace internal {
+  
+template<typename _MatrixType, typename Derived, typename Rhs>
+struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
+  : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
+{
+  typedef SuperLUBase<_MatrixType,Derived> Dec;
+  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    dec().derived()._solve(rhs(),dst);
+  }
+};
+
+template<typename _MatrixType, typename Derived, typename Rhs>
+struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
+  : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
+{
+  typedef SuperLUBase<_MatrixType,Derived> Dec;
+  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    dec().derived()._solve(rhs(),dst);
+  }
+};
+
+}
+
 #endif // EIGEN_SUPERLUSUPPORT_H

unsupported/Eigen/src/SparseExtra/SuperLUSupportLegacy.h

@@ -0,0 +1,404 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_SUPERLUSUPPORT_LEGACY_H
+#define EIGEN_SUPERLUSUPPORT_LEGACY_H
+
+template<typename MatrixType>
+class SparseLU<MatrixType,SuperLULegacy> : public SparseLU<MatrixType>
+{
+  protected:
+    typedef SparseLU<MatrixType> Base;
+    typedef typename Base::Scalar Scalar;
+    typedef typename Base::RealScalar RealScalar;
+    typedef Matrix<Scalar,Dynamic,1> Vector;
+    typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
+    typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
+    typedef SparseMatrix<Scalar,Lower|UnitDiag> LMatrixType;
+    typedef SparseMatrix<Scalar,Upper> UMatrixType;
+    using Base::m_flags;
+    using Base::m_status;
+
+  public:
+
+    SparseLU(int flags = NaturalOrdering)
+      : Base(flags)
+    {
+    }
+
+    SparseLU(const MatrixType& matrix, int flags = NaturalOrdering)
+      : Base(flags)
+    {
+      compute(matrix);
+    }
+
+    ~SparseLU()
+    {
+      Destroy_SuperNode_Matrix(&m_sluL);
+      Destroy_CompCol_Matrix(&m_sluU);
+    }
+
+    inline const LMatrixType& matrixL() const
+    {
+      if (m_extractedDataAreDirty) extractData();
+      return m_l;
+    }
+
+    inline const UMatrixType& matrixU() const
+    {
+      if (m_extractedDataAreDirty) extractData();
+      return m_u;
+    }
+
+    inline const IntColVectorType& permutationP() const
+    {
+      if (m_extractedDataAreDirty) extractData();
+      return m_p;
+    }
+
+    inline const IntRowVectorType& permutationQ() const
+    {
+      if (m_extractedDataAreDirty) extractData();
+      return m_q;
+    }
+
+    Scalar determinant() const;
+
+    template<typename BDerived, typename XDerived>
+    bool solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>* x, const int transposed = SvNoTrans) const;
+
+    void compute(const MatrixType& matrix);
+
+  protected:
+
+    void extractData() const;
+
+  protected:
+    // cached data to reduce reallocation, etc.
+    mutable LMatrixType m_l;
+    mutable UMatrixType m_u;
+    mutable IntColVectorType m_p;
+    mutable IntRowVectorType m_q;
+
+    mutable SparseMatrix<Scalar> m_matrix;
+    mutable SluMatrix m_sluA;
+    mutable SuperMatrix m_sluL, m_sluU;
+    mutable SluMatrix m_sluB, m_sluX;
+    mutable SuperLUStat_t m_sluStat;
+    mutable superlu_options_t m_sluOptions;
+    mutable std::vector<int> m_sluEtree;
+    mutable std::vector<RealScalar> m_sluRscale, m_sluCscale;
+    mutable std::vector<RealScalar> m_sluFerr, m_sluBerr;
+    mutable char m_sluEqued;
+    mutable bool m_extractedDataAreDirty;
+};
+
+template<typename MatrixType>
+void SparseLU<MatrixType,SuperLULegacy>::compute(const MatrixType& a)
+{
+  const int size = a.rows();
+  m_matrix = a;
+
+  set_default_options(&m_sluOptions);
+  m_sluOptions.ColPerm = NATURAL;
+  m_sluOptions.PrintStat = NO;
+  m_sluOptions.ConditionNumber = NO;
+  m_sluOptions.Trans = NOTRANS;
+  // m_sluOptions.Equil = NO;
+
+  switch (Base::orderingMethod())
+  {
+      case NaturalOrdering          : m_sluOptions.ColPerm = NATURAL; break;
+      case MinimumDegree_AT_PLUS_A  : m_sluOptions.ColPerm = MMD_AT_PLUS_A; break;
+      case MinimumDegree_ATA        : m_sluOptions.ColPerm = MMD_ATA; break;
+      case ColApproxMinimumDegree   : m_sluOptions.ColPerm = COLAMD; break;
+      default:
+        //std::cerr << "Eigen: ordering method \"" << Base::orderingMethod() << "\" not supported by the SuperLU backend\n";
+        m_sluOptions.ColPerm = NATURAL;
+  };
+
+  m_sluA = internal::asSluMatrix(m_matrix);
+  memset(&m_sluL,0,sizeof m_sluL);
+  memset(&m_sluU,0,sizeof m_sluU);
+  m_sluEqued = 'N';
+  int info = 0;
+
+  m_p.resize(size);
+  m_q.resize(size);
+  m_sluRscale.resize(size);
+  m_sluCscale.resize(size);
+  m_sluEtree.resize(size);
+
+  RealScalar recip_pivot_gross, rcond;
+  RealScalar ferr, berr;
+
+  // set empty B and X
+  m_sluB.setStorageType(SLU_DN);
+  m_sluB.setScalarType<Scalar>();
+  m_sluB.Mtype = SLU_GE;
+  m_sluB.storage.values = 0;
+  m_sluB.nrow = m_sluB.ncol = 0;
+  m_sluB.storage.lda = size;
+  m_sluX = m_sluB;
+
+  StatInit(&m_sluStat);
+  if (m_flags&IncompleteFactorization)
+  {
+    #ifdef EIGEN_SUPERLU_HAS_ILU
+    ilu_set_default_options(&m_sluOptions);
+
+    // no attempt to preserve column sum
+    m_sluOptions.ILU_MILU = SILU;
+
+    // only basic ILU(k) support -- no direct control over memory consumption
+    // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
+    // and set ILU_FillFactor to max memory growth
+    m_sluOptions.ILU_DropRule = DROP_BASIC;
+    m_sluOptions.ILU_DropTol = Base::m_precision;
+
+    SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
+      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
+      &m_sluL, &m_sluU,
+      NULL, 0,
+      &m_sluB, &m_sluX,
+      &recip_pivot_gross, &rcond,
+      &m_sluStat, &info, Scalar());
+    #else
+    //std::cerr << "Incomplete factorization is only available in SuperLU v4\n";
+    Base::m_succeeded = false;
+    return;
+    #endif
+  }
+  else
+  {
+    SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
+      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
+      &m_sluL, &m_sluU,
+      NULL, 0,
+      &m_sluB, &m_sluX,
+      &recip_pivot_gross, &rcond,
+      &ferr, &berr,
+      &m_sluStat, &info, Scalar());
+  }
+  StatFree(&m_sluStat);
+
+  m_extractedDataAreDirty = true;
+
+  // FIXME how to better check for errors ???
+  Base::m_succeeded = (info == 0);
+}
+
+template<typename MatrixType>
+template<typename BDerived,typename XDerived>
+bool SparseLU<MatrixType,SuperLULegacy>::solve(const MatrixBase<BDerived> &b,
+                        MatrixBase<XDerived> *x, const int transposed) const
+{
+  const int size = m_matrix.rows();
+  const int rhsCols = b.cols();
+  eigen_assert(size==b.rows());
+
+  switch (transposed) {
+      case SvNoTrans    :  m_sluOptions.Trans = NOTRANS; break;
+      case SvTranspose  :  m_sluOptions.Trans = TRANS;   break;
+      case SvAdjoint    :  m_sluOptions.Trans = CONJ;    break;
+      default:
+        //std::cerr << "Eigen: transposition  option \"" << transposed << "\" not supported by the SuperLU backend\n";
+        m_sluOptions.Trans = NOTRANS;
+  }
+
+  m_sluOptions.Fact = FACTORED;
+  m_sluOptions.IterRefine = NOREFINE;
+
+  m_sluFerr.resize(rhsCols);
+  m_sluBerr.resize(rhsCols);
+  m_sluB = SluMatrix::Map(b.const_cast_derived());
+  m_sluX = SluMatrix::Map(x->derived());
+  
+  typename BDerived::PlainObject b_cpy;
+  if(m_sluEqued!='N')
+  {
+    b_cpy = b;
+    m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());  
+  }
+
+  StatInit(&m_sluStat);
+  int info = 0;
+  RealScalar recip_pivot_gross, rcond;
+
+  if (m_flags&IncompleteFactorization)
+  {
+    #ifdef EIGEN_SUPERLU_HAS_ILU
+    SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
+      &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
+      &m_sluL, &m_sluU,
+      NULL, 0,
+      &m_sluB, &m_sluX,
+      &recip_pivot_gross, &rcond,
+      &m_sluStat, &info, Scalar());
+    #else
+    //std::cerr << "Incomplete factorization is only available in SuperLU v4\n";
+    return false;
+    #endif
+  }
+  else
+  {
+    SuperLU_gssvx(
+      &m_sluOptions, &m_sluA,
+      m_q.data(), m_p.data(),
+      &m_sluEtree[0], &m_sluEqued,
+      &m_sluRscale[0], &m_sluCscale[0],
+      &m_sluL, &m_sluU,
+      NULL, 0,
+      &m_sluB, &m_sluX,
+      &recip_pivot_gross, &rcond,
+      &m_sluFerr[0], &m_sluBerr[0],
+      &m_sluStat, &info, Scalar());
+  }
+  StatFree(&m_sluStat);
+
+  // reset to previous state
+  m_sluOptions.Trans = NOTRANS;
+  return info==0;
+}
+
+//
+// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
+//
+//  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+//
+//  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+//  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+//
+template<typename MatrixType>
+void SparseLU<MatrixType,SuperLULegacy>::extractData() const
+{
+  if (m_extractedDataAreDirty)
+  {
+    int         upper;
+    int         fsupc, istart, nsupr;
+    int         lastl = 0, lastu = 0;
+    SCformat    *Lstore = static_cast<SCformat*>(m_sluL.Store);
+    NCformat    *Ustore = static_cast<NCformat*>(m_sluU.Store);
+    Scalar      *SNptr;
+
+    const int size = m_matrix.rows();
+    m_l.resize(size,size);
+    m_l.resizeNonZeros(Lstore->nnz);
+    m_u.resize(size,size);
+    m_u.resizeNonZeros(Ustore->nnz);
+
+    int* Lcol = m_l._outerIndexPtr();
+    int* Lrow = m_l._innerIndexPtr();
+    Scalar* Lval = m_l._valuePtr();
+
+    int* Ucol = m_u._outerIndexPtr();
+    int* Urow = m_u._innerIndexPtr();
+    Scalar* Uval = m_u._valuePtr();
+
+    Ucol[0] = 0;
+    Ucol[0] = 0;
+
+    /* for each supernode */
+    for (int k = 0; k <= Lstore->nsuper; ++k)
+    {
+      fsupc   = L_FST_SUPC(k);
+      istart  = L_SUB_START(fsupc);
+      nsupr   = L_SUB_START(fsupc+1) - istart;
+      upper   = 1;
+
+      /* for each column in the supernode */
+      for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
+      {
+        SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
+
+        /* Extract U */
+        for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
+        {
+          Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
+          /* Matlab doesn't like explicit zero. */
+          if (Uval[lastu] != 0.0)
+            Urow[lastu++] = U_SUB(i);
+        }
+        for (int i = 0; i < upper; ++i)
+        {
+          /* upper triangle in the supernode */
+          Uval[lastu] = SNptr[i];
+          /* Matlab doesn't like explicit zero. */
+          if (Uval[lastu] != 0.0)
+            Urow[lastu++] = L_SUB(istart+i);
+        }
+        Ucol[j+1] = lastu;
+
+        /* Extract L */
+        Lval[lastl] = 1.0; /* unit diagonal */
+        Lrow[lastl++] = L_SUB(istart + upper - 1);
+        for (int i = upper; i < nsupr; ++i)
+        {
+          Lval[lastl] = SNptr[i];
+          /* Matlab doesn't like explicit zero. */
+          if (Lval[lastl] != 0.0)
+            Lrow[lastl++] = L_SUB(istart+i);
+        }
+        Lcol[j+1] = lastl;
+
+        ++upper;
+      } /* for j ... */
+
+    } /* for k ... */
+
+    // squeeze the matrices :
+    m_l.resizeNonZeros(lastl);
+    m_u.resizeNonZeros(lastu);
+
+    m_extractedDataAreDirty = false;
+  }
+}
+
+template<typename MatrixType>
+typename SparseLU<MatrixType,SuperLULegacy>::Scalar SparseLU<MatrixType,SuperLULegacy>::determinant() const
+{
+  assert((!NumTraits<Scalar>::IsComplex) && "This function is not implemented for complex yet");
+  if (m_extractedDataAreDirty)
+    extractData();
+
+  // TODO this code could be moved to the default/base backend
+  // FIXME perhaps we have to take into account the scale factors m_sluRscale and m_sluCscale ???
+  Scalar det = Scalar(1);
+  for (int j=0; j<m_u.cols(); ++j)
+  {
+    if (m_u._outerIndexPtr()[j+1]-m_u._outerIndexPtr()[j] > 0)
+    {
+      int lastId = m_u._outerIndexPtr()[j+1]-1;
+      eigen_assert(m_u._innerIndexPtr()[lastId]<=j);
+      if (m_u._innerIndexPtr()[lastId]==j)
+      {
+        det *= m_u._valuePtr()[lastId];
+      }
+    }
+//       std::cout << m_sluRscale[j] << " " << m_sluCscale[j] << "   \n";
+  }
+  return det;
+}
+
+#endif // EIGEN_SUPERLUSUPPORT_LEGACY_H

unsupported/test/sparse_lu.cpp

@@ -76,27 +76,62 @@ template<typename Scalar> void sparse_lu(int rows, int cols)
     #endif
     
     #ifdef EIGEN_SUPERLU_SUPPORT
+    // legacy, deprecated API
     {
       x.setZero();
-      SparseLU<SparseMatrix<Scalar>,SuperLU> slu(m2);
+      SparseLU<SparseMatrix<Scalar>,SuperLULegacy> slu(m2);
       if (slu.succeeded())
       {
+        DenseVector oldb = b;
         if (slu.solve(b,&x)) {
           VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU");
         }
+        else
+          std::cerr << "super lu solving failed\n";
+        VERIFY(oldb.isApprox(b) && "the rhs should not be modified!");
+        
         // std::cerr << refDet << " == " << slu.determinant() << "\n";
         if (slu.solve(b, &x, SvTranspose)) {
           VERIFY(b.isApprox(m2.transpose() * x, test_precision<Scalar>()));
         }
+        else
+          std::cerr << "super lu solving failed\n";
 
         if (slu.solve(b, &x, SvAdjoint)) {
          VERIFY(b.isApprox(m2.adjoint() * x, test_precision<Scalar>()));
         }
+        else
+          std::cerr << "super lu solving failed\n";
 
         if (!NumTraits<Scalar>::IsComplex) {
           VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
         }
       }
+      else
+        std::cerr << "super lu factorize failed\n";
+    }
+    
+    // New API
+    {
+      x.setZero();
+      SuperLU<SparseMatrix<Scalar> > slu(m2);
+      if (slu.info() == Success)
+      {
+        DenseVector oldb = b;
+        x = slu.solve(b);
+        VERIFY(oldb.isApprox(b) && "the rhs should not be modified!");
+        if (slu.info() == Success) {
+          VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "SuperLU");
+        }
+        else
+          std::cerr << "super lu solving failed\n";
+
+        if (!NumTraits<Scalar>::IsComplex) {
+          VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
+        }
+      }
+      else
+        std::cerr << "super lu factorize failed\n";
     }
     #endif