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