text: // This file is part of Eigen, a lightweight C++ template library
  // for linear algebra.
  //
  // Copyright (C) 2008-2015 Gael Guennebaud <gael.guennebaud@inria.fr>
  //
  // This Source Code Form is subject to the terms of the Mozilla
  // Public License v. 2.0. If a copy of the MPL was not distributed
  // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
  
  #ifndef EIGEN_SUPERLUSUPPORT_H
  #define EIGEN_SUPERLUSUPPORT_H
  
  namespace Eigen {
  
  #if defined(SUPERLU_MAJOR_VERSION) && (SUPERLU_MAJOR_VERSION >= 5)
  #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE)		\
      extern "C" {                                                                                          \
        extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,                  \
                                  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,           \
                                  void *, int, SuperMatrix *, SuperMatrix *,                                \
                                  FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *,                       \
                                  GlobalLU_t *, mem_usage_t *, SuperLUStat_t *, int *);                     \
      }                                                                                                     \
      inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,                                \
           int *perm_c, int *perm_r, int *etree, char *equed,                                               \
           FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                                      \
           SuperMatrix *U, void *work, int lwork,                                                           \
           SuperMatrix *B, SuperMatrix *X,                                                                  \
           FLOATTYPE *recip_pivot_growth,                                                                   \
           FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr,                                              \
           SuperLUStat_t *stats, int *info, KEYTYPE) {                                                      \
      mem_usage_t mem_usage;                                                                                \
      GlobalLU_t gLU;                                                                                       \
      PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L,                                      \
           U, work, lwork, B, X, recip_pivot_growth, rcond,                                                 \
           ferr, berr, &gLU, &mem_usage, stats, info);                                                      \
      return mem_usage.for_lu; /* bytes used by the factor storage */                                       \
    }
  #else // version < 5.0
  #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE)		\
      extern "C" {                                                                                          \
        extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,                  \
                                  char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,           \
                                  void *, int, SuperMatrix *, SuperMatrix *,                                \
                                  FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *,                       \
                                  mem_usage_t *, SuperLUStat_t *, int *);                                   \
      }                                                                                                     \
      inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,                                \
           int *perm_c, int *perm_r, int *etree, char *equed,                                               \
           FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                                      \
           SuperMatrix *U, void *work, int lwork,                                                           \
           SuperMatrix *B, SuperMatrix *X,                                                                  \
           FLOATTYPE *recip_pivot_growth,                                                                   \
           FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr,                                              \
           SuperLUStat_t *stats, int *info, KEYTYPE) {                                                      \
      mem_usage_t mem_usage;                                                                                \
      PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L,                                      \
           U, work, lwork, B, X, recip_pivot_growth, rcond,                                                 \
           ferr, berr, &mem_usage, stats, info);                                                            \
      return mem_usage.for_lu; /* bytes used by the factor storage */                                       \
    }
  #endif
  
  DECL_GSSVX(s,float,float)
  DECL_GSSVX(c,float,std::complex<float>)
  DECL_GSSVX(d,double,double)
  DECL_GSSVX(z,double,std::complex<double>)
  
  #ifdef MILU_ALPHA
  #define EIGEN_SUPERLU_HAS_ILU
  #endif
  
  #ifdef EIGEN_SUPERLU_HAS_ILU
  
  // similarly for the incomplete factorization using gsisx
  #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE)                                                    \
      extern "C" {                                                                                \
        extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *,        \
                           char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,        \
                           void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *,   \
                           mem_usage_t *, SuperLUStat_t *, int *);                        \
      }                                                                                           \
      inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A,                      \
           int *perm_c, int *perm_r, int *etree, char *equed,                                     \
           FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                            \
           SuperMatrix *U, void *work, int lwork,                                                 \
           SuperMatrix *B, SuperMatrix *X,                                                        \
           FLOATTYPE *recip_pivot_growth,                                                         \
           FLOATTYPE *rcond,                                                                      \
           SuperLUStat_t *stats, int *info, KEYTYPE) {                                            \
      mem_usage_t mem_usage;                                                              \
      PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L,                            \
           U, work, lwork, B, X, recip_pivot_growth, rcond,                                       \
           &mem_usage, stats, info);                                                              \
      return mem_usage.for_lu; /* bytes used by the factor storage */                             \
    }
  
  DECL_GSISX(s,float,float)
  DECL_GSISX(c,float,std::complex<float>)
  DECL_GSISX(d,double,double)
  DECL_GSISX(z,double,std::complex<double>)
  
  #endif
  
  template<typename MatrixType>
  struct SluMatrixMapHelper;
  
  /** \internal
    *
    * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
    * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
    *
    * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
    */
  struct SluMatrix : SuperMatrix
  {
    SluMatrix()
    {
      Store = &storage;
    }
  
    SluMatrix(const SluMatrix& other)
      : SuperMatrix(other)
    {
      Store = &storage;
      storage = other.storage;
    }
  
    SluMatrix& operator=(const SluMatrix& other)
    {
      SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
      Store = &storage;
      storage = other.storage;
      return *this;
    }
  
    struct
    {
      union {int nnz;int lda;};
      void *values;
      int *innerInd;
      int *outerInd;
    } storage;
  
    void setStorageType(Stype_t t)
    {
      Stype = t;
      if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
        Store = &storage;
      else
      {
        eigen_assert(false && "storage type not supported");
        Store = 0;
      }
    }
  
    template<typename Scalar>
    void setScalarType()
    {
      if (internal::is_same<Scalar,float>::value)
        Dtype = SLU_S;
      else if (internal::is_same<Scalar,double>::value)
        Dtype = SLU_D;
      else if (internal::is_same<Scalar,std::complex<float> >::value)
        Dtype = SLU_C;
      else if (internal::is_same<Scalar,std::complex<double> >::value)
        Dtype = SLU_Z;
      else
      {
        eigen_assert(false && "Scalar type not supported by SuperLU");
      }
    }
  
    template<typename MatrixType>
    static SluMatrix Map(MatrixBase<MatrixType>& _mat)
    {
      MatrixType& mat(_mat.derived());
      eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
      SluMatrix res;
      res.setStorageType(SLU_DN);
      res.setScalarType<typename MatrixType::Scalar>();
      res.Mtype     = SLU_GE;
  
      res.nrow      = internal::convert_index<int>(mat.rows());
      res.ncol      = internal::convert_index<int>(mat.cols());
  
      res.storage.lda       = internal::convert_index<int>(MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride());
      res.storage.values    = (void*)(mat.data());
      return res;
    }
  
    template<typename MatrixType>
    static SluMatrix Map(SparseMatrixBase<MatrixType>& a_mat)
    {
      MatrixType &mat(a_mat.derived());
      SluMatrix res;
      if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
      {
        res.setStorageType(SLU_NR);
        res.nrow      = internal::convert_index<int>(mat.cols());
        res.ncol      = internal::convert_index<int>(mat.rows());
      }
      else
      {
        res.setStorageType(SLU_NC);
        res.nrow      = internal::convert_index<int>(mat.rows());
        res.ncol      = internal::convert_index<int>(mat.cols());
      }
  
      res.Mtype       = SLU_GE;
  
      res.storage.nnz       = internal::convert_index<int>(mat.nonZeros());
      res.storage.values    = mat.valuePtr();
      res.storage.innerInd  = mat.innerIndexPtr();
      res.storage.outerInd  = mat.outerIndexPtr();
  
      res.setScalarType<typename MatrixType::Scalar>();
  
      // FIXME the following is not very accurate
      if (int(MatrixType::Flags) & int(Upper))
        res.Mtype = SLU_TRU;
      if (int(MatrixType::Flags) & int(Lower))
        res.Mtype = SLU_TRL;
  
      eigen_assert(((int(MatrixType::Flags) & int(SelfAdjoint))==0) && "SelfAdjoint matrix shape not supported by SuperLU");
  
      return res;
    }
  };
  
  template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
  struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
  {
    typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
    static void run(MatrixType& mat, SluMatrix& res)
    {
      eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
      res.setStorageType(SLU_DN);
      res.setScalarType<Scalar>();
      res.Mtype     = SLU_GE;
  
      res.nrow      = mat.rows();
      res.ncol      = mat.cols();
  
      res.storage.lda       = mat.outerStride();
      res.storage.values    = mat.data();
    }
  };
  
  template<typename Derived>
  struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
  {
    typedef Derived MatrixType;
    static void run(MatrixType& mat, SluMatrix& res)
    {
      if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
      {
        res.setStorageType(SLU_NR);
        res.nrow      = mat.cols();
        res.ncol      = mat.rows();
      }
      else
      {
        res.setStorageType(SLU_NC);
        res.nrow      = mat.rows();
        res.ncol      = mat.cols();
      }
  
      res.Mtype       = SLU_GE;
  
      res.storage.nnz       = mat.nonZeros();
      res.storage.values    = mat.valuePtr();
      res.storage.innerInd  = mat.innerIndexPtr();
      res.storage.outerInd  = mat.outerIndexPtr();
  
      res.setScalarType<typename MatrixType::Scalar>();
  
      // FIXME the following is not very accurate
      if (MatrixType::Flags & Upper)
        res.Mtype = SLU_TRU;
      if (MatrixType::Flags & Lower)
        res.Mtype = SLU_TRL;
  
      eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
    }
  };
  
  namespace internal {
  
  template<typename MatrixType>
  SluMatrix asSluMatrix(MatrixType& mat)
  {
    return SluMatrix::Map(mat);
  }
  
  /** View a Super LU matrix as an Eigen expression */
  template<typename Scalar, int Flags, typename Index>
  MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
  {
    eigen_assert(((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR)
           || ((Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC));
  
    Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
  
    return MappedSparseMatrix<Scalar,Flags,Index>(
      sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
      sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
  }
  
  } // end namespace internal
  
  /** \ingroup SuperLUSupport_Module
    * \class SuperLUBase
    * \brief The base class for the direct and incomplete LU factorization of SuperLU
    */
  template<typename _MatrixType, typename Derived>
  class SuperLUBase : public SparseSolverBase<Derived>
  {
    protected:
      typedef SparseSolverBase<Derived> Base;
      using Base::derived;
      using Base::m_isInitialized;
    public:
      typedef _MatrixType MatrixType;
      typedef typename MatrixType::Scalar Scalar;
      typedef typename MatrixType::RealScalar RealScalar;
      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 Map<PermutationMatrix<Dynamic,Dynamic,int> > PermutationMap;
      typedef SparseMatrix<Scalar> LUMatrixType;
      enum {
        ColsAtCompileTime = MatrixType::ColsAtCompileTime,
        MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
      };
  
    public:
  
      SuperLUBase() {}
  
      ~SuperLUBase()
      {
        clearFactors();
      }
      
      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 successful,
        *          \c NumericalIssue if the matrix.appears to be negative.
        */
      ComputationInfo info() const
      {
        eigen_assert(m_isInitialized && "Decomposition is not initialized.");
        return m_info;
      }
  
      /** Computes the sparse Cholesky decomposition of \a matrix */
      void compute(const MatrixType& matrix)
      {
        derived().analyzePattern(matrix);
        derived().factorize(matrix);
      }
  
      /** 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*/)
      {
        m_isInitialized = true;
        m_info = Success;
        m_analysisIsOk = true;
        m_factorizationIsOk = false;
      }
      
      template<typename Stream>
      void dumpMemory(Stream& /*s*/)
      {}
      
    protected:
      
      void initFactorization(const MatrixType& a)
      {
        set_default_options(&this->m_sluOptions);
        
        const Index size = a.rows();
        m_matrix = a;
  
        m_sluA = internal::asSluMatrix(m_matrix);
        clearFactors();
  
        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    = internal::convert_index<int>(size);
        m_sluX                = m_sluB;
        
        m_extractedDataAreDirty = true;
      }
      
      void init()
      {
        m_info = InvalidInput;
        m_isInitialized = false;
        m_sluL.Store = 0;
        m_sluU.Store = 0;
      }
      
      void extractData() const;
  
      void clearFactors()
      {
        if(m_sluL.Store)
          Destroy_SuperNode_Matrix(&m_sluL);
        if(m_sluU.Store)
          Destroy_CompCol_Matrix(&m_sluU);
  
        m_sluL.Store = 0;
        m_sluU.Store = 0;
  
        memset(&m_sluL,0,sizeof m_sluL);
        memset(&m_sluU,0,sizeof m_sluU);
      }
  
      // cached data to reduce reallocation, etc.
      mutable LUMatrixType m_l;
      mutable LUMatrixType m_u;
      mutable IntColVectorType m_p;
      mutable IntRowVectorType m_q;
  
      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;
      mutable SuperLUStat_t m_sluStat;
      mutable superlu_options_t m_sluOptions;
      mutable std::vector<int> m_sluEtree;
      mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
      mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
      mutable char m_sluEqued;
  
      mutable ComputationInfo m_info;
      int m_factorizationIsOk;
      int m_analysisIsOk;
      mutable bool m_extractedDataAreDirty;
      
    private:
      SuperLUBase(SuperLUBase& ) { }
  };
  
  
  /** \ingroup SuperLUSupport_Module
    * \class SuperLU
    * \brief A sparse direct LU factorization and solver based on the SuperLU library
    *
    * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
    * using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices
    * X and B can be either dense or sparse.
    *
    * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
    *
    * \warning This class is only for the 4.x versions of SuperLU. The 3.x and 5.x versions are not supported.
    *
    * \implsparsesolverconcept
    *
    * \sa \ref TutorialSparseSolverConcept, class SparseLU
    */
  template<typename _MatrixType>
  class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
  {
    public:
      typedef SuperLUBase<_MatrixType,SuperLU> Base;
      typedef _MatrixType MatrixType;
      typedef typename Base::Scalar Scalar;
      typedef typename Base::RealScalar RealScalar;
      typedef typename Base::StorageIndex StorageIndex;
      typedef typename Base::IntRowVectorType IntRowVectorType;
      typedef typename Base::IntColVectorType IntColVectorType;   
      typedef typename Base::PermutationMap PermutationMap;
      typedef typename Base::LUMatrixType LUMatrixType;
      typedef TriangularView<LUMatrixType, Lower|UnitDiag>  LMatrixType;
      typedef TriangularView<LUMatrixType,  Upper>          UMatrixType;
  
    public:
      using Base::_solve_impl;
  
      SuperLU() : Base() { init(); }
  
      explicit SuperLU(const MatrixType& matrix) : Base()
      {
        init();
        Base::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)
      {
        m_info = InvalidInput;
        m_isInitialized = false;
        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);
      
      /** \internal */
      template<typename Rhs,typename Dest>
      void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
      
      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;
      }
      
      
    private:
      SuperLU(SuperLU& ) { }
  };
  
  template<typename MatrixType>
  void SuperLU<MatrixType>::factorize(const MatrixType& a)
  {
    eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
    if(!m_analysisIsOk)
    {
      m_info = InvalidInput;
      return;
    }
    
    this->initFactorization(a);
    
    m_sluOptions.ColPerm = COLAMD;
    int info = 0;
    RealScalar recip_pivot_growth, rcond;
    RealScalar ferr, berr;
  
    StatInit(&m_sluStat);
    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_growth, &rcond,
                  &ferr, &berr,
                  &m_sluStat, &info, Scalar());
    StatFree(&m_sluStat);
  
    m_extractedDataAreDirty = true;
  
    // 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 SuperLU<MatrixType>::_solve_impl(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 Index rhsCols = b.cols();
    eigen_assert(m_matrix.rows()==b.rows());
  
    m_sluOptions.Trans = NOTRANS;
    m_sluOptions.Fact = FACTORED;
    m_sluOptions.IterRefine = NOREFINE;
    
  
    m_sluFerr.resize(rhsCols);
    m_sluBerr.resize(rhsCols);
    
    Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b);
    Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(x);
    
    m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
    m_sluX = SluMatrix::Map(x_ref.const_cast_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_growth, rcond;
    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_growth, &rcond,
                  &m_sluFerr[0], &m_sluBerr[0],
                  &m_sluStat, &info, Scalar());
    StatFree(&m_sluStat);
    
    if(x.derived().data() != x_ref.data())
      x = x_ref;
    
    m_info = info==0 ? Success : NumericalIssue;
  }
  
  // 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, typename Derived>
  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;
      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 Index 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 SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
  {
    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)
      this->extractData();
  
    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];
      }
    }
    if(PermutationMap(m_p.data(),m_p.size()).determinant()*PermutationMap(m_q.data(),m_q.size()).determinant()<0)
      det = -det;
    if(m_sluEqued!='N')
      return det/m_sluRscale.prod()/m_sluCscale.prod();
    else
      return det;
  }
  
  #ifdef EIGEN_PARSED_BY_DOXYGEN
  #define EIGEN_SUPERLU_HAS_ILU
  #endif
  
  #ifdef EIGEN_SUPERLU_HAS_ILU
  
  /** \ingroup SuperLUSupport_Module
    * \class SuperILU
    * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library
    *
    * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization
    * using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers.
    *
    * \warning This class is only for the 4.x versions of SuperLU. The 3.x and 5.x versions are not supported.
    *
    * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
    *
    * \implsparsesolverconcept
    *
    * \sa \ref TutorialSparseSolverConcept, class IncompleteLUT, class ConjugateGradient, class BiCGSTAB
    */
  
  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;
  
    public:
      using Base::_solve_impl;
  
      SuperILU() : Base() { init(); }
  
      SuperILU(const MatrixType& matrix) : Base()
      {
        init();
        Base::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_impl(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;
      }
      
    private:
      SuperILU(SuperILU& ) { }
  };
  
  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;
  }
  
  #ifndef EIGEN_PARSED_BY_DOXYGEN
  template<typename MatrixType>
  template<typename Rhs,typename Dest>
  void SuperILU<MatrixType>::_solve_impl(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 rhsCols = b.cols();
    eigen_assert(m_matrix.rows()==b.rows());
  
    m_sluOptions.Trans = NOTRANS;
    m_sluOptions.Fact = FACTORED;
    m_sluOptions.IterRefine = NOREFINE;
  
    m_sluFerr.resize(rhsCols);
    m_sluBerr.resize(rhsCols);
    
    Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b);
    Ref<const Matrix<typename Dest::Scalar,Dynamic,Dynamic,ColMajor> > x_ref(x);
    
    m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
    m_sluX = SluMatrix::Map(x_ref.const_cast_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);
    
    if(x.derived().data() != x_ref.data())
      x = x_ref;
  
    m_info = info==0 ? Success : NumericalIssue;
  }
  #endif
  
  #endif
  
  } // end namespace Eigen
  
  #endif // EIGEN_SUPERLUSUPPORT_H
  