Home | History | Annotate | Download | only in SuperLUSupport
      1 // This file is part of Eigen, a lightweight C++ template library
      2 // for linear algebra.
      3 //
      4 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud (at) inria.fr>
      5 //
      6 // This Source Code Form is subject to the terms of the Mozilla
      7 // Public License v. 2.0. If a copy of the MPL was not distributed
      8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
      9 
     10 #ifndef EIGEN_SUPERLUSUPPORT_H
     11 #define EIGEN_SUPERLUSUPPORT_H
     12 
     13 namespace Eigen {
     14 
     15 #define DECL_GSSVX(PREFIX,FLOATTYPE,KEYTYPE)		\
     16     extern "C" {                                                                                          \
     17       typedef struct { FLOATTYPE for_lu; FLOATTYPE total_needed; int expansions; } PREFIX##mem_usage_t;   \
     18       extern void PREFIX##gssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,                  \
     19                                 char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,           \
     20                                 void *, int, SuperMatrix *, SuperMatrix *,                                \
     21                                 FLOATTYPE *, FLOATTYPE *, FLOATTYPE *, FLOATTYPE *,                       \
     22                                 PREFIX##mem_usage_t *, SuperLUStat_t *, int *);                           \
     23     }                                                                                                     \
     24     inline float SuperLU_gssvx(superlu_options_t *options, SuperMatrix *A,                                \
     25          int *perm_c, int *perm_r, int *etree, char *equed,                                               \
     26          FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                                      \
     27          SuperMatrix *U, void *work, int lwork,                                                           \
     28          SuperMatrix *B, SuperMatrix *X,                                                                  \
     29          FLOATTYPE *recip_pivot_growth,                                                                   \
     30          FLOATTYPE *rcond, FLOATTYPE *ferr, FLOATTYPE *berr,                                              \
     31          SuperLUStat_t *stats, int *info, KEYTYPE) {                                                      \
     32     PREFIX##mem_usage_t mem_usage;                                                                        \
     33     PREFIX##gssvx(options, A, perm_c, perm_r, etree, equed, R, C, L,                                      \
     34          U, work, lwork, B, X, recip_pivot_growth, rcond,                                                 \
     35          ferr, berr, &mem_usage, stats, info);                                                            \
     36     return mem_usage.for_lu; /* bytes used by the factor storage */                                       \
     37   }
     38 
     39 DECL_GSSVX(s,float,float)
     40 DECL_GSSVX(c,float,std::complex<float>)
     41 DECL_GSSVX(d,double,double)
     42 DECL_GSSVX(z,double,std::complex<double>)
     43 
     44 #ifdef MILU_ALPHA
     45 #define EIGEN_SUPERLU_HAS_ILU
     46 #endif
     47 
     48 #ifdef EIGEN_SUPERLU_HAS_ILU
     49 
     50 // similarly for the incomplete factorization using gsisx
     51 #define DECL_GSISX(PREFIX,FLOATTYPE,KEYTYPE)                                                    \
     52     extern "C" {                                                                                \
     53       extern void PREFIX##gsisx(superlu_options_t *, SuperMatrix *, int *, int *, int *,        \
     54                          char *, FLOATTYPE *, FLOATTYPE *, SuperMatrix *, SuperMatrix *,        \
     55                          void *, int, SuperMatrix *, SuperMatrix *, FLOATTYPE *, FLOATTYPE *,   \
     56                          PREFIX##mem_usage_t *, SuperLUStat_t *, int *);                        \
     57     }                                                                                           \
     58     inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A,                      \
     59          int *perm_c, int *perm_r, int *etree, char *equed,                                     \
     60          FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L,                                            \
     61          SuperMatrix *U, void *work, int lwork,                                                 \
     62          SuperMatrix *B, SuperMatrix *X,                                                        \
     63          FLOATTYPE *recip_pivot_growth,                                                         \
     64          FLOATTYPE *rcond,                                                                      \
     65          SuperLUStat_t *stats, int *info, KEYTYPE) {                                            \
     66     PREFIX##mem_usage_t mem_usage;                                                              \
     67     PREFIX##gsisx(options, A, perm_c, perm_r, etree, equed, R, C, L,                            \
     68          U, work, lwork, B, X, recip_pivot_growth, rcond,                                       \
     69          &mem_usage, stats, info);                                                              \
     70     return mem_usage.for_lu; /* bytes used by the factor storage */                             \
     71   }
     72 
     73 DECL_GSISX(s,float,float)
     74 DECL_GSISX(c,float,std::complex<float>)
     75 DECL_GSISX(d,double,double)
     76 DECL_GSISX(z,double,std::complex<double>)
     77 
     78 #endif
     79 
     80 template<typename MatrixType>
     81 struct SluMatrixMapHelper;
     82 
     83 /** \internal
     84   *
     85   * A wrapper class for SuperLU matrices. It supports only compressed sparse matrices
     86   * and dense matrices. Supernodal and other fancy format are not supported by this wrapper.
     87   *
     88   * This wrapper class mainly aims to avoids the need of dynamic allocation of the storage structure.
     89   */
     90 struct SluMatrix : SuperMatrix
     91 {
     92   SluMatrix()
     93   {
     94     Store = &storage;
     95   }
     96 
     97   SluMatrix(const SluMatrix& other)
     98     : SuperMatrix(other)
     99   {
    100     Store = &storage;
    101     storage = other.storage;
    102   }
    103 
    104   SluMatrix& operator=(const SluMatrix& other)
    105   {
    106     SuperMatrix::operator=(static_cast<const SuperMatrix&>(other));
    107     Store = &storage;
    108     storage = other.storage;
    109     return *this;
    110   }
    111 
    112   struct
    113   {
    114     union {int nnz;int lda;};
    115     void *values;
    116     int *innerInd;
    117     int *outerInd;
    118   } storage;
    119 
    120   void setStorageType(Stype_t t)
    121   {
    122     Stype = t;
    123     if (t==SLU_NC || t==SLU_NR || t==SLU_DN)
    124       Store = &storage;
    125     else
    126     {
    127       eigen_assert(false && "storage type not supported");
    128       Store = 0;
    129     }
    130   }
    131 
    132   template<typename Scalar>
    133   void setScalarType()
    134   {
    135     if (internal::is_same<Scalar,float>::value)
    136       Dtype = SLU_S;
    137     else if (internal::is_same<Scalar,double>::value)
    138       Dtype = SLU_D;
    139     else if (internal::is_same<Scalar,std::complex<float> >::value)
    140       Dtype = SLU_C;
    141     else if (internal::is_same<Scalar,std::complex<double> >::value)
    142       Dtype = SLU_Z;
    143     else
    144     {
    145       eigen_assert(false && "Scalar type not supported by SuperLU");
    146     }
    147   }
    148 
    149   template<typename MatrixType>
    150   static SluMatrix Map(MatrixBase<MatrixType>& _mat)
    151   {
    152     MatrixType& mat(_mat.derived());
    153     eigen_assert( ((MatrixType::Flags&RowMajorBit)!=RowMajorBit) && "row-major dense matrices are not supported by SuperLU");
    154     SluMatrix res;
    155     res.setStorageType(SLU_DN);
    156     res.setScalarType<typename MatrixType::Scalar>();
    157     res.Mtype     = SLU_GE;
    158 
    159     res.nrow      = mat.rows();
    160     res.ncol      = mat.cols();
    161 
    162     res.storage.lda       = MatrixType::IsVectorAtCompileTime ? mat.size() : mat.outerStride();
    163     res.storage.values    = mat.data();
    164     return res;
    165   }
    166 
    167   template<typename MatrixType>
    168   static SluMatrix Map(SparseMatrixBase<MatrixType>& mat)
    169   {
    170     SluMatrix res;
    171     if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
    172     {
    173       res.setStorageType(SLU_NR);
    174       res.nrow      = mat.cols();
    175       res.ncol      = mat.rows();
    176     }
    177     else
    178     {
    179       res.setStorageType(SLU_NC);
    180       res.nrow      = mat.rows();
    181       res.ncol      = mat.cols();
    182     }
    183 
    184     res.Mtype       = SLU_GE;
    185 
    186     res.storage.nnz       = mat.nonZeros();
    187     res.storage.values    = mat.derived().valuePtr();
    188     res.storage.innerInd  = mat.derived().innerIndexPtr();
    189     res.storage.outerInd  = mat.derived().outerIndexPtr();
    190 
    191     res.setScalarType<typename MatrixType::Scalar>();
    192 
    193     // FIXME the following is not very accurate
    194     if (MatrixType::Flags & Upper)
    195       res.Mtype = SLU_TRU;
    196     if (MatrixType::Flags & Lower)
    197       res.Mtype = SLU_TRL;
    198 
    199     eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
    200 
    201     return res;
    202   }
    203 };
    204 
    205 template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
    206 struct SluMatrixMapHelper<Matrix<Scalar,Rows,Cols,Options,MRows,MCols> >
    207 {
    208   typedef Matrix<Scalar,Rows,Cols,Options,MRows,MCols> MatrixType;
    209   static void run(MatrixType& mat, SluMatrix& res)
    210   {
    211     eigen_assert( ((Options&RowMajor)!=RowMajor) && "row-major dense matrices is not supported by SuperLU");
    212     res.setStorageType(SLU_DN);
    213     res.setScalarType<Scalar>();
    214     res.Mtype     = SLU_GE;
    215 
    216     res.nrow      = mat.rows();
    217     res.ncol      = mat.cols();
    218 
    219     res.storage.lda       = mat.outerStride();
    220     res.storage.values    = mat.data();
    221   }
    222 };
    223 
    224 template<typename Derived>
    225 struct SluMatrixMapHelper<SparseMatrixBase<Derived> >
    226 {
    227   typedef Derived MatrixType;
    228   static void run(MatrixType& mat, SluMatrix& res)
    229   {
    230     if ((MatrixType::Flags&RowMajorBit)==RowMajorBit)
    231     {
    232       res.setStorageType(SLU_NR);
    233       res.nrow      = mat.cols();
    234       res.ncol      = mat.rows();
    235     }
    236     else
    237     {
    238       res.setStorageType(SLU_NC);
    239       res.nrow      = mat.rows();
    240       res.ncol      = mat.cols();
    241     }
    242 
    243     res.Mtype       = SLU_GE;
    244 
    245     res.storage.nnz       = mat.nonZeros();
    246     res.storage.values    = mat.valuePtr();
    247     res.storage.innerInd  = mat.innerIndexPtr();
    248     res.storage.outerInd  = mat.outerIndexPtr();
    249 
    250     res.setScalarType<typename MatrixType::Scalar>();
    251 
    252     // FIXME the following is not very accurate
    253     if (MatrixType::Flags & Upper)
    254       res.Mtype = SLU_TRU;
    255     if (MatrixType::Flags & Lower)
    256       res.Mtype = SLU_TRL;
    257 
    258     eigen_assert(((MatrixType::Flags & SelfAdjoint)==0) && "SelfAdjoint matrix shape not supported by SuperLU");
    259   }
    260 };
    261 
    262 namespace internal {
    263 
    264 template<typename MatrixType>
    265 SluMatrix asSluMatrix(MatrixType& mat)
    266 {
    267   return SluMatrix::Map(mat);
    268 }
    269 
    270 /** View a Super LU matrix as an Eigen expression */
    271 template<typename Scalar, int Flags, typename Index>
    272 MappedSparseMatrix<Scalar,Flags,Index> map_superlu(SluMatrix& sluMat)
    273 {
    274   eigen_assert((Flags&RowMajor)==RowMajor && sluMat.Stype == SLU_NR
    275          || (Flags&ColMajor)==ColMajor && sluMat.Stype == SLU_NC);
    276 
    277   Index outerSize = (Flags&RowMajor)==RowMajor ? sluMat.ncol : sluMat.nrow;
    278 
    279   return MappedSparseMatrix<Scalar,Flags,Index>(
    280     sluMat.nrow, sluMat.ncol, sluMat.storage.outerInd[outerSize],
    281     sluMat.storage.outerInd, sluMat.storage.innerInd, reinterpret_cast<Scalar*>(sluMat.storage.values) );
    282 }
    283 
    284 } // end namespace internal
    285 
    286 /** \ingroup SuperLUSupport_Module
    287   * \class SuperLUBase
    288   * \brief The base class for the direct and incomplete LU factorization of SuperLU
    289   */
    290 template<typename _MatrixType, typename Derived>
    291 class SuperLUBase : internal::noncopyable
    292 {
    293   public:
    294     typedef _MatrixType MatrixType;
    295     typedef typename MatrixType::Scalar Scalar;
    296     typedef typename MatrixType::RealScalar RealScalar;
    297     typedef typename MatrixType::Index Index;
    298     typedef Matrix<Scalar,Dynamic,1> Vector;
    299     typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
    300     typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
    301     typedef SparseMatrix<Scalar> LUMatrixType;
    302 
    303   public:
    304 
    305     SuperLUBase() {}
    306 
    307     ~SuperLUBase()
    308     {
    309       clearFactors();
    310     }
    311 
    312     Derived& derived() { return *static_cast<Derived*>(this); }
    313     const Derived& derived() const { return *static_cast<const Derived*>(this); }
    314 
    315     inline Index rows() const { return m_matrix.rows(); }
    316     inline Index cols() const { return m_matrix.cols(); }
    317 
    318     /** \returns a reference to the Super LU option object to configure the  Super LU algorithms. */
    319     inline superlu_options_t& options() { return m_sluOptions; }
    320 
    321     /** \brief Reports whether previous computation was successful.
    322       *
    323       * \returns \c Success if computation was succesful,
    324       *          \c NumericalIssue if the matrix.appears to be negative.
    325       */
    326     ComputationInfo info() const
    327     {
    328       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
    329       return m_info;
    330     }
    331 
    332     /** Computes the sparse Cholesky decomposition of \a matrix */
    333     void compute(const MatrixType& matrix)
    334     {
    335       derived().analyzePattern(matrix);
    336       derived().factorize(matrix);
    337     }
    338 
    339     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
    340       *
    341       * \sa compute()
    342       */
    343     template<typename Rhs>
    344     inline const internal::solve_retval<SuperLUBase, Rhs> solve(const MatrixBase<Rhs>& b) const
    345     {
    346       eigen_assert(m_isInitialized && "SuperLU is not initialized.");
    347       eigen_assert(rows()==b.rows()
    348                 && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
    349       return internal::solve_retval<SuperLUBase, Rhs>(*this, b.derived());
    350     }
    351 
    352     /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
    353       *
    354       * \sa compute()
    355       */
    356 //     template<typename Rhs>
    357 //     inline const internal::sparse_solve_retval<SuperLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
    358 //     {
    359 //       eigen_assert(m_isInitialized && "SuperLU is not initialized.");
    360 //       eigen_assert(rows()==b.rows()
    361 //                 && "SuperLU::solve(): invalid number of rows of the right hand side matrix b");
    362 //       return internal::sparse_solve_retval<SuperLU, Rhs>(*this, b.derived());
    363 //     }
    364 
    365     /** Performs a symbolic decomposition on the sparcity of \a matrix.
    366       *
    367       * This function is particularly useful when solving for several problems having the same structure.
    368       *
    369       * \sa factorize()
    370       */
    371     void analyzePattern(const MatrixType& /*matrix*/)
    372     {
    373       m_isInitialized = true;
    374       m_info = Success;
    375       m_analysisIsOk = true;
    376       m_factorizationIsOk = false;
    377     }
    378 
    379     template<typename Stream>
    380     void dumpMemory(Stream& s)
    381     {}
    382 
    383   protected:
    384 
    385     void initFactorization(const MatrixType& a)
    386     {
    387       set_default_options(&this->m_sluOptions);
    388 
    389       const int size = a.rows();
    390       m_matrix = a;
    391 
    392       m_sluA = internal::asSluMatrix(m_matrix);
    393       clearFactors();
    394 
    395       m_p.resize(size);
    396       m_q.resize(size);
    397       m_sluRscale.resize(size);
    398       m_sluCscale.resize(size);
    399       m_sluEtree.resize(size);
    400 
    401       // set empty B and X
    402       m_sluB.setStorageType(SLU_DN);
    403       m_sluB.setScalarType<Scalar>();
    404       m_sluB.Mtype          = SLU_GE;
    405       m_sluB.storage.values = 0;
    406       m_sluB.nrow           = 0;
    407       m_sluB.ncol           = 0;
    408       m_sluB.storage.lda    = size;
    409       m_sluX                = m_sluB;
    410 
    411       m_extractedDataAreDirty = true;
    412     }
    413 
    414     void init()
    415     {
    416       m_info = InvalidInput;
    417       m_isInitialized = false;
    418       m_sluL.Store = 0;
    419       m_sluU.Store = 0;
    420     }
    421 
    422     void extractData() const;
    423 
    424     void clearFactors()
    425     {
    426       if(m_sluL.Store)
    427         Destroy_SuperNode_Matrix(&m_sluL);
    428       if(m_sluU.Store)
    429         Destroy_CompCol_Matrix(&m_sluU);
    430 
    431       m_sluL.Store = 0;
    432       m_sluU.Store = 0;
    433 
    434       memset(&m_sluL,0,sizeof m_sluL);
    435       memset(&m_sluU,0,sizeof m_sluU);
    436     }
    437 
    438     // cached data to reduce reallocation, etc.
    439     mutable LUMatrixType m_l;
    440     mutable LUMatrixType m_u;
    441     mutable IntColVectorType m_p;
    442     mutable IntRowVectorType m_q;
    443 
    444     mutable LUMatrixType m_matrix;  // copy of the factorized matrix
    445     mutable SluMatrix m_sluA;
    446     mutable SuperMatrix m_sluL, m_sluU;
    447     mutable SluMatrix m_sluB, m_sluX;
    448     mutable SuperLUStat_t m_sluStat;
    449     mutable superlu_options_t m_sluOptions;
    450     mutable std::vector<int> m_sluEtree;
    451     mutable Matrix<RealScalar,Dynamic,1> m_sluRscale, m_sluCscale;
    452     mutable Matrix<RealScalar,Dynamic,1> m_sluFerr, m_sluBerr;
    453     mutable char m_sluEqued;
    454 
    455     mutable ComputationInfo m_info;
    456     bool m_isInitialized;
    457     int m_factorizationIsOk;
    458     int m_analysisIsOk;
    459     mutable bool m_extractedDataAreDirty;
    460 
    461   private:
    462     SuperLUBase(SuperLUBase& ) { }
    463 };
    464 
    465 
    466 /** \ingroup SuperLUSupport_Module
    467   * \class SuperLU
    468   * \brief A sparse direct LU factorization and solver based on the SuperLU library
    469   *
    470   * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
    471   * using the SuperLU library. The sparse matrix A must be squared and invertible. The vectors or matrices
    472   * X and B can be either dense or sparse.
    473   *
    474   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
    475   *
    476   * \sa \ref TutorialSparseDirectSolvers
    477   */
    478 template<typename _MatrixType>
    479 class SuperLU : public SuperLUBase<_MatrixType,SuperLU<_MatrixType> >
    480 {
    481   public:
    482     typedef SuperLUBase<_MatrixType,SuperLU> Base;
    483     typedef _MatrixType MatrixType;
    484     typedef typename Base::Scalar Scalar;
    485     typedef typename Base::RealScalar RealScalar;
    486     typedef typename Base::Index Index;
    487     typedef typename Base::IntRowVectorType IntRowVectorType;
    488     typedef typename Base::IntColVectorType IntColVectorType;
    489     typedef typename Base::LUMatrixType LUMatrixType;
    490     typedef TriangularView<LUMatrixType, Lower|UnitDiag>  LMatrixType;
    491     typedef TriangularView<LUMatrixType,  Upper>           UMatrixType;
    492 
    493   public:
    494 
    495     SuperLU() : Base() { init(); }
    496 
    497     SuperLU(const MatrixType& matrix) : Base()
    498     {
    499       Base::init();
    500       compute(matrix);
    501     }
    502 
    503     ~SuperLU()
    504     {
    505     }
    506 
    507     /** Performs a symbolic decomposition on the sparcity of \a matrix.
    508       *
    509       * This function is particularly useful when solving for several problems having the same structure.
    510       *
    511       * \sa factorize()
    512       */
    513     void analyzePattern(const MatrixType& matrix)
    514     {
    515       m_info = InvalidInput;
    516       m_isInitialized = false;
    517       Base::analyzePattern(matrix);
    518     }
    519 
    520     /** Performs a numeric decomposition of \a matrix
    521       *
    522       * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
    523       *
    524       * \sa analyzePattern()
    525       */
    526     void factorize(const MatrixType& matrix);
    527 
    528     #ifndef EIGEN_PARSED_BY_DOXYGEN
    529     /** \internal */
    530     template<typename Rhs,typename Dest>
    531     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
    532     #endif // EIGEN_PARSED_BY_DOXYGEN
    533 
    534     inline const LMatrixType& matrixL() const
    535     {
    536       if (m_extractedDataAreDirty) this->extractData();
    537       return m_l;
    538     }
    539 
    540     inline const UMatrixType& matrixU() const
    541     {
    542       if (m_extractedDataAreDirty) this->extractData();
    543       return m_u;
    544     }
    545 
    546     inline const IntColVectorType& permutationP() const
    547     {
    548       if (m_extractedDataAreDirty) this->extractData();
    549       return m_p;
    550     }
    551 
    552     inline const IntRowVectorType& permutationQ() const
    553     {
    554       if (m_extractedDataAreDirty) this->extractData();
    555       return m_q;
    556     }
    557 
    558     Scalar determinant() const;
    559 
    560   protected:
    561 
    562     using Base::m_matrix;
    563     using Base::m_sluOptions;
    564     using Base::m_sluA;
    565     using Base::m_sluB;
    566     using Base::m_sluX;
    567     using Base::m_p;
    568     using Base::m_q;
    569     using Base::m_sluEtree;
    570     using Base::m_sluEqued;
    571     using Base::m_sluRscale;
    572     using Base::m_sluCscale;
    573     using Base::m_sluL;
    574     using Base::m_sluU;
    575     using Base::m_sluStat;
    576     using Base::m_sluFerr;
    577     using Base::m_sluBerr;
    578     using Base::m_l;
    579     using Base::m_u;
    580 
    581     using Base::m_analysisIsOk;
    582     using Base::m_factorizationIsOk;
    583     using Base::m_extractedDataAreDirty;
    584     using Base::m_isInitialized;
    585     using Base::m_info;
    586 
    587     void init()
    588     {
    589       Base::init();
    590 
    591       set_default_options(&this->m_sluOptions);
    592       m_sluOptions.PrintStat        = NO;
    593       m_sluOptions.ConditionNumber  = NO;
    594       m_sluOptions.Trans            = NOTRANS;
    595       m_sluOptions.ColPerm          = COLAMD;
    596     }
    597 
    598 
    599   private:
    600     SuperLU(SuperLU& ) { }
    601 };
    602 
    603 template<typename MatrixType>
    604 void SuperLU<MatrixType>::factorize(const MatrixType& a)
    605 {
    606   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
    607   if(!m_analysisIsOk)
    608   {
    609     m_info = InvalidInput;
    610     return;
    611   }
    612 
    613   this->initFactorization(a);
    614 
    615   int info = 0;
    616   RealScalar recip_pivot_growth, rcond;
    617   RealScalar ferr, berr;
    618 
    619   StatInit(&m_sluStat);
    620   SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
    621                 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
    622                 &m_sluL, &m_sluU,
    623                 NULL, 0,
    624                 &m_sluB, &m_sluX,
    625                 &recip_pivot_growth, &rcond,
    626                 &ferr, &berr,
    627                 &m_sluStat, &info, Scalar());
    628   StatFree(&m_sluStat);
    629 
    630   m_extractedDataAreDirty = true;
    631 
    632   // FIXME how to better check for errors ???
    633   m_info = info == 0 ? Success : NumericalIssue;
    634   m_factorizationIsOk = true;
    635 }
    636 
    637 template<typename MatrixType>
    638 template<typename Rhs,typename Dest>
    639 void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
    640 {
    641   eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
    642 
    643   const int size = m_matrix.rows();
    644   const int rhsCols = b.cols();
    645   eigen_assert(size==b.rows());
    646 
    647   m_sluOptions.Trans = NOTRANS;
    648   m_sluOptions.Fact = FACTORED;
    649   m_sluOptions.IterRefine = NOREFINE;
    650 
    651 
    652   m_sluFerr.resize(rhsCols);
    653   m_sluBerr.resize(rhsCols);
    654   m_sluB = SluMatrix::Map(b.const_cast_derived());
    655   m_sluX = SluMatrix::Map(x.derived());
    656 
    657   typename Rhs::PlainObject b_cpy;
    658   if(m_sluEqued!='N')
    659   {
    660     b_cpy = b;
    661     m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
    662   }
    663 
    664   StatInit(&m_sluStat);
    665   int info = 0;
    666   RealScalar recip_pivot_growth, rcond;
    667   SuperLU_gssvx(&m_sluOptions, &m_sluA,
    668                 m_q.data(), m_p.data(),
    669                 &m_sluEtree[0], &m_sluEqued,
    670                 &m_sluRscale[0], &m_sluCscale[0],
    671                 &m_sluL, &m_sluU,
    672                 NULL, 0,
    673                 &m_sluB, &m_sluX,
    674                 &recip_pivot_growth, &rcond,
    675                 &m_sluFerr[0], &m_sluBerr[0],
    676                 &m_sluStat, &info, Scalar());
    677   StatFree(&m_sluStat);
    678   m_info = info==0 ? Success : NumericalIssue;
    679 }
    680 
    681 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
    682 //
    683 //  Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
    684 //
    685 //  THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
    686 //  EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
    687 //
    688 template<typename MatrixType, typename Derived>
    689 void SuperLUBase<MatrixType,Derived>::extractData() const
    690 {
    691   eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
    692   if (m_extractedDataAreDirty)
    693   {
    694     int         upper;
    695     int         fsupc, istart, nsupr;
    696     int         lastl = 0, lastu = 0;
    697     SCformat    *Lstore = static_cast<SCformat*>(m_sluL.Store);
    698     NCformat    *Ustore = static_cast<NCformat*>(m_sluU.Store);
    699     Scalar      *SNptr;
    700 
    701     const int size = m_matrix.rows();
    702     m_l.resize(size,size);
    703     m_l.resizeNonZeros(Lstore->nnz);
    704     m_u.resize(size,size);
    705     m_u.resizeNonZeros(Ustore->nnz);
    706 
    707     int* Lcol = m_l.outerIndexPtr();
    708     int* Lrow = m_l.innerIndexPtr();
    709     Scalar* Lval = m_l.valuePtr();
    710 
    711     int* Ucol = m_u.outerIndexPtr();
    712     int* Urow = m_u.innerIndexPtr();
    713     Scalar* Uval = m_u.valuePtr();
    714 
    715     Ucol[0] = 0;
    716     Ucol[0] = 0;
    717 
    718     /* for each supernode */
    719     for (int k = 0; k <= Lstore->nsuper; ++k)
    720     {
    721       fsupc   = L_FST_SUPC(k);
    722       istart  = L_SUB_START(fsupc);
    723       nsupr   = L_SUB_START(fsupc+1) - istart;
    724       upper   = 1;
    725 
    726       /* for each column in the supernode */
    727       for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
    728       {
    729         SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
    730 
    731         /* Extract U */
    732         for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
    733         {
    734           Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
    735           /* Matlab doesn't like explicit zero. */
    736           if (Uval[lastu] != 0.0)
    737             Urow[lastu++] = U_SUB(i);
    738         }
    739         for (int i = 0; i < upper; ++i)
    740         {
    741           /* upper triangle in the supernode */
    742           Uval[lastu] = SNptr[i];
    743           /* Matlab doesn't like explicit zero. */
    744           if (Uval[lastu] != 0.0)
    745             Urow[lastu++] = L_SUB(istart+i);
    746         }
    747         Ucol[j+1] = lastu;
    748 
    749         /* Extract L */
    750         Lval[lastl] = 1.0; /* unit diagonal */
    751         Lrow[lastl++] = L_SUB(istart + upper - 1);
    752         for (int i = upper; i < nsupr; ++i)
    753         {
    754           Lval[lastl] = SNptr[i];
    755           /* Matlab doesn't like explicit zero. */
    756           if (Lval[lastl] != 0.0)
    757             Lrow[lastl++] = L_SUB(istart+i);
    758         }
    759         Lcol[j+1] = lastl;
    760 
    761         ++upper;
    762       } /* for j ... */
    763 
    764     } /* for k ... */
    765 
    766     // squeeze the matrices :
    767     m_l.resizeNonZeros(lastl);
    768     m_u.resizeNonZeros(lastu);
    769 
    770     m_extractedDataAreDirty = false;
    771   }
    772 }
    773 
    774 template<typename MatrixType>
    775 typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
    776 {
    777   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()");
    778 
    779   if (m_extractedDataAreDirty)
    780     this->extractData();
    781 
    782   Scalar det = Scalar(1);
    783   for (int j=0; j<m_u.cols(); ++j)
    784   {
    785     if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
    786     {
    787       int lastId = m_u.outerIndexPtr()[j+1]-1;
    788       eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
    789       if (m_u.innerIndexPtr()[lastId]==j)
    790         det *= m_u.valuePtr()[lastId];
    791     }
    792   }
    793   if(m_sluEqued!='N')
    794     return det/m_sluRscale.prod()/m_sluCscale.prod();
    795   else
    796     return det;
    797 }
    798 
    799 #ifdef EIGEN_PARSED_BY_DOXYGEN
    800 #define EIGEN_SUPERLU_HAS_ILU
    801 #endif
    802 
    803 #ifdef EIGEN_SUPERLU_HAS_ILU
    804 
    805 /** \ingroup SuperLUSupport_Module
    806   * \class SuperILU
    807   * \brief A sparse direct \b incomplete LU factorization and solver based on the SuperLU library
    808   *
    809   * This class allows to solve for an approximate solution of A.X = B sparse linear problems via an incomplete LU factorization
    810   * using the SuperLU library. This class is aimed to be used as a preconditioner of the iterative linear solvers.
    811   *
    812   * \warning This class requires SuperLU 4 or later.
    813   *
    814   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
    815   *
    816   * \sa \ref TutorialSparseDirectSolvers, class ConjugateGradient, class BiCGSTAB
    817   */
    818 
    819 template<typename _MatrixType>
    820 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
    821 {
    822   public:
    823     typedef SuperLUBase<_MatrixType,SuperILU> Base;
    824     typedef _MatrixType MatrixType;
    825     typedef typename Base::Scalar Scalar;
    826     typedef typename Base::RealScalar RealScalar;
    827     typedef typename Base::Index Index;
    828 
    829   public:
    830 
    831     SuperILU() : Base() { init(); }
    832 
    833     SuperILU(const MatrixType& matrix) : Base()
    834     {
    835       init();
    836       compute(matrix);
    837     }
    838 
    839     ~SuperILU()
    840     {
    841     }
    842 
    843     /** Performs a symbolic decomposition on the sparcity of \a matrix.
    844       *
    845       * This function is particularly useful when solving for several problems having the same structure.
    846       *
    847       * \sa factorize()
    848       */
    849     void analyzePattern(const MatrixType& matrix)
    850     {
    851       Base::analyzePattern(matrix);
    852     }
    853 
    854     /** Performs a numeric decomposition of \a matrix
    855       *
    856       * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
    857       *
    858       * \sa analyzePattern()
    859       */
    860     void factorize(const MatrixType& matrix);
    861 
    862     #ifndef EIGEN_PARSED_BY_DOXYGEN
    863     /** \internal */
    864     template<typename Rhs,typename Dest>
    865     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
    866     #endif // EIGEN_PARSED_BY_DOXYGEN
    867 
    868   protected:
    869 
    870     using Base::m_matrix;
    871     using Base::m_sluOptions;
    872     using Base::m_sluA;
    873     using Base::m_sluB;
    874     using Base::m_sluX;
    875     using Base::m_p;
    876     using Base::m_q;
    877     using Base::m_sluEtree;
    878     using Base::m_sluEqued;
    879     using Base::m_sluRscale;
    880     using Base::m_sluCscale;
    881     using Base::m_sluL;
    882     using Base::m_sluU;
    883     using Base::m_sluStat;
    884     using Base::m_sluFerr;
    885     using Base::m_sluBerr;
    886     using Base::m_l;
    887     using Base::m_u;
    888 
    889     using Base::m_analysisIsOk;
    890     using Base::m_factorizationIsOk;
    891     using Base::m_extractedDataAreDirty;
    892     using Base::m_isInitialized;
    893     using Base::m_info;
    894 
    895     void init()
    896     {
    897       Base::init();
    898 
    899       ilu_set_default_options(&m_sluOptions);
    900       m_sluOptions.PrintStat        = NO;
    901       m_sluOptions.ConditionNumber  = NO;
    902       m_sluOptions.Trans            = NOTRANS;
    903       m_sluOptions.ColPerm          = MMD_AT_PLUS_A;
    904 
    905       // no attempt to preserve column sum
    906       m_sluOptions.ILU_MILU = SILU;
    907       // only basic ILU(k) support -- no direct control over memory consumption
    908       // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
    909       // and set ILU_FillFactor to max memory growth
    910       m_sluOptions.ILU_DropRule = DROP_BASIC;
    911       m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
    912     }
    913 
    914   private:
    915     SuperILU(SuperILU& ) { }
    916 };
    917 
    918 template<typename MatrixType>
    919 void SuperILU<MatrixType>::factorize(const MatrixType& a)
    920 {
    921   eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
    922   if(!m_analysisIsOk)
    923   {
    924     m_info = InvalidInput;
    925     return;
    926   }
    927 
    928   this->initFactorization(a);
    929 
    930   int info = 0;
    931   RealScalar recip_pivot_growth, rcond;
    932 
    933   StatInit(&m_sluStat);
    934   SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
    935                 &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
    936                 &m_sluL, &m_sluU,
    937                 NULL, 0,
    938                 &m_sluB, &m_sluX,
    939                 &recip_pivot_growth, &rcond,
    940                 &m_sluStat, &info, Scalar());
    941   StatFree(&m_sluStat);
    942 
    943   // FIXME how to better check for errors ???
    944   m_info = info == 0 ? Success : NumericalIssue;
    945   m_factorizationIsOk = true;
    946 }
    947 
    948 template<typename MatrixType>
    949 template<typename Rhs,typename Dest>
    950 void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
    951 {
    952   eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
    953 
    954   const int size = m_matrix.rows();
    955   const int rhsCols = b.cols();
    956   eigen_assert(size==b.rows());
    957 
    958   m_sluOptions.Trans = NOTRANS;
    959   m_sluOptions.Fact = FACTORED;
    960   m_sluOptions.IterRefine = NOREFINE;
    961 
    962   m_sluFerr.resize(rhsCols);
    963   m_sluBerr.resize(rhsCols);
    964   m_sluB = SluMatrix::Map(b.const_cast_derived());
    965   m_sluX = SluMatrix::Map(x.derived());
    966 
    967   typename Rhs::PlainObject b_cpy;
    968   if(m_sluEqued!='N')
    969   {
    970     b_cpy = b;
    971     m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
    972   }
    973 
    974   int info = 0;
    975   RealScalar recip_pivot_growth, rcond;
    976 
    977   StatInit(&m_sluStat);
    978   SuperLU_gsisx(&m_sluOptions, &m_sluA,
    979                 m_q.data(), m_p.data(),
    980                 &m_sluEtree[0], &m_sluEqued,
    981                 &m_sluRscale[0], &m_sluCscale[0],
    982                 &m_sluL, &m_sluU,
    983                 NULL, 0,
    984                 &m_sluB, &m_sluX,
    985                 &recip_pivot_growth, &rcond,
    986                 &m_sluStat, &info, Scalar());
    987   StatFree(&m_sluStat);
    988 
    989   m_info = info==0 ? Success : NumericalIssue;
    990 }
    991 #endif
    992 
    993 namespace internal {
    994 
    995 template<typename _MatrixType, typename Derived, typename Rhs>
    996 struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
    997   : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
    998 {
    999   typedef SuperLUBase<_MatrixType,Derived> Dec;
   1000   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
   1001 
   1002   template<typename Dest> void evalTo(Dest& dst) const
   1003   {
   1004     dec().derived()._solve(rhs(),dst);
   1005   }
   1006 };
   1007 
   1008 template<typename _MatrixType, typename Derived, typename Rhs>
   1009 struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
   1010   : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
   1011 {
   1012   typedef SuperLUBase<_MatrixType,Derived> Dec;
   1013   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
   1014 
   1015   template<typename Dest> void evalTo(Dest& dst) const
   1016   {
   1017     dec().derived()._solve(rhs(),dst);
   1018   }
   1019 };
   1020 
   1021 } // end namespace internal
   1022 
   1023 } // end namespace Eigen
   1024 
   1025 #endif // EIGEN_SUPERLUSUPPORT_H
   1026