SuperLUSupport.h

Go to the documentation of this file.

// 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
 
// IWYU pragma: private
#include "./InternalHeaderCheck.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;
 
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);
}
 
template <typename Scalar, int Flags, typename Index>
Map<SparseMatrix<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 Map<SparseMatrix<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
 
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(); }
 
  inline superlu_options_t &options() { return m_sluOptions; }
 
  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "Decomposition is not initialized.");
    return m_info;
  }
 
  void compute(const MatrixType &matrix) {
    derived().analyzePattern(matrix);
    derived().factorize(matrix);
  }
 
  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 &) {}
};
 
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() {}
 
  void analyzePattern(const MatrixType &matrix) {
    m_info = InvalidInput;
    m_isInitialized = false;
    Base::analyzePattern(matrix);
  }
 
  void factorize(const MatrixType &matrix);
 
  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_l;
  using Base::m_matrix;
  using Base::m_p;
  using Base::m_q;
  using Base::m_sluA;
  using Base::m_sluB;
  using Base::m_sluBerr;
  using Base::m_sluCscale;
  using Base::m_sluEqued;
  using Base::m_sluEtree;
  using Base::m_sluFerr;
  using Base::m_sluL;
  using Base::m_sluOptions;
  using Base::m_sluRscale;
  using Base::m_sluStat;
  using Base::m_sluU;
  using Base::m_sluX;
  using Base::m_u;
 
  using Base::m_analysisIsOk;
  using Base::m_extractedDataAreDirty;
  using Base::m_factorizationIsOk;
  using Base::m_info;
  using Base::m_isInitialized;
 
  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
 
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() {}
 
  void analyzePattern(const MatrixType &matrix) { Base::analyzePattern(matrix); }
 
  void factorize(const MatrixType &matrix);
 
#ifndef EIGEN_PARSED_BY_DOXYGEN
  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_l;
  using Base::m_matrix;
  using Base::m_p;
  using Base::m_q;
  using Base::m_sluA;
  using Base::m_sluB;
  using Base::m_sluBerr;
  using Base::m_sluCscale;
  using Base::m_sluEqued;
  using Base::m_sluEtree;
  using Base::m_sluFerr;
  using Base::m_sluL;
  using Base::m_sluOptions;
  using Base::m_sluRscale;
  using Base::m_sluStat;
  using Base::m_sluU;
  using Base::m_sluX;
  using Base::m_u;
 
  using Base::m_analysisIsOk;
  using Base::m_extractedDataAreDirty;
  using Base::m_factorizationIsOk;
  using Base::m_info;
  using Base::m_isInitialized;
 
  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