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