Eigen  5.0.1-dev+60122df6
 
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
495 void factorize(const MatrixType &matrix);
496
498 template <typename Rhs, typename Dest>
499 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
500
501 inline const LMatrixType &matrixL() const {
502 if (m_extractedDataAreDirty) this->extractData();
503 return m_l;
504 }
505
506 inline const UMatrixType &matrixU() const {
507 if (m_extractedDataAreDirty) this->extractData();
508 return m_u;
509 }
510
511 inline const IntColVectorType &permutationP() const {
512 if (m_extractedDataAreDirty) this->extractData();
513 return m_p;
514 }
515
516 inline const IntRowVectorType &permutationQ() const {
517 if (m_extractedDataAreDirty) this->extractData();
518 return m_q;
519 }
520
521 Scalar determinant() const;
522
523 protected:
524 using Base::m_l;
525 using Base::m_matrix;
526 using Base::m_p;
527 using Base::m_q;
528 using Base::m_sluA;
529 using Base::m_sluB;
530 using Base::m_sluBerr;
531 using Base::m_sluCscale;
532 using Base::m_sluEqued;
533 using Base::m_sluEtree;
534 using Base::m_sluFerr;
535 using Base::m_sluL;
536 using Base::m_sluOptions;
537 using Base::m_sluRscale;
538 using Base::m_sluStat;
539 using Base::m_sluU;
540 using Base::m_sluX;
541 using Base::m_u;
542
543 using Base::m_analysisIsOk;
544 using Base::m_extractedDataAreDirty;
545 using Base::m_factorizationIsOk;
546 using Base::m_info;
547 using Base::m_isInitialized;
548
549 void init() {
550 Base::init();
551
552 set_default_options(&this->m_sluOptions);
553 m_sluOptions.PrintStat = NO;
554 m_sluOptions.ConditionNumber = NO;
555 m_sluOptions.Trans = NOTRANS;
556 m_sluOptions.ColPerm = COLAMD;
557 }
558
559 private:
560 SuperLU(SuperLU &) {}
561};
562
563template <typename MatrixType>
564void SuperLU<MatrixType>::factorize(const MatrixType &a) {
565 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
566 if (!m_analysisIsOk) {
567 m_info = InvalidInput;
568 return;
569 }
570
571 this->initFactorization(a);
572
573 m_sluOptions.ColPerm = COLAMD;
574 int info = 0;
575 RealScalar recip_pivot_growth, rcond;
576 RealScalar ferr, berr;
577
578 StatInit(&m_sluStat);
579 SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], &m_sluEqued, &m_sluRscale[0],
580 &m_sluCscale[0], &m_sluL, &m_sluU, NULL, 0, &m_sluB, &m_sluX, &recip_pivot_growth, &rcond, &ferr, &berr,
581 &m_sluStat, &info, Scalar());
582 StatFree(&m_sluStat);
583
584 m_extractedDataAreDirty = true;
585
586 // FIXME how to better check for errors ???
587 m_info = info == 0 ? Success : NumericalIssue;
588 m_factorizationIsOk = true;
589}
590
591template <typename MatrixType>
592template <typename Rhs, typename Dest>
593void SuperLU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const {
594 eigen_assert(m_factorizationIsOk &&
595 "The decomposition is not in a valid state for solving, you must first call either compute() or "
596 "analyzePattern()/factorize()");
597
598 const Index rhsCols = b.cols();
599 eigen_assert(m_matrix.rows() == b.rows());
600
601 m_sluOptions.Trans = NOTRANS;
602 m_sluOptions.Fact = FACTORED;
603 m_sluOptions.IterRefine = NOREFINE;
604
605 m_sluFerr.resize(rhsCols);
606 m_sluBerr.resize(rhsCols);
607
608 Ref<const Matrix<typename Rhs::Scalar, Dynamic, Dynamic, ColMajor> > b_ref(b);
609 Ref<const Matrix<typename Dest::Scalar, Dynamic, Dynamic, ColMajor> > x_ref(x);
610
611 m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
612 m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
613
614 typename Rhs::PlainObject b_cpy;
615 if (m_sluEqued != 'N') {
616 b_cpy = b;
617 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
618 }
619
620 StatInit(&m_sluStat);
621 int info = 0;
622 RealScalar recip_pivot_growth, rcond;
623 SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], &m_sluEqued, &m_sluRscale[0],
624 &m_sluCscale[0], &m_sluL, &m_sluU, NULL, 0, &m_sluB, &m_sluX, &recip_pivot_growth, &rcond,
625 &m_sluFerr[0], &m_sluBerr[0], &m_sluStat, &info, Scalar());
626 StatFree(&m_sluStat);
627
628 if (x.derived().data() != x_ref.data()) x = x_ref;
629
630 m_info = info == 0 ? Success : NumericalIssue;
631}
632
633// the code of this extractData() function has been adapted from the SuperLU's Matlab support code,
634//
635// Copyright (c) 1994 by Xerox Corporation. All rights reserved.
636//
637// THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
638// EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
639//
640template <typename MatrixType, typename Derived>
641void SuperLUBase<MatrixType, Derived>::extractData() const {
642 eigen_assert(m_factorizationIsOk &&
643 "The decomposition is not in a valid state for extracting factors, you must first call either compute() "
644 "or analyzePattern()/factorize()");
645 if (m_extractedDataAreDirty) {
646 int upper;
647 int fsupc, istart, nsupr;
648 int lastl = 0, lastu = 0;
649 SCformat *Lstore = static_cast<SCformat *>(m_sluL.Store);
650 NCformat *Ustore = static_cast<NCformat *>(m_sluU.Store);
651 Scalar *SNptr;
652
653 const Index size = m_matrix.rows();
654 m_l.resize(size, size);
655 m_l.resizeNonZeros(Lstore->nnz);
656 m_u.resize(size, size);
657 m_u.resizeNonZeros(Ustore->nnz);
658
659 int *Lcol = m_l.outerIndexPtr();
660 int *Lrow = m_l.innerIndexPtr();
661 Scalar *Lval = m_l.valuePtr();
662
663 int *Ucol = m_u.outerIndexPtr();
664 int *Urow = m_u.innerIndexPtr();
665 Scalar *Uval = m_u.valuePtr();
666
667 Ucol[0] = 0;
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
783 ~SuperILU() {}
784
791 void analyzePattern(const MatrixType &matrix) { Base::analyzePattern(matrix); }
792
800 void factorize(const MatrixType &matrix);
801
802#ifndef EIGEN_PARSED_BY_DOXYGEN
804 template <typename Rhs, typename Dest>
805 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const;
806#endif // EIGEN_PARSED_BY_DOXYGEN
807
808 protected:
809 using Base::m_l;
810 using Base::m_matrix;
811 using Base::m_p;
812 using Base::m_q;
813 using Base::m_sluA;
814 using Base::m_sluB;
815 using Base::m_sluBerr;
816 using Base::m_sluCscale;
817 using Base::m_sluEqued;
818 using Base::m_sluEtree;
819 using Base::m_sluFerr;
820 using Base::m_sluL;
821 using Base::m_sluOptions;
822 using Base::m_sluRscale;
823 using Base::m_sluStat;
824 using Base::m_sluU;
825 using Base::m_sluX;
826 using Base::m_u;
827
828 using Base::m_analysisIsOk;
829 using Base::m_extractedDataAreDirty;
830 using Base::m_factorizationIsOk;
831 using Base::m_info;
832 using Base::m_isInitialized;
833
834 void init() {
835 Base::init();
836
837 ilu_set_default_options(&m_sluOptions);
838 m_sluOptions.PrintStat = NO;
839 m_sluOptions.ConditionNumber = NO;
840 m_sluOptions.Trans = NOTRANS;
841 m_sluOptions.ColPerm = MMD_AT_PLUS_A;
842
843 // no attempt to preserve column sum
844 m_sluOptions.ILU_MILU = SILU;
845 // only basic ILU(k) support -- no direct control over memory consumption
846 // better to use ILU_DropRule = DROP_BASIC | DROP_AREA
847 // and set ILU_FillFactor to max memory growth
848 m_sluOptions.ILU_DropRule = DROP_BASIC;
849 m_sluOptions.ILU_DropTol = NumTraits<Scalar>::dummy_precision() * 10;
850 }
851
852 private:
853 SuperILU(SuperILU &) {}
854};
855
856template <typename MatrixType>
857void SuperILU<MatrixType>::factorize(const MatrixType &a) {
858 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
859 if (!m_analysisIsOk) {
860 m_info = InvalidInput;
861 return;
862 }
863
864 this->initFactorization(a);
865
866 int info = 0;
867 RealScalar recip_pivot_growth, rcond;
868
869 StatInit(&m_sluStat);
870 SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], &m_sluEqued, &m_sluRscale[0],
871 &m_sluCscale[0], &m_sluL, &m_sluU, NULL, 0, &m_sluB, &m_sluX, &recip_pivot_growth, &rcond, &m_sluStat,
872 &info, Scalar());
873 StatFree(&m_sluStat);
874
875 // FIXME how to better check for errors ???
876 m_info = info == 0 ? Success : NumericalIssue;
877 m_factorizationIsOk = true;
878}
879
880#ifndef EIGEN_PARSED_BY_DOXYGEN
881template <typename MatrixType>
882template <typename Rhs, typename Dest>
883void SuperILU<MatrixType>::_solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const {
884 eigen_assert(m_factorizationIsOk &&
885 "The decomposition is not in a valid state for solving, you must first call either compute() or "
886 "analyzePattern()/factorize()");
887
888 const int rhsCols = b.cols();
889 eigen_assert(m_matrix.rows() == b.rows());
890
891 m_sluOptions.Trans = NOTRANS;
892 m_sluOptions.Fact = FACTORED;
893 m_sluOptions.IterRefine = NOREFINE;
894
895 m_sluFerr.resize(rhsCols);
896 m_sluBerr.resize(rhsCols);
897
898 Ref<const Matrix<typename Rhs::Scalar, Dynamic, Dynamic, ColMajor> > b_ref(b);
899 Ref<const Matrix<typename Dest::Scalar, Dynamic, Dynamic, ColMajor> > x_ref(x);
900
901 m_sluB = SluMatrix::Map(b_ref.const_cast_derived());
902 m_sluX = SluMatrix::Map(x_ref.const_cast_derived());
903
904 typename Rhs::PlainObject b_cpy;
905 if (m_sluEqued != 'N') {
906 b_cpy = b;
907 m_sluB = SluMatrix::Map(b_cpy.const_cast_derived());
908 }
909
910 int info = 0;
911 RealScalar recip_pivot_growth, rcond;
912
913 StatInit(&m_sluStat);
914 SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], &m_sluEqued, &m_sluRscale[0],
915 &m_sluCscale[0], &m_sluL, &m_sluU, NULL, 0, &m_sluB, &m_sluX, &recip_pivot_growth, &rcond, &m_sluStat,
916 &info, Scalar());
917 StatFree(&m_sluStat);
918
919 if (x.derived().data() != x_ref.data()) x = x_ref;
920
921 m_info = info == 0 ? Success : NumericalIssue;
922}
923#endif
924
925#endif
926
927} // end namespace Eigen
928
929#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:766
Eigen::SuperILU::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition SuperLUSupport.h:791
Eigen::SuperILU::factorize
void factorize(const MatrixType &matrix)
Definition SuperLUSupport.h:857
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:564
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