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

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

Doc
File size: 20.6 KB
Line 
1/*
2 Copyright (c) 2011, Intel Corporation. All rights reserved.
3
4 Redistribution and use in source and binary forms, with or without modification,
5 are permitted provided that the following conditions are met:
6
7 * Redistributions of source code must retain the above copyright notice, this
8 list of conditions and the following disclaimer.
9 * Redistributions in binary form must reproduce the above copyright notice,
10 this list of conditions and the following disclaimer in the documentation
11 and/or other materials provided with the distribution.
12 * Neither the name of Intel Corporation nor the names of its contributors may
13 be used to endorse or promote products derived from this software without
14 specific prior written permission.
15
16 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
17 ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
18 WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
19 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
20 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
21 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
22 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
23 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
25 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26
27 ********************************************************************************
28 * Content : Eigen bindings to Intel(R) MKL PARDISO
29 ********************************************************************************
30*/
31
32#ifndef EIGEN_PARDISOSUPPORT_H
33#define EIGEN_PARDISOSUPPORT_H
34
35namespace Eigen {
36
37template<typename _MatrixType> class PardisoLU;
38template<typename _MatrixType, int Options=Upper> class PardisoLLT;
39template<typename _MatrixType, int Options=Upper> class PardisoLDLT;
40
41namespace internal
42{
43 template<typename Index>
44 struct pardiso_run_selector
45 {
46 static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a,
47 Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
48 {
49 Index error = 0;
50 ::pardiso(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
51 return error;
52 }
53 };
54 template<>
55 struct pardiso_run_selector<long long int>
56 {
57 typedef long long int Index;
58 static Index run( _MKL_DSS_HANDLE_t pt, Index maxfct, Index mnum, Index type, Index phase, Index n, void *a,
59 Index *ia, Index *ja, Index *perm, Index nrhs, Index *iparm, Index msglvl, void *b, void *x)
60 {
61 Index error = 0;
62 ::pardiso_64(pt, &maxfct, &mnum, &type, &phase, &n, a, ia, ja, perm, &nrhs, iparm, &msglvl, b, x, &error);
63 return error;
64 }
65 };
66
67 template<class Pardiso> struct pardiso_traits;
68
69 template<typename _MatrixType>
70 struct pardiso_traits< PardisoLU<_MatrixType> >
71 {
72 typedef _MatrixType MatrixType;
73 typedef typename _MatrixType::Scalar Scalar;
74 typedef typename _MatrixType::RealScalar RealScalar;
75 typedef typename _MatrixType::Index Index;
76 };
77
78 template<typename _MatrixType, int Options>
79 struct pardiso_traits< PardisoLLT<_MatrixType, Options> >
80 {
81 typedef _MatrixType MatrixType;
82 typedef typename _MatrixType::Scalar Scalar;
83 typedef typename _MatrixType::RealScalar RealScalar;
84 typedef typename _MatrixType::Index Index;
85 };
86
87 template<typename _MatrixType, int Options>
88 struct pardiso_traits< PardisoLDLT<_MatrixType, Options> >
89 {
90 typedef _MatrixType MatrixType;
91 typedef typename _MatrixType::Scalar Scalar;
92 typedef typename _MatrixType::RealScalar RealScalar;
93 typedef typename _MatrixType::Index Index;
94 };
95
96}
97
98template<class Derived>
99class PardisoImpl
100{
101 typedef internal::pardiso_traits<Derived> Traits;
102 public:
103 typedef typename Traits::MatrixType MatrixType;
104 typedef typename Traits::Scalar Scalar;
105 typedef typename Traits::RealScalar RealScalar;
106 typedef typename Traits::Index Index;
107 typedef SparseMatrix<Scalar,RowMajor,Index> SparseMatrixType;
108 typedef Matrix<Scalar,Dynamic,1> VectorType;
109 typedef Matrix<Index, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
110 typedef Matrix<Index, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
111 typedef Array<Index,64,1,DontAlign> ParameterType;
112 enum {
113 ScalarIsComplex = NumTraits<Scalar>::IsComplex
114 };
115
116 PardisoImpl()
117 {
118 eigen_assert((sizeof(Index) >= sizeof(_INTEGER_t) && sizeof(Index) <= 8) && "Non-supported index type");
119 m_iparm.setZero();
120 m_msglvl = 0; // No output
121 m_initialized = false;
122 }
123
124 ~PardisoImpl()
125 {
126 pardisoRelease();
127 }
128
129 inline Index cols() const { return m_size; }
130 inline Index rows() const { return m_size; }
131
132 /** \brief Reports whether previous computation was successful.
133 *
134 * \returns \c Success if computation was succesful,
135 * \c NumericalIssue if the matrix appears to be negative.
136 */
137 ComputationInfo info() const
138 {
139 eigen_assert(m_initialized && "Decomposition is not initialized.");
140 return m_info;
141 }
142
143 /** \warning for advanced usage only.
144 * \returns a reference to the parameter array controlling PARDISO.
145 * See the PARDISO manual to know how to use it. */
146 ParameterType& pardisoParameterArray()
147 {
148 return m_iparm;
149 }
150
151 /** Performs a symbolic decomposition on the sparcity of \a matrix.
152 *
153 * This function is particularly useful when solving for several problems having the same structure.
154 *
155 * \sa factorize()
156 */
157 Derived& analyzePattern(const MatrixType& matrix);
158
159 /** Performs a numeric decomposition of \a matrix
160 *
161 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
162 *
163 * \sa analyzePattern()
164 */
165 Derived& factorize(const MatrixType& matrix);
166
167 Derived& compute(const MatrixType& matrix);
168
169 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
170 *
171 * \sa compute()
172 */
173 template<typename Rhs>
174 inline const internal::solve_retval<PardisoImpl, Rhs>
175 solve(const MatrixBase<Rhs>& b) const
176 {
177 eigen_assert(m_initialized && "Pardiso solver is not initialized.");
178 eigen_assert(rows()==b.rows()
179 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
180 return internal::solve_retval<PardisoImpl, Rhs>(*this, b.derived());
181 }
182
183 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
184 *
185 * \sa compute()
186 */
187 template<typename Rhs>
188 inline const internal::sparse_solve_retval<PardisoImpl, Rhs>
189 solve(const SparseMatrixBase<Rhs>& b) const
190 {
191 eigen_assert(m_initialized && "Pardiso solver is not initialized.");
192 eigen_assert(rows()==b.rows()
193 && "PardisoImpl::solve(): invalid number of rows of the right hand side matrix b");
194 return internal::sparse_solve_retval<PardisoImpl, Rhs>(*this, b.derived());
195 }
196
197 Derived& derived()
198 {
199 return *static_cast<Derived*>(this);
200 }
201 const Derived& derived() const
202 {
203 return *static_cast<const Derived*>(this);
204 }
205
206 template<typename BDerived, typename XDerived>
207 bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const;
208
209 protected:
210 void pardisoRelease()
211 {
212 if(m_initialized) // Factorization ran at least once
213 {
214 internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, -1, m_size, 0, 0, 0, m_perm.data(), 0,
215 m_iparm.data(), m_msglvl, 0, 0);
216 }
217 }
218
219 void pardisoInit(int type)
220 {
221 m_type = type;
222 bool symmetric = std::abs(m_type) < 10;
223 m_iparm[0] = 1; // No solver default
224 m_iparm[1] = 2; // use Metis for the ordering
225 m_iparm[2] = 0; // Reserved. Set to zero. (??Numbers of processors, value of OMP_NUM_THREADS??)
226 m_iparm[3] = 0; // No iterative-direct algorithm
227 m_iparm[4] = 0; // No user fill-in reducing permutation
228 m_iparm[5] = 0; // Write solution into x, b is left unchanged
229 m_iparm[6] = 0; // Not in use
230 m_iparm[7] = 2; // Max numbers of iterative refinement steps
231 m_iparm[8] = 0; // Not in use
232 m_iparm[9] = 13; // Perturb the pivot elements with 1E-13
233 m_iparm[10] = symmetric ? 0 : 1; // Use nonsymmetric permutation and scaling MPS
234 m_iparm[11] = 0; // Not in use
235 m_iparm[12] = symmetric ? 0 : 1; // Maximum weighted matching algorithm is switched-off (default for symmetric).
236 // Try m_iparm[12] = 1 in case of inappropriate accuracy
237 m_iparm[13] = 0; // Output: Number of perturbed pivots
238 m_iparm[14] = 0; // Not in use
239 m_iparm[15] = 0; // Not in use
240 m_iparm[16] = 0; // Not in use
241 m_iparm[17] = -1; // Output: Number of nonzeros in the factor LU
242 m_iparm[18] = -1; // Output: Mflops for LU factorization
243 m_iparm[19] = 0; // Output: Numbers of CG Iterations
244
245 m_iparm[20] = 0; // 1x1 pivoting
246 m_iparm[26] = 0; // No matrix checker
247 m_iparm[27] = (sizeof(RealScalar) == 4) ? 1 : 0;
248 m_iparm[34] = 1; // C indexing
249 m_iparm[36] = 0; // CSR
250 m_iparm[59] = 0; // 0 - In-Core ; 1 - Automatic switch between In-Core and Out-of-Core modes ; 2 - Out-of-Core
251
252 memset(m_pt, 0, sizeof(m_pt));
253 }
254
255 protected:
256 // cached data to reduce reallocation, etc.
257
258 void manageErrorCode(Index error)
259 {
260 switch(error)
261 {
262 case 0:
263 m_info = Success;
264 break;
265 case -4:
266 case -7:
267 m_info = NumericalIssue;
268 break;
269 default:
270 m_info = InvalidInput;
271 }
272 }
273
274 mutable SparseMatrixType m_matrix;
275 ComputationInfo m_info;
276 bool m_initialized, m_analysisIsOk, m_factorizationIsOk;
277 Index m_type, m_msglvl;
278 mutable void *m_pt[64];
279 mutable ParameterType m_iparm;
280 mutable IntColVectorType m_perm;
281 Index m_size;
282
283 private:
284 PardisoImpl(PardisoImpl &) {}
285};
286
287template<class Derived>
288Derived& PardisoImpl<Derived>::compute(const MatrixType& a)
289{
290 m_size = a.rows();
291 eigen_assert(a.rows() == a.cols());
292
293 pardisoRelease();
294 memset(m_pt, 0, sizeof(m_pt));
295 m_perm.setZero(m_size);
296 derived().getMatrix(a);
297
298 Index error;
299 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 12, m_size,
300 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
301 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
302
303 manageErrorCode(error);
304 m_analysisIsOk = true;
305 m_factorizationIsOk = true;
306 m_initialized = true;
307 return derived();
308}
309
310template<class Derived>
311Derived& PardisoImpl<Derived>::analyzePattern(const MatrixType& a)
312{
313 m_size = a.rows();
314 eigen_assert(m_size == a.cols());
315
316 pardisoRelease();
317 memset(m_pt, 0, sizeof(m_pt));
318 m_perm.setZero(m_size);
319 derived().getMatrix(a);
320
321 Index error;
322 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 11, m_size,
323 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
324 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
325
326 manageErrorCode(error);
327 m_analysisIsOk = true;
328 m_factorizationIsOk = false;
329 m_initialized = true;
330 return derived();
331}
332
333template<class Derived>
334Derived& PardisoImpl<Derived>::factorize(const MatrixType& a)
335{
336 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
337 eigen_assert(m_size == a.rows() && m_size == a.cols());
338
339 derived().getMatrix(a);
340
341 Index error;
342 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 22, m_size,
343 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
344 m_perm.data(), 0, m_iparm.data(), m_msglvl, NULL, NULL);
345
346 manageErrorCode(error);
347 m_factorizationIsOk = true;
348 return derived();
349}
350
351template<class Base>
352template<typename BDerived,typename XDerived>
353bool PardisoImpl<Base>::_solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived>& x) const
354{
355 if(m_iparm[0] == 0) // Factorization was not computed
356 return false;
357
358 //Index n = m_matrix.rows();
359 Index nrhs = Index(b.cols());
360 eigen_assert(m_size==b.rows());
361 eigen_assert(((MatrixBase<BDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major right hand sides are not supported");
362 eigen_assert(((MatrixBase<XDerived>::Flags & RowMajorBit) == 0 || nrhs == 1) && "Row-major matrices of unknowns are not supported");
363 eigen_assert(((nrhs == 1) || b.outerStride() == b.rows()));
364
365
366// switch (transposed) {
367// case SvNoTrans : m_iparm[11] = 0 ; break;
368// case SvTranspose : m_iparm[11] = 2 ; break;
369// case SvAdjoint : m_iparm[11] = 1 ; break;
370// default:
371// //std::cerr << "Eigen: transposition option \"" << transposed << "\" not supported by the PARDISO backend\n";
372// m_iparm[11] = 0;
373// }
374
375 Scalar* rhs_ptr = const_cast<Scalar*>(b.derived().data());
376 Matrix<Scalar,Dynamic,Dynamic,ColMajor> tmp;
377
378 // Pardiso cannot solve in-place
379 if(rhs_ptr == x.derived().data())
380 {
381 tmp = b;
382 rhs_ptr = tmp.data();
383 }
384
385 Index error;
386 error = internal::pardiso_run_selector<Index>::run(m_pt, 1, 1, m_type, 33, m_size,
387 m_matrix.valuePtr(), m_matrix.outerIndexPtr(), m_matrix.innerIndexPtr(),
388 m_perm.data(), nrhs, m_iparm.data(), m_msglvl,
389 rhs_ptr, x.derived().data());
390 return error==0;
391}
392
393
394/** \ingroup PardisoSupport_Module
395 * \class PardisoLU
396 * \brief A sparse direct LU factorization and solver based on the PARDISO library
397 *
398 * This class allows to solve for A.X = B sparse linear problems via a direct LU factorization
399 * using the Intel MKL PARDISO library. The sparse matrix A must be squared and invertible.
400 * The vectors or matrices X and B can be either dense or sparse.
401 *
402 * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set:
403 * \code solver.pardisoParameterArray()[59] = 1; \endcode
404 *
405 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
406 *
407 * \sa \ref TutorialSparseDirectSolvers
408 */
409template<typename MatrixType>
410class PardisoLU : public PardisoImpl< PardisoLU<MatrixType> >
411{
412 protected:
413 typedef PardisoImpl< PardisoLU<MatrixType> > Base;
414 typedef typename Base::Scalar Scalar;
415 typedef typename Base::RealScalar RealScalar;
416 using Base::pardisoInit;
417 using Base::m_matrix;
418 friend class PardisoImpl< PardisoLU<MatrixType> >;
419
420 public:
421
422 using Base::compute;
423 using Base::solve;
424
425 PardisoLU()
426 : Base()
427 {
428 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
429 }
430
431 PardisoLU(const MatrixType& matrix)
432 : Base()
433 {
434 pardisoInit(Base::ScalarIsComplex ? 13 : 11);
435 compute(matrix);
436 }
437 protected:
438 void getMatrix(const MatrixType& matrix)
439 {
440 m_matrix = matrix;
441 }
442
443 private:
444 PardisoLU(PardisoLU& ) {}
445};
446
447/** \ingroup PardisoSupport_Module
448 * \class PardisoLLT
449 * \brief A sparse direct Cholesky (LLT) factorization and solver based on the PARDISO library
450 *
451 * This class allows to solve for A.X = B sparse linear problems via a LL^T Cholesky factorization
452 * using the Intel MKL PARDISO library. The sparse matrix A must be selfajoint and positive definite.
453 * The vectors or matrices X and B can be either dense or sparse.
454 *
455 * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set:
456 * \code solver.pardisoParameterArray()[59] = 1; \endcode
457 *
458 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
459 * \tparam UpLo can be any bitwise combination of Upper, Lower. The default is Upper, meaning only the upper triangular part has to be used.
460 * Upper|Lower can be used to tell both triangular parts can be used as input.
461 *
462 * \sa \ref TutorialSparseDirectSolvers
463 */
464template<typename MatrixType, int _UpLo>
465class PardisoLLT : public PardisoImpl< PardisoLLT<MatrixType,_UpLo> >
466{
467 protected:
468 typedef PardisoImpl< PardisoLLT<MatrixType,_UpLo> > Base;
469 typedef typename Base::Scalar Scalar;
470 typedef typename Base::Index Index;
471 typedef typename Base::RealScalar RealScalar;
472 using Base::pardisoInit;
473 using Base::m_matrix;
474 friend class PardisoImpl< PardisoLLT<MatrixType,_UpLo> >;
475
476 public:
477
478 enum { UpLo = _UpLo };
479 using Base::compute;
480 using Base::solve;
481
482 PardisoLLT()
483 : Base()
484 {
485 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
486 }
487
488 PardisoLLT(const MatrixType& matrix)
489 : Base()
490 {
491 pardisoInit(Base::ScalarIsComplex ? 4 : 2);
492 compute(matrix);
493 }
494
495 protected:
496
497 void getMatrix(const MatrixType& matrix)
498 {
499 // PARDISO supports only upper, row-major matrices
500 PermutationMatrix<Dynamic,Dynamic,Index> p_null;
501 m_matrix.resize(matrix.rows(), matrix.cols());
502 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
503 }
504
505 private:
506 PardisoLLT(PardisoLLT& ) {}
507};
508
509/** \ingroup PardisoSupport_Module
510 * \class PardisoLDLT
511 * \brief A sparse direct Cholesky (LDLT) factorization and solver based on the PARDISO library
512 *
513 * This class allows to solve for A.X = B sparse linear problems via a LDL^T Cholesky factorization
514 * using the Intel MKL PARDISO library. The sparse matrix A is assumed to be selfajoint and positive definite.
515 * For complex matrices, A can also be symmetric only, see the \a Options template parameter.
516 * The vectors or matrices X and B can be either dense or sparse.
517 *
518 * By default, it runs in in-core mode. To enable PARDISO's out-of-core feature, set:
519 * \code solver.pardisoParameterArray()[59] = 1; \endcode
520 *
521 * \tparam MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
522 * \tparam Options can be any bitwise combination of Upper, Lower, and Symmetric. The default is Upper, meaning only the upper triangular part has to be used.
523 * Symmetric can be used for symmetric, non-selfadjoint complex matrices, the default being to assume a selfadjoint matrix.
524 * Upper|Lower can be used to tell both triangular parts can be used as input.
525 *
526 * \sa \ref TutorialSparseDirectSolvers
527 */
528template<typename MatrixType, int Options>
529class PardisoLDLT : public PardisoImpl< PardisoLDLT<MatrixType,Options> >
530{
531 protected:
532 typedef PardisoImpl< PardisoLDLT<MatrixType,Options> > Base;
533 typedef typename Base::Scalar Scalar;
534 typedef typename Base::Index Index;
535 typedef typename Base::RealScalar RealScalar;
536 using Base::pardisoInit;
537 using Base::m_matrix;
538 friend class PardisoImpl< PardisoLDLT<MatrixType,Options> >;
539
540 public:
541
542 using Base::compute;
543 using Base::solve;
544 enum { UpLo = Options&(Upper|Lower) };
545
546 PardisoLDLT()
547 : Base()
548 {
549 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
550 }
551
552 PardisoLDLT(const MatrixType& matrix)
553 : Base()
554 {
555 pardisoInit(Base::ScalarIsComplex ? ( bool(Options&Symmetric) ? 6 : -4 ) : -2);
556 compute(matrix);
557 }
558
559 void getMatrix(const MatrixType& matrix)
560 {
561 // PARDISO supports only upper, row-major matrices
562 PermutationMatrix<Dynamic,Dynamic,Index> p_null;
563 m_matrix.resize(matrix.rows(), matrix.cols());
564 m_matrix.template selfadjointView<Upper>() = matrix.template selfadjointView<UpLo>().twistedBy(p_null);
565 }
566
567 private:
568 PardisoLDLT(PardisoLDLT& ) {}
569};
570
571namespace internal {
572
573template<typename _Derived, typename Rhs>
574struct solve_retval<PardisoImpl<_Derived>, Rhs>
575 : solve_retval_base<PardisoImpl<_Derived>, Rhs>
576{
577 typedef PardisoImpl<_Derived> Dec;
578 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
579
580 template<typename Dest> void evalTo(Dest& dst) const
581 {
582 dec()._solve(rhs(),dst);
583 }
584};
585
586template<typename Derived, typename Rhs>
587struct sparse_solve_retval<PardisoImpl<Derived>, Rhs>
588 : sparse_solve_retval_base<PardisoImpl<Derived>, Rhs>
589{
590 typedef PardisoImpl<Derived> Dec;
591 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
592
593 template<typename Dest> void evalTo(Dest& dst) const
594 {
595 this->defaultEvalTo(dst);
596 }
597};
598
599} // end namespace internal
600
601} // end namespace Eigen
602
603#endif // EIGEN_PARDISOSUPPORT_H
Note: See TracBrowser for help on using the repository browser.