source: pacpussensors/trunk/Vislab/lib3dv/eigen/Eigen/src/PaStiXSupport/PaStiXSupport.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 23.0 KB
Line 
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@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_PASTIXSUPPORT_H
11#define EIGEN_PASTIXSUPPORT_H
12
13#if defined(DCOMPLEX)
14 #define PASTIX_COMPLEX COMPLEX
15 #define PASTIX_DCOMPLEX DCOMPLEX
16#else
17 #define PASTIX_COMPLEX std::complex<float>
18 #define PASTIX_DCOMPLEX std::complex<double>
19#endif
20
21namespace Eigen {
22
23/** \ingroup PaStiXSupport_Module
24 * \brief Interface to the PaStix solver
25 *
26 * This class is used to solve the linear systems A.X = B via the PaStix library.
27 * The matrix can be either real or complex, symmetric or not.
28 *
29 * \sa TutorialSparseDirectSolvers
30 */
31template<typename _MatrixType, bool IsStrSym = false> class PastixLU;
32template<typename _MatrixType, int Options> class PastixLLT;
33template<typename _MatrixType, int Options> class PastixLDLT;
34
35namespace internal
36{
37
38 template<class Pastix> struct pastix_traits;
39
40 template<typename _MatrixType>
41 struct pastix_traits< PastixLU<_MatrixType> >
42 {
43 typedef _MatrixType MatrixType;
44 typedef typename _MatrixType::Scalar Scalar;
45 typedef typename _MatrixType::RealScalar RealScalar;
46 typedef typename _MatrixType::Index Index;
47 };
48
49 template<typename _MatrixType, int Options>
50 struct pastix_traits< PastixLLT<_MatrixType,Options> >
51 {
52 typedef _MatrixType MatrixType;
53 typedef typename _MatrixType::Scalar Scalar;
54 typedef typename _MatrixType::RealScalar RealScalar;
55 typedef typename _MatrixType::Index Index;
56 };
57
58 template<typename _MatrixType, int Options>
59 struct pastix_traits< PastixLDLT<_MatrixType,Options> >
60 {
61 typedef _MatrixType MatrixType;
62 typedef typename _MatrixType::Scalar Scalar;
63 typedef typename _MatrixType::RealScalar RealScalar;
64 typedef typename _MatrixType::Index Index;
65 };
66
67 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, float *vals, int *perm, int * invp, float *x, int nbrhs, int *iparm, double *dparm)
68 {
69 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
70 if (nbrhs == 0) {x = NULL; nbrhs=1;}
71 s_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
72 }
73
74 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, double *vals, int *perm, int * invp, double *x, int nbrhs, int *iparm, double *dparm)
75 {
76 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
77 if (nbrhs == 0) {x = NULL; nbrhs=1;}
78 d_pastix(pastix_data, pastix_comm, n, ptr, idx, vals, perm, invp, x, nbrhs, iparm, dparm);
79 }
80
81 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<float> *vals, int *perm, int * invp, std::complex<float> *x, int nbrhs, int *iparm, double *dparm)
82 {
83 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
84 if (nbrhs == 0) {x = NULL; nbrhs=1;}
85 c_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_COMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_COMPLEX*>(x), nbrhs, iparm, dparm);
86 }
87
88 void eigen_pastix(pastix_data_t **pastix_data, int pastix_comm, int n, int *ptr, int *idx, std::complex<double> *vals, int *perm, int * invp, std::complex<double> *x, int nbrhs, int *iparm, double *dparm)
89 {
90 if (n == 0) { ptr = NULL; idx = NULL; vals = NULL; }
91 if (nbrhs == 0) {x = NULL; nbrhs=1;}
92 z_pastix(pastix_data, pastix_comm, n, ptr, idx, reinterpret_cast<PASTIX_DCOMPLEX*>(vals), perm, invp, reinterpret_cast<PASTIX_DCOMPLEX*>(x), nbrhs, iparm, dparm);
93 }
94
95 // Convert the matrix to Fortran-style Numbering
96 template <typename MatrixType>
97 void c_to_fortran_numbering (MatrixType& mat)
98 {
99 if ( !(mat.outerIndexPtr()[0]) )
100 {
101 int i;
102 for(i = 0; i <= mat.rows(); ++i)
103 ++mat.outerIndexPtr()[i];
104 for(i = 0; i < mat.nonZeros(); ++i)
105 ++mat.innerIndexPtr()[i];
106 }
107 }
108
109 // Convert to C-style Numbering
110 template <typename MatrixType>
111 void fortran_to_c_numbering (MatrixType& mat)
112 {
113 // Check the Numbering
114 if ( mat.outerIndexPtr()[0] == 1 )
115 { // Convert to C-style numbering
116 int i;
117 for(i = 0; i <= mat.rows(); ++i)
118 --mat.outerIndexPtr()[i];
119 for(i = 0; i < mat.nonZeros(); ++i)
120 --mat.innerIndexPtr()[i];
121 }
122 }
123}
124
125// This is the base class to interface with PaStiX functions.
126// Users should not used this class directly.
127template <class Derived>
128class PastixBase : internal::noncopyable
129{
130 public:
131 typedef typename internal::pastix_traits<Derived>::MatrixType _MatrixType;
132 typedef _MatrixType MatrixType;
133 typedef typename MatrixType::Scalar Scalar;
134 typedef typename MatrixType::RealScalar RealScalar;
135 typedef typename MatrixType::Index Index;
136 typedef Matrix<Scalar,Dynamic,1> Vector;
137 typedef SparseMatrix<Scalar, ColMajor> ColSpMatrix;
138
139 public:
140
141 PastixBase() : m_initisOk(false), m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false), m_pastixdata(0), m_size(0)
142 {
143 init();
144 }
145
146 ~PastixBase()
147 {
148 clean();
149 }
150
151 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
152 *
153 * \sa compute()
154 */
155 template<typename Rhs>
156 inline const internal::solve_retval<PastixBase, Rhs>
157 solve(const MatrixBase<Rhs>& b) const
158 {
159 eigen_assert(m_isInitialized && "Pastix solver is not initialized.");
160 eigen_assert(rows()==b.rows()
161 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
162 return internal::solve_retval<PastixBase, Rhs>(*this, b.derived());
163 }
164
165 template<typename Rhs,typename Dest>
166 bool _solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const;
167
168 Derived& derived()
169 {
170 return *static_cast<Derived*>(this);
171 }
172 const Derived& derived() const
173 {
174 return *static_cast<const Derived*>(this);
175 }
176
177 /** Returns a reference to the integer vector IPARM of PaStiX parameters
178 * to modify the default parameters.
179 * The statistics related to the different phases of factorization and solve are saved here as well
180 * \sa analyzePattern() factorize()
181 */
182 Array<Index,IPARM_SIZE,1>& iparm()
183 {
184 return m_iparm;
185 }
186
187 /** Return a reference to a particular index parameter of the IPARM vector
188 * \sa iparm()
189 */
190
191 int& iparm(int idxparam)
192 {
193 return m_iparm(idxparam);
194 }
195
196 /** Returns a reference to the double vector DPARM of PaStiX parameters
197 * The statistics related to the different phases of factorization and solve are saved here as well
198 * \sa analyzePattern() factorize()
199 */
200 Array<RealScalar,IPARM_SIZE,1>& dparm()
201 {
202 return m_dparm;
203 }
204
205
206 /** Return a reference to a particular index parameter of the DPARM vector
207 * \sa dparm()
208 */
209 double& dparm(int idxparam)
210 {
211 return m_dparm(idxparam);
212 }
213
214 inline Index cols() const { return m_size; }
215 inline Index rows() const { return m_size; }
216
217 /** \brief Reports whether previous computation was successful.
218 *
219 * \returns \c Success if computation was succesful,
220 * \c NumericalIssue if the PaStiX reports a problem
221 * \c InvalidInput if the input matrix is invalid
222 *
223 * \sa iparm()
224 */
225 ComputationInfo info() const
226 {
227 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
228 return m_info;
229 }
230
231 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
232 *
233 * \sa compute()
234 */
235 template<typename Rhs>
236 inline const internal::sparse_solve_retval<PastixBase, Rhs>
237 solve(const SparseMatrixBase<Rhs>& b) const
238 {
239 eigen_assert(m_isInitialized && "Pastix LU, LLT or LDLT is not initialized.");
240 eigen_assert(rows()==b.rows()
241 && "PastixBase::solve(): invalid number of rows of the right hand side matrix b");
242 return internal::sparse_solve_retval<PastixBase, Rhs>(*this, b.derived());
243 }
244
245 protected:
246
247 // Initialize the Pastix data structure, check the matrix
248 void init();
249
250 // Compute the ordering and the symbolic factorization
251 void analyzePattern(ColSpMatrix& mat);
252
253 // Compute the numerical factorization
254 void factorize(ColSpMatrix& mat);
255
256 // Free all the data allocated by Pastix
257 void clean()
258 {
259 eigen_assert(m_initisOk && "The Pastix structure should be allocated first");
260 m_iparm(IPARM_START_TASK) = API_TASK_CLEAN;
261 m_iparm(IPARM_END_TASK) = API_TASK_CLEAN;
262 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
263 m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
264 }
265
266 void compute(ColSpMatrix& mat);
267
268 int m_initisOk;
269 int m_analysisIsOk;
270 int m_factorizationIsOk;
271 bool m_isInitialized;
272 mutable ComputationInfo m_info;
273 mutable pastix_data_t *m_pastixdata; // Data structure for pastix
274 mutable int m_comm; // The MPI communicator identifier
275 mutable Matrix<int,IPARM_SIZE,1> m_iparm; // integer vector for the input parameters
276 mutable Matrix<double,DPARM_SIZE,1> m_dparm; // Scalar vector for the input parameters
277 mutable Matrix<Index,Dynamic,1> m_perm; // Permutation vector
278 mutable Matrix<Index,Dynamic,1> m_invp; // Inverse permutation vector
279 mutable int m_size; // Size of the matrix
280};
281
282 /** Initialize the PaStiX data structure.
283 *A first call to this function fills iparm and dparm with the default PaStiX parameters
284 * \sa iparm() dparm()
285 */
286template <class Derived>
287void PastixBase<Derived>::init()
288{
289 m_size = 0;
290 m_iparm.setZero(IPARM_SIZE);
291 m_dparm.setZero(DPARM_SIZE);
292
293 m_iparm(IPARM_MODIFY_PARAMETER) = API_NO;
294 pastix(&m_pastixdata, MPI_COMM_WORLD,
295 0, 0, 0, 0,
296 0, 0, 0, 1, m_iparm.data(), m_dparm.data());
297
298 m_iparm[IPARM_MATRIX_VERIFICATION] = API_NO;
299 m_iparm[IPARM_VERBOSE] = 2;
300 m_iparm[IPARM_ORDERING] = API_ORDER_SCOTCH;
301 m_iparm[IPARM_INCOMPLETE] = API_NO;
302 m_iparm[IPARM_OOC_LIMIT] = 2000;
303 m_iparm[IPARM_RHS_MAKING] = API_RHS_B;
304 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
305
306 m_iparm(IPARM_START_TASK) = API_TASK_INIT;
307 m_iparm(IPARM_END_TASK) = API_TASK_INIT;
308 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, 0, 0, 0, (Scalar*)0,
309 0, 0, 0, 0, m_iparm.data(), m_dparm.data());
310
311 // Check the returned error
312 if(m_iparm(IPARM_ERROR_NUMBER)) {
313 m_info = InvalidInput;
314 m_initisOk = false;
315 }
316 else {
317 m_info = Success;
318 m_initisOk = true;
319 }
320}
321
322template <class Derived>
323void PastixBase<Derived>::compute(ColSpMatrix& mat)
324{
325 eigen_assert(mat.rows() == mat.cols() && "The input matrix should be squared");
326
327 analyzePattern(mat);
328 factorize(mat);
329
330 m_iparm(IPARM_MATRIX_VERIFICATION) = API_NO;
331 m_isInitialized = m_factorizationIsOk;
332}
333
334
335template <class Derived>
336void PastixBase<Derived>::analyzePattern(ColSpMatrix& mat)
337{
338 eigen_assert(m_initisOk && "The initialization of PaSTiX failed");
339
340 // clean previous calls
341 if(m_size>0)
342 clean();
343
344 m_size = mat.rows();
345 m_perm.resize(m_size);
346 m_invp.resize(m_size);
347
348 m_iparm(IPARM_START_TASK) = API_TASK_ORDERING;
349 m_iparm(IPARM_END_TASK) = API_TASK_ANALYSE;
350 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
351 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
352
353 // Check the returned error
354 if(m_iparm(IPARM_ERROR_NUMBER))
355 {
356 m_info = NumericalIssue;
357 m_analysisIsOk = false;
358 }
359 else
360 {
361 m_info = Success;
362 m_analysisIsOk = true;
363 }
364}
365
366template <class Derived>
367void PastixBase<Derived>::factorize(ColSpMatrix& mat)
368{
369// if(&m_cpyMat != &mat) m_cpyMat = mat;
370 eigen_assert(m_analysisIsOk && "The analysis phase should be called before the factorization phase");
371 m_iparm(IPARM_START_TASK) = API_TASK_NUMFACT;
372 m_iparm(IPARM_END_TASK) = API_TASK_NUMFACT;
373 m_size = mat.rows();
374
375 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, m_size, mat.outerIndexPtr(), mat.innerIndexPtr(),
376 mat.valuePtr(), m_perm.data(), m_invp.data(), 0, 0, m_iparm.data(), m_dparm.data());
377
378 // Check the returned error
379 if(m_iparm(IPARM_ERROR_NUMBER))
380 {
381 m_info = NumericalIssue;
382 m_factorizationIsOk = false;
383 m_isInitialized = false;
384 }
385 else
386 {
387 m_info = Success;
388 m_factorizationIsOk = true;
389 m_isInitialized = true;
390 }
391}
392
393/* Solve the system */
394template<typename Base>
395template<typename Rhs,typename Dest>
396bool PastixBase<Base>::_solve (const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
397{
398 eigen_assert(m_isInitialized && "The matrix should be factorized first");
399 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
400 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
401 int rhs = 1;
402
403 x = b; /* on return, x is overwritten by the computed solution */
404
405 for (int i = 0; i < b.cols(); i++){
406 m_iparm[IPARM_START_TASK] = API_TASK_SOLVE;
407 m_iparm[IPARM_END_TASK] = API_TASK_REFINE;
408
409 internal::eigen_pastix(&m_pastixdata, MPI_COMM_WORLD, x.rows(), 0, 0, 0,
410 m_perm.data(), m_invp.data(), &x(0, i), rhs, m_iparm.data(), m_dparm.data());
411 }
412
413 // Check the returned error
414 m_info = m_iparm(IPARM_ERROR_NUMBER)==0 ? Success : NumericalIssue;
415
416 return m_iparm(IPARM_ERROR_NUMBER)==0;
417}
418
419/** \ingroup PaStiXSupport_Module
420 * \class PastixLU
421 * \brief Sparse direct LU solver based on PaStiX library
422 *
423 * This class is used to solve the linear systems A.X = B with a supernodal LU
424 * factorization in the PaStiX library. The matrix A should be squared and nonsingular
425 * PaStiX requires that the matrix A has a symmetric structural pattern.
426 * This interface can symmetrize the input matrix otherwise.
427 * The vectors or matrices X and B can be either dense or sparse.
428 *
429 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
430 * \tparam IsStrSym Indicates if the input matrix has a symmetric pattern, default is false
431 * NOTE : Note that if the analysis and factorization phase are called separately,
432 * the input matrix will be symmetrized at each call, hence it is advised to
433 * symmetrize the matrix in a end-user program and set \p IsStrSym to true
434 *
435 * \sa \ref TutorialSparseDirectSolvers
436 *
437 */
438template<typename _MatrixType, bool IsStrSym>
439class PastixLU : public PastixBase< PastixLU<_MatrixType> >
440{
441 public:
442 typedef _MatrixType MatrixType;
443 typedef PastixBase<PastixLU<MatrixType> > Base;
444 typedef typename Base::ColSpMatrix ColSpMatrix;
445 typedef typename MatrixType::Index Index;
446
447 public:
448 PastixLU() : Base()
449 {
450 init();
451 }
452
453 PastixLU(const MatrixType& matrix):Base()
454 {
455 init();
456 compute(matrix);
457 }
458 /** Compute the LU supernodal factorization of \p matrix.
459 * iparm and dparm can be used to tune the PaStiX parameters.
460 * see the PaStiX user's manual
461 * \sa analyzePattern() factorize()
462 */
463 void compute (const MatrixType& matrix)
464 {
465 m_structureIsUptodate = false;
466 ColSpMatrix temp;
467 grabMatrix(matrix, temp);
468 Base::compute(temp);
469 }
470 /** Compute the LU symbolic factorization of \p matrix using its sparsity pattern.
471 * Several ordering methods can be used at this step. See the PaStiX user's manual.
472 * The result of this operation can be used with successive matrices having the same pattern as \p matrix
473 * \sa factorize()
474 */
475 void analyzePattern(const MatrixType& matrix)
476 {
477 m_structureIsUptodate = false;
478 ColSpMatrix temp;
479 grabMatrix(matrix, temp);
480 Base::analyzePattern(temp);
481 }
482
483 /** Compute the LU supernodal factorization of \p matrix
484 * WARNING The matrix \p matrix should have the same structural pattern
485 * as the same used in the analysis phase.
486 * \sa analyzePattern()
487 */
488 void factorize(const MatrixType& matrix)
489 {
490 ColSpMatrix temp;
491 grabMatrix(matrix, temp);
492 Base::factorize(temp);
493 }
494 protected:
495
496 void init()
497 {
498 m_structureIsUptodate = false;
499 m_iparm(IPARM_SYM) = API_SYM_NO;
500 m_iparm(IPARM_FACTORIZATION) = API_FACT_LU;
501 }
502
503 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
504 {
505 if(IsStrSym)
506 out = matrix;
507 else
508 {
509 if(!m_structureIsUptodate)
510 {
511 // update the transposed structure
512 m_transposedStructure = matrix.transpose();
513
514 // Set the elements of the matrix to zero
515 for (Index j=0; j<m_transposedStructure.outerSize(); ++j)
516 for(typename ColSpMatrix::InnerIterator it(m_transposedStructure, j); it; ++it)
517 it.valueRef() = 0.0;
518
519 m_structureIsUptodate = true;
520 }
521
522 out = m_transposedStructure + matrix;
523 }
524 internal::c_to_fortran_numbering(out);
525 }
526
527 using Base::m_iparm;
528 using Base::m_dparm;
529
530 ColSpMatrix m_transposedStructure;
531 bool m_structureIsUptodate;
532};
533
534/** \ingroup PaStiXSupport_Module
535 * \class PastixLLT
536 * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library
537 *
538 * This class is used to solve the linear systems A.X = B via a LL^T supernodal Cholesky factorization
539 * available in the PaStiX library. The matrix A should be symmetric and positive definite
540 * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX
541 * The vectors or matrices X and B can be either dense or sparse
542 *
543 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
544 * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX
545 *
546 * \sa \ref TutorialSparseDirectSolvers
547 */
548template<typename _MatrixType, int _UpLo>
549class PastixLLT : public PastixBase< PastixLLT<_MatrixType, _UpLo> >
550{
551 public:
552 typedef _MatrixType MatrixType;
553 typedef PastixBase<PastixLLT<MatrixType, _UpLo> > Base;
554 typedef typename Base::ColSpMatrix ColSpMatrix;
555
556 public:
557 enum { UpLo = _UpLo };
558 PastixLLT() : Base()
559 {
560 init();
561 }
562
563 PastixLLT(const MatrixType& matrix):Base()
564 {
565 init();
566 compute(matrix);
567 }
568
569 /** Compute the L factor of the LL^T supernodal factorization of \p matrix
570 * \sa analyzePattern() factorize()
571 */
572 void compute (const MatrixType& matrix)
573 {
574 ColSpMatrix temp;
575 grabMatrix(matrix, temp);
576 Base::compute(temp);
577 }
578
579 /** Compute the LL^T symbolic factorization of \p matrix using its sparsity pattern
580 * The result of this operation can be used with successive matrices having the same pattern as \p matrix
581 * \sa factorize()
582 */
583 void analyzePattern(const MatrixType& matrix)
584 {
585 ColSpMatrix temp;
586 grabMatrix(matrix, temp);
587 Base::analyzePattern(temp);
588 }
589 /** Compute the LL^T supernodal numerical factorization of \p matrix
590 * \sa analyzePattern()
591 */
592 void factorize(const MatrixType& matrix)
593 {
594 ColSpMatrix temp;
595 grabMatrix(matrix, temp);
596 Base::factorize(temp);
597 }
598 protected:
599 using Base::m_iparm;
600
601 void init()
602 {
603 m_iparm(IPARM_SYM) = API_SYM_YES;
604 m_iparm(IPARM_FACTORIZATION) = API_FACT_LLT;
605 }
606
607 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
608 {
609 // Pastix supports only lower, column-major matrices
610 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
611 internal::c_to_fortran_numbering(out);
612 }
613};
614
615/** \ingroup PaStiXSupport_Module
616 * \class PastixLDLT
617 * \brief A sparse direct supernodal Cholesky (LLT) factorization and solver based on the PaStiX library
618 *
619 * This class is used to solve the linear systems A.X = B via a LDL^T supernodal Cholesky factorization
620 * available in the PaStiX library. The matrix A should be symmetric and positive definite
621 * WARNING Selfadjoint complex matrices are not supported in the current version of PaStiX
622 * The vectors or matrices X and B can be either dense or sparse
623 *
624 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
625 * \tparam UpLo The part of the matrix to use : Lower or Upper. The default is Lower as required by PaStiX
626 *
627 * \sa \ref TutorialSparseDirectSolvers
628 */
629template<typename _MatrixType, int _UpLo>
630class PastixLDLT : public PastixBase< PastixLDLT<_MatrixType, _UpLo> >
631{
632 public:
633 typedef _MatrixType MatrixType;
634 typedef PastixBase<PastixLDLT<MatrixType, _UpLo> > Base;
635 typedef typename Base::ColSpMatrix ColSpMatrix;
636
637 public:
638 enum { UpLo = _UpLo };
639 PastixLDLT():Base()
640 {
641 init();
642 }
643
644 PastixLDLT(const MatrixType& matrix):Base()
645 {
646 init();
647 compute(matrix);
648 }
649
650 /** Compute the L and D factors of the LDL^T factorization of \p matrix
651 * \sa analyzePattern() factorize()
652 */
653 void compute (const MatrixType& matrix)
654 {
655 ColSpMatrix temp;
656 grabMatrix(matrix, temp);
657 Base::compute(temp);
658 }
659
660 /** Compute the LDL^T symbolic factorization of \p matrix using its sparsity pattern
661 * The result of this operation can be used with successive matrices having the same pattern as \p matrix
662 * \sa factorize()
663 */
664 void analyzePattern(const MatrixType& matrix)
665 {
666 ColSpMatrix temp;
667 grabMatrix(matrix, temp);
668 Base::analyzePattern(temp);
669 }
670 /** Compute the LDL^T supernodal numerical factorization of \p matrix
671 *
672 */
673 void factorize(const MatrixType& matrix)
674 {
675 ColSpMatrix temp;
676 grabMatrix(matrix, temp);
677 Base::factorize(temp);
678 }
679
680 protected:
681 using Base::m_iparm;
682
683 void init()
684 {
685 m_iparm(IPARM_SYM) = API_SYM_YES;
686 m_iparm(IPARM_FACTORIZATION) = API_FACT_LDLT;
687 }
688
689 void grabMatrix(const MatrixType& matrix, ColSpMatrix& out)
690 {
691 // Pastix supports only lower, column-major matrices
692 out.template selfadjointView<Lower>() = matrix.template selfadjointView<UpLo>();
693 internal::c_to_fortran_numbering(out);
694 }
695};
696
697namespace internal {
698
699template<typename _MatrixType, typename Rhs>
700struct solve_retval<PastixBase<_MatrixType>, Rhs>
701 : solve_retval_base<PastixBase<_MatrixType>, Rhs>
702{
703 typedef PastixBase<_MatrixType> Dec;
704 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
705
706 template<typename Dest> void evalTo(Dest& dst) const
707 {
708 dec()._solve(rhs(),dst);
709 }
710};
711
712template<typename _MatrixType, typename Rhs>
713struct sparse_solve_retval<PastixBase<_MatrixType>, Rhs>
714 : sparse_solve_retval_base<PastixBase<_MatrixType>, Rhs>
715{
716 typedef PastixBase<_MatrixType> Dec;
717 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
718
719 template<typename Dest> void evalTo(Dest& dst) const
720 {
721 this->defaultEvalTo(dst);
722 }
723};
724
725} // end namespace internal
726
727} // end namespace Eigen
728
729#endif
Note: See TracBrowser for help on using the repository browser.