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