Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
15namespace 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
66DECL_GSSVX(s,float,float)
67DECL_GSSVX(c,float,std::complex<float>)
68DECL_GSSVX(d,double,double)
69DECL_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
100DECL_GSISX(s,float,float)
101DECL_GSISX(c,float,std::complex<float>)
102DECL_GSISX(d,double,double)
103DECL_GSISX(z,double,std::complex<double>)
104
105#endif
106
107template<typename MatrixType>
108struct SluMatrixMapHelper;
109
117struct 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
233template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
234struct 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
252template<typename Derived>
253struct 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
290namespace internal {
291
292template<typename MatrixType>
293SluMatrix asSluMatrix(MatrixType& mat)
294{
295 return SluMatrix::Map(mat);
296}
297
299template<typename Scalar, int Flags, typename Index>
300Map<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
318template<typename MatrixType_, typename Derived>
319class SuperLUBase : public SparseSolverBase<Derived>
320{
321 protected:
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;
335 enum {
336 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
337 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
338 };
339
340 public:
341
342 SuperLUBase() {}
343
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
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
489template<typename MatrixType_>
490class SuperLU : public SuperLUBase<MatrixType_,SuperLU<MatrixType_> >
491{
492 public:
494 typedef MatrixType_ MatrixType;
495 typedef typename Base::Scalar Scalar;
496 typedef typename Base::RealScalar RealScalar;
497 typedef typename Base::StorageIndex StorageIndex;
500 typedef typename Base::PermutationMap PermutationMap;
501 typedef typename Base::LUMatrixType LUMatrixType;
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
614template<typename MatrixType>
615void 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
649template<typename MatrixType>
650template<typename Rhs,typename Dest>
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
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//
707template<typename MatrixType, typename Derived>
708void 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
793template<typename MatrixType>
794typename 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
842template<typename MatrixType_>
843class SuperILU : public SuperLUBase<MatrixType_,SuperILU<MatrixType_> >
844{
845 public:
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
941template<typename MatrixType>
942void 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
972template<typename MatrixType>
973template<typename Rhs,typename Dest>
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
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_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:65
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:98
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:283
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:285
Index cols() const
Definition: SparseMatrix.h:142
Index rows() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
A sparse direct incomplete LU factorization and solver based on the SuperLU library.
Definition: SuperLUSupport.h:844
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:872
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:942
The base class for the direct and incomplete LU factorization of SuperLU.
Definition: SuperLUSupport.h:320
void compute(const MatrixType &matrix)
Definition: SuperLUSupport.h:367
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SuperLUSupport.h:360
void analyzePattern(const MatrixType &)
Definition: SuperLUSupport.h:379
superlu_options_t & options()
Definition: SuperLUSupport.h:353
A sparse direct LU factorization and solver based on the SuperLU library.
Definition: SuperLUSupport.h:491
void factorize(const MatrixType &matrix)
Definition: SuperLUSupport.h:615
void analyzePattern(const MatrixType &matrix)
Definition: SuperLUSupport.h:526
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:191
ComputationInfo
Definition: Constants.h:442
@ SelfAdjoint
Definition: Constants.h:227
@ Lower
Definition: Constants.h:211
@ Upper
Definition: Constants.h:213
@ NumericalIssue
Definition: Constants.h:446
@ InvalidInput
Definition: Constants.h:451
@ Success
Definition: Constants.h:444
@ ColMajor
Definition: Constants.h:321
@ RowMajor
Definition: Constants.h:323
const unsigned int RowMajorBit
Definition: Constants.h:68
Namespace containing all symbols from the Eigen library.
Definition: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:235