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
13namespace 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
39DECL_GSSVX(s,float,float)
40DECL_GSSVX(c,float,std::complex<float>)
41DECL_GSSVX(d,double,double)
42DECL_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
73DECL_GSISX(s,float,float)
74DECL_GSISX(c,float,std::complex<float>)
75DECL_GSISX(d,double,double)
76DECL_GSISX(z,double,std::complex<double>)
77
78#endif
79
80template<typename MatrixType>
81struct SluMatrixMapHelper;
82
90struct 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
205template<typename Scalar, int Rows, int Cols, int Options, int MRows, int MCols>
206struct 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
224template<typename Derived>
225struct 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
262namespace internal {
263
264template<typename MatrixType>
265SluMatrix asSluMatrix(MatrixType& mat)
266{
267 return SluMatrix::Map(mat);
268}
269
271template<typename Scalar, int Flags, typename Index>
272MappedSparseMatrix<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
290template<typename _MatrixType, typename Derived>
291class 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<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
478template<typename _MatrixType>
479class 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
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
603template<typename MatrixType>
604void 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
637template<typename MatrixType>
638template<typename Rhs,typename Dest>
639void SuperLU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
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//
688template<typename MatrixType, typename Derived>
689void 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
774template<typename MatrixType>
775typename 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
818
819template<typename _MatrixType>
820class SuperILU : public SuperLUBase<_MatrixType,SuperILU<_MatrixType> >
821{
822 public:
823 typedef SuperLUBase<_MatrixType,SuperILU> Base;
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
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
918template<typename MatrixType>
919void 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
948template<typename MatrixType>
949template<typename Rhs,typename Dest>
950void SuperILU<MatrixType>::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest>& x) const
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
993namespace internal {
994
995template<typename _MatrixType, typename Derived, typename Rhs>
996struct 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
1008template<typename _MatrixType, typename Derived, typename Rhs>
1009struct 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
Eigen::MappedSparseMatrix
Sparse matrix.
Definition MappedSparseMatrix.h:33
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:50
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:129
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::resize
void resize(Index rows, Index cols)
Definition PlainObjectBase.h:219
Eigen::SparseMatrixBase
Base class of any sparse matrices or sparse expressions.
Definition SparseMatrixBase.h:27
Eigen::SparseMatrix
A versatible sparse matrix representation.
Definition SparseMatrix.h:87
Eigen::SuperILU
A sparse direct incomplete LU factorization and solver based on the SuperLU library.
Definition SuperLUSupport.h:821
Eigen::SuperILU::factorize
void factorize(const MatrixType &matrix)
Definition SuperLUSupport.h:919
Eigen::SuperILU::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition SuperLUSupport.h:849
Eigen::SuperLUBase::compute
void compute(const MatrixType &matrix)
Definition SuperLUSupport.h:333
Eigen::SuperLUBase::solve
const internal::solve_retval< SuperLUBase, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition SuperLUSupport.h:344
Eigen::SuperLUBase::options
superlu_options_t & options()
Definition SuperLUSupport.h:319
Eigen::SuperLUBase::analyzePattern
void analyzePattern(const MatrixType &)
Definition SuperLUSupport.h:371
Eigen::SuperLUBase::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition SuperLUSupport.h:326
Eigen::SuperLU
A sparse direct LU factorization and solver based on the SuperLU library.
Definition SuperLUSupport.h:480
Eigen::SuperLU::factorize
void factorize(const MatrixType &matrix)
Definition SuperLUSupport.h:604
Eigen::SuperLU::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition SuperLUSupport.h:513
Eigen::TriangularView
Base class for triangular part in a matrix.
Definition TriangularMatrix.h:160
Eigen::ComputationInfo
ComputationInfo
Definition Constants.h:367
Eigen::RowMajor
@ RowMajor
Definition Constants.h:259
Eigen::ColMajor
@ ColMajor
Definition Constants.h:257
Eigen::NumericalIssue
@ NumericalIssue
Definition Constants.h:371
Eigen::InvalidInput
@ InvalidInput
Definition Constants.h:376
Eigen::Success
@ Success
Definition Constants.h:369
Eigen::SelfAdjoint
@ SelfAdjoint
Definition Constants.h:178
Eigen::Upper
@ Upper
Definition Constants.h:164
Eigen::Lower
@ Lower
Definition Constants.h:162
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition Constants.h:48
Eigen::internal
Definition LDLT.h:18
Eigen::internal::map_superlu
MappedSparseMatrix< Scalar, Flags, Index > map_superlu(SluMatrix &sluMat)
Definition SuperLUSupport.h:272