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