Eigen  3.2.7
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SuperLUSupport.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2011 Gael Guennebaud <[email protected]>
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 
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 = (void*)(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 
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 
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 
319  inline superlu_options_t& options() { return m_sluOptions; }
320 
326  ComputationInfo info() const
327  {
328  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
329  return m_info;
330  }
331 
333  void compute(const MatrixType& matrix)
334  {
335  derived().analyzePattern(matrix);
336  derived().factorize(matrix);
337  }
338 
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 
356  template<typename Rhs>
357  inline const internal::sparse_solve_retval<SuperLUBase, 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<SuperLUBase, Rhs>(*this, b.derived());
363  }
364 
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 
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  init();
500  Base::compute(matrix);
501  }
502 
503  ~SuperLU()
504  {
505  }
506 
513  void analyzePattern(const MatrixType& matrix)
514  {
515  m_info = InvalidInput;
516  m_isInitialized = false;
517  Base::analyzePattern(matrix);
518  }
519 
526  void factorize(const MatrixType& matrix);
527 
528  #ifndef EIGEN_PARSED_BY_DOXYGEN
529 
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  m_sluOptions.ColPerm = COLAMD;
616  int info = 0;
617  RealScalar recip_pivot_growth, rcond;
618  RealScalar ferr, berr;
619 
620  StatInit(&m_sluStat);
621  SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
622  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
623  &m_sluL, &m_sluU,
624  NULL, 0,
625  &m_sluB, &m_sluX,
626  &recip_pivot_growth, &rcond,
627  &ferr, &berr,
628  &m_sluStat, &info, Scalar());
629  StatFree(&m_sluStat);
630 
631  m_extractedDataAreDirty = true;
632 
633  // FIXME how to better check for errors ???
634  m_info = info == 0 ? Success : NumericalIssue;
635  m_factorizationIsOk = true;
636 }
637 
638 template<typename MatrixType>
639 template<typename Rhs,typename Dest>
640 void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
641 {
642  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
643 
644  const int size = m_matrix.rows();
645  const int rhsCols = b.cols();
646  eigen_assert(size==b.rows());
647 
648  m_sluOptions.Trans = NOTRANS;
649  m_sluOptions.Fact = FACTORED;
650  m_sluOptions.IterRefine = NOREFINE;
651 
652 
653  m_sluFerr.resize(rhsCols);
654  m_sluBerr.resize(rhsCols);
655  m_sluB = SluMatrix::Map(b.const_cast_derived());
656  m_sluX = SluMatrix::Map(x.derived());
657 
658  typename Rhs::PlainObject b_cpy;
659  if(m_sluEqued!='N')
660  {
661  b_cpy = b;
662  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
663  }
664 
665  StatInit(&m_sluStat);
666  int info = 0;
667  RealScalar recip_pivot_growth, rcond;
668  SuperLU_gssvx(&m_sluOptions, &m_sluA,
669  m_q.data(), m_p.data(),
670  &m_sluEtree[0], &m_sluEqued,
671  &m_sluRscale[0], &m_sluCscale[0],
672  &m_sluL, &m_sluU,
673  NULL, 0,
674  &m_sluB, &m_sluX,
675  &recip_pivot_growth, &rcond,
676  &m_sluFerr[0], &m_sluBerr[0],
677  &m_sluStat, &info, Scalar());
678  StatFree(&m_sluStat);
679  m_info = info==0 ? Success : NumericalIssue;
680 }
681 
682 // the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
683 //
684 // Copyright (c) 1994 by Xerox Corporation. All rights reserved.
685 //
686 // THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
687 // EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
688 //
689 template<typename MatrixType, typename Derived>
690 void SuperLUBase<MatrixType,Derived>::extractData() const
691 {
692  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for extracting factors, you must first call either compute() or analyzePattern()/factorize()");
693  if (m_extractedDataAreDirty)
694  {
695  int upper;
696  int fsupc, istart, nsupr;
697  int lastl = 0, lastu = 0;
698  SCformat *Lstore = static_cast<SCformat*>(m_sluL.Store);
699  NCformat *Ustore = static_cast<NCformat*>(m_sluU.Store);
700  Scalar *SNptr;
701 
702  const int size = m_matrix.rows();
703  m_l.resize(size,size);
704  m_l.resizeNonZeros(Lstore->nnz);
705  m_u.resize(size,size);
706  m_u.resizeNonZeros(Ustore->nnz);
707 
708  int* Lcol = m_l.outerIndexPtr();
709  int* Lrow = m_l.innerIndexPtr();
710  Scalar* Lval = m_l.valuePtr();
711 
712  int* Ucol = m_u.outerIndexPtr();
713  int* Urow = m_u.innerIndexPtr();
714  Scalar* Uval = m_u.valuePtr();
715 
716  Ucol[0] = 0;
717  Ucol[0] = 0;
718 
719  /* for each supernode */
720  for (int k = 0; k <= Lstore->nsuper; ++k)
721  {
722  fsupc = L_FST_SUPC(k);
723  istart = L_SUB_START(fsupc);
724  nsupr = L_SUB_START(fsupc+1) - istart;
725  upper = 1;
726 
727  /* for each column in the supernode */
728  for (int j = fsupc; j < L_FST_SUPC(k+1); ++j)
729  {
730  SNptr = &((Scalar*)Lstore->nzval)[L_NZ_START(j)];
731 
732  /* Extract U */
733  for (int i = U_NZ_START(j); i < U_NZ_START(j+1); ++i)
734  {
735  Uval[lastu] = ((Scalar*)Ustore->nzval)[i];
736  /* Matlab doesn't like explicit zero. */
737  if (Uval[lastu] != 0.0)
738  Urow[lastu++] = U_SUB(i);
739  }
740  for (int i = 0; i < upper; ++i)
741  {
742  /* upper triangle in the supernode */
743  Uval[lastu] = SNptr[i];
744  /* Matlab doesn't like explicit zero. */
745  if (Uval[lastu] != 0.0)
746  Urow[lastu++] = L_SUB(istart+i);
747  }
748  Ucol[j+1] = lastu;
749 
750  /* Extract L */
751  Lval[lastl] = 1.0; /* unit diagonal */
752  Lrow[lastl++] = L_SUB(istart + upper - 1);
753  for (int i = upper; i < nsupr; ++i)
754  {
755  Lval[lastl] = SNptr[i];
756  /* Matlab doesn't like explicit zero. */
757  if (Lval[lastl] != 0.0)
758  Lrow[lastl++] = L_SUB(istart+i);
759  }
760  Lcol[j+1] = lastl;
761 
762  ++upper;
763  } /* for j ... */
764 
765  } /* for k ... */
766 
767  // squeeze the matrices :
768  m_l.resizeNonZeros(lastl);
769  m_u.resizeNonZeros(lastu);
770 
771  m_extractedDataAreDirty = false;
772  }
773 }
774 
775 template<typename MatrixType>
776 typename SuperLU<MatrixType>::Scalar SuperLU<MatrixType>::determinant() const
777 {
778  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()");
779 
780  if (m_extractedDataAreDirty)
781  this->extractData();
782 
783  Scalar det = Scalar(1);
784  for (int j=0; j<m_u.cols(); ++j)
785  {
786  if (m_u.outerIndexPtr()[j+1]-m_u.outerIndexPtr()[j] > 0)
787  {
788  int lastId = m_u.outerIndexPtr()[j+1]-1;
789  eigen_assert(m_u.innerIndexPtr()[lastId]<=j);
790  if (m_u.innerIndexPtr()[lastId]==j)
791  det *= m_u.valuePtr()[lastId];
792  }
793  }
794  if(m_sluEqued!='N')
795  return det/m_sluRscale.prod()/m_sluCscale.prod();
796  else
797  return det;
798 }
799 
800 #ifdef EIGEN_PARSED_BY_DOXYGEN
801 #define EIGEN_SUPERLU_HAS_ILU
802 #endif
803 
804 #ifdef EIGEN_SUPERLU_HAS_ILU
805 
820 template<typename _MatrixType>
821 class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
822 {
823  public:
824  typedef SuperLUBase<_MatrixType,SuperILU> Base;
825  typedef _MatrixType MatrixType;
826  typedef typename Base::Scalar Scalar;
827  typedef typename Base::RealScalar RealScalar;
828  typedef typename Base::Index Index;
829 
830  public:
831 
832  SuperILU() : Base() { init(); }
833 
834  SuperILU(const MatrixType& matrix) : Base()
835  {
836  init();
837  Base::compute(matrix);
838  }
839 
840  ~SuperILU()
841  {
842  }
843 
850  void analyzePattern(const MatrixType& matrix)
851  {
852  Base::analyzePattern(matrix);
853  }
854 
861  void factorize(const MatrixType& matrix);
862 
863  #ifndef EIGEN_PARSED_BY_DOXYGEN
864 
865  template<typename Rhs,typename Dest>
866  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
867  #endif // EIGEN_PARSED_BY_DOXYGEN
868 
869  protected:
870 
871  using Base::m_matrix;
872  using Base::m_sluOptions;
873  using Base::m_sluA;
874  using Base::m_sluB;
875  using Base::m_sluX;
876  using Base::m_p;
877  using Base::m_q;
878  using Base::m_sluEtree;
879  using Base::m_sluEqued;
880  using Base::m_sluRscale;
881  using Base::m_sluCscale;
882  using Base::m_sluL;
883  using Base::m_sluU;
884  using Base::m_sluStat;
885  using Base::m_sluFerr;
886  using Base::m_sluBerr;
887  using Base::m_l;
888  using Base::m_u;
889 
890  using Base::m_analysisIsOk;
891  using Base::m_factorizationIsOk;
892  using Base::m_extractedDataAreDirty;
893  using Base::m_isInitialized;
894  using Base::m_info;
895 
896  void init()
897  {
898  Base::init();
899 
900  ilu_set_default_options(&m_sluOptions);
901  m_sluOptions.PrintStat = NO;
902  m_sluOptions.ConditionNumber = NO;
903  m_sluOptions.Trans = NOTRANS;
904  m_sluOptions.ColPerm = MMD_AT_PLUS_A;
905 
906  // no attempt to preserve column sum
907  m_sluOptions.ILU_MILU = SILU;
908  // only basic ILU(k) support -- no direct control over memory consumption
909  // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
910  // and set ILU_FillFactor to max memory growth
911  m_sluOptions.ILU_DropRule = DROP_BASIC;
912  m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision()*10;
913  }
914 
915  private:
916  SuperILU(SuperILU& ) { }
917 };
918 
919 template<typename MatrixType>
920 void SuperILU<MatrixType>::factorize(const MatrixType& a)
921 {
922  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
923  if(!m_analysisIsOk)
924  {
925  m_info = InvalidInput;
926  return;
927  }
928 
929  this->initFactorization(a);
930 
931  int info = 0;
932  RealScalar recip_pivot_growth, rcond;
933 
934  StatInit(&m_sluStat);
935  SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
936  &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
937  &m_sluL, &m_sluU,
938  NULL, 0,
939  &m_sluB, &m_sluX,
940  &recip_pivot_growth, &rcond,
941  &m_sluStat, &info, Scalar());
942  StatFree(&m_sluStat);
943 
944  // FIXME how to better check for errors ???
945  m_info = info == 0 ? Success : NumericalIssue;
946  m_factorizationIsOk = true;
947 }
948 
949 template<typename MatrixType>
950 template<typename Rhs,typename Dest>
951 void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
952 {
953  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or analyzePattern()/factorize()");
954 
955  const int size = m_matrix.rows();
956  const int rhsCols = b.cols();
957  eigen_assert(size==b.rows());
958 
959  m_sluOptions.Trans = NOTRANS;
960  m_sluOptions.Fact = FACTORED;
961  m_sluOptions.IterRefine = NOREFINE;
962 
963  m_sluFerr.resize(rhsCols);
964  m_sluBerr.resize(rhsCols);
965  m_sluB = SluMatrix::Map(b.const_cast_derived());
966  m_sluX = SluMatrix::Map(x.derived());
967 
968  typename Rhs::PlainObject b_cpy;
969  if(m_sluEqued!='N')
970  {
971  b_cpy = b;
972  m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
973  }
974 
975  int info = 0;
976  RealScalar recip_pivot_growth, rcond;
977 
978  StatInit(&m_sluStat);
979  SuperLU_gsisx(&m_sluOptions, &m_sluA,
980  m_q.data(), m_p.data(),
981  &m_sluEtree[0], &m_sluEqued,
982  &m_sluRscale[0], &m_sluCscale[0],
983  &m_sluL, &m_sluU,
984  NULL, 0,
985  &m_sluB, &m_sluX,
986  &recip_pivot_growth, &rcond,
987  &m_sluStat, &info, Scalar());
988  StatFree(&m_sluStat);
989 
990  m_info = info==0 ? Success : NumericalIssue;
991 }
992 #endif
993 
994 namespace internal {
995 
996 template<typename _MatrixType, typename Derived, typename Rhs>
997 struct solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
998  : solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
999 {
1000  typedef SuperLUBase<_MatrixType,Derived> Dec;
1001  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
1002 
1003  template<typename Dest> void evalTo(Dest& dst) const
1004  {
1005  dec().derived()._solve(rhs(),dst);
1006  }
1007 };
1008 
1009 template<typename _MatrixType, typename Derived, typename Rhs>
1010 struct sparse_solve_retval<SuperLUBase<_MatrixType,Derived>, Rhs>
1011  : sparse_solve_retval_base<SuperLUBase<_MatrixType,Derived>, Rhs>
1012 {
1013  typedef SuperLUBase<_MatrixType,Derived> Dec;
1014  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
1015 
1016  template<typename Dest> void evalTo(Dest& dst) const
1017  {
1018  this->defaultEvalTo(dst);
1019  }
1020 };
1021 
1022 } // end namespace internal
1023 
1024 } // end namespace Eigen
1025 
1026 #endif // EIGEN_SUPERLUSUPPORT_H
Eigen::SparseMatrix::rows
Index rows() const
Definition: SparseMatrix.h:119
Eigen::internal::map_superlu
MappedSparseMatrix< Scalar, Flags, Index > map_superlu(SluMatrix &sluMat)
Definition: SuperLUSupport.h:272
Eigen::SparseMatrix::cols
Index cols() const
Definition: SparseMatrix.h:121
Eigen::Lower
Definition: Constants.h:167
Eigen::SuperLUBase::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuperLUSupport.h:326
Eigen::SuperILU::factorize
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:920
Eigen::SuperLUBase::analyzePattern
void analyzePattern(const MatrixType &)
Definition: SuperLUSupport.h:371
Eigen::SuperLU
A sparse direct LU factorization and solver based on the SuperLU library.
Definition: SuperLUSupport.h:479
Eigen::SuperILU
A sparse direct incomplete LU factorization and solver based on the SuperLU library.
Definition: SuperLUSupport.h:821
Eigen::NumericalIssue
Definition: Constants.h:378
Eigen::SuperLU::factorize
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:604
Eigen::SuperLUBase::solve
const internal::solve_retval< SuperLUBase, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SuperLUSupport.h:344
Eigen::NumTraits
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:88
Eigen::ColMajor
Definition: Constants.h:264
Eigen::Upper
Definition: Constants.h:169
Eigen::SuperLUBase::options
superlu_options_t & options()
Definition: SuperLUSupport.h:319
Eigen::SparseMatrixBase
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:239
Eigen::EigenBase::derived
Derived & derived()
Definition: EigenBase.h:34
Eigen::SuperLUBase
The base class for the direct and incomplete LU factorization of SuperLU.
Definition: SuperLUSupport.h:291
Eigen::SuperLUBase::compute
void compute(const MatrixType &matrix)
Definition: SuperLUSupport.h:333
Eigen::InvalidInput
Definition: Constants.h:383
Eigen::SuperILU::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:850
Eigen::SelfAdjoint
Definition: Constants.h:183
Eigen::TriangularView
Base class for triangular part in a matrix.
Definition: TriangularMatrix.h:158
Eigen::RowMajor
Definition: Constants.h:266
Eigen::Success
Definition: Constants.h:376
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:53
Eigen::SuperLU::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:513
Eigen::PlainObjectBase::resize
void resize(Index nbRows, Index nbCols)
Definition: PlainObjectBase.h:235
Eigen::SparseMatrixBase::rows
Index rows() const
Definition: SparseMatrixBase.h:159
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:374
Eigen::MappedSparseMatrix
Sparse matrix.
Definition: MappedSparseMatrix.h:31
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Eigen::SuperLUBase::solve
const internal::sparse_solve_retval< SuperLUBase, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
Definition: SuperLUSupport.h:357