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