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