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