1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2008-2010 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_CHOLMODSUPPORT_H
|
---|
11 | #define EIGEN_CHOLMODSUPPORT_H
|
---|
12 |
|
---|
13 | namespace Eigen {
|
---|
14 |
|
---|
15 | namespace internal {
|
---|
16 |
|
---|
17 | template<typename Scalar, typename CholmodType>
|
---|
18 | void cholmod_configure_matrix(CholmodType& mat)
|
---|
19 | {
|
---|
20 | if (internal::is_same<Scalar,float>::value)
|
---|
21 | {
|
---|
22 | mat.xtype = CHOLMOD_REAL;
|
---|
23 | mat.dtype = CHOLMOD_SINGLE;
|
---|
24 | }
|
---|
25 | else if (internal::is_same<Scalar,double>::value)
|
---|
26 | {
|
---|
27 | mat.xtype = CHOLMOD_REAL;
|
---|
28 | mat.dtype = CHOLMOD_DOUBLE;
|
---|
29 | }
|
---|
30 | else if (internal::is_same<Scalar,std::complex<float> >::value)
|
---|
31 | {
|
---|
32 | mat.xtype = CHOLMOD_COMPLEX;
|
---|
33 | mat.dtype = CHOLMOD_SINGLE;
|
---|
34 | }
|
---|
35 | else if (internal::is_same<Scalar,std::complex<double> >::value)
|
---|
36 | {
|
---|
37 | mat.xtype = CHOLMOD_COMPLEX;
|
---|
38 | mat.dtype = CHOLMOD_DOUBLE;
|
---|
39 | }
|
---|
40 | else
|
---|
41 | {
|
---|
42 | eigen_assert(false && "Scalar type not supported by CHOLMOD");
|
---|
43 | }
|
---|
44 | }
|
---|
45 |
|
---|
46 | } // namespace internal
|
---|
47 |
|
---|
48 | /** Wraps the Eigen sparse matrix \a mat into a Cholmod sparse matrix object.
|
---|
49 | * Note that the data are shared.
|
---|
50 | */
|
---|
51 | template<typename _Scalar, int _Options, typename _Index>
|
---|
52 | cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
|
---|
53 | {
|
---|
54 | cholmod_sparse res;
|
---|
55 | res.nzmax = mat.nonZeros();
|
---|
56 | res.nrow = mat.rows();;
|
---|
57 | res.ncol = mat.cols();
|
---|
58 | res.p = mat.outerIndexPtr();
|
---|
59 | res.i = mat.innerIndexPtr();
|
---|
60 | res.x = mat.valuePtr();
|
---|
61 | res.z = 0;
|
---|
62 | res.sorted = 1;
|
---|
63 | if(mat.isCompressed())
|
---|
64 | {
|
---|
65 | res.packed = 1;
|
---|
66 | res.nz = 0;
|
---|
67 | }
|
---|
68 | else
|
---|
69 | {
|
---|
70 | res.packed = 0;
|
---|
71 | res.nz = mat.innerNonZeroPtr();
|
---|
72 | }
|
---|
73 |
|
---|
74 | res.dtype = 0;
|
---|
75 | res.stype = -1;
|
---|
76 |
|
---|
77 | if (internal::is_same<_Index,int>::value)
|
---|
78 | {
|
---|
79 | res.itype = CHOLMOD_INT;
|
---|
80 | }
|
---|
81 | else if (internal::is_same<_Index,SuiteSparse_long>::value)
|
---|
82 | {
|
---|
83 | res.itype = CHOLMOD_LONG;
|
---|
84 | }
|
---|
85 | else
|
---|
86 | {
|
---|
87 | eigen_assert(false && "Index type not supported yet");
|
---|
88 | }
|
---|
89 |
|
---|
90 | // setup res.xtype
|
---|
91 | internal::cholmod_configure_matrix<_Scalar>(res);
|
---|
92 |
|
---|
93 | res.stype = 0;
|
---|
94 |
|
---|
95 | return res;
|
---|
96 | }
|
---|
97 |
|
---|
98 | template<typename _Scalar, int _Options, typename _Index>
|
---|
99 | const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
|
---|
100 | {
|
---|
101 | cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
|
---|
102 | return res;
|
---|
103 | }
|
---|
104 |
|
---|
105 | /** Returns a view of the Eigen sparse matrix \a mat as Cholmod sparse matrix.
|
---|
106 | * The data are not copied but shared. */
|
---|
107 | template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
|
---|
108 | cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
|
---|
109 | {
|
---|
110 | cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
|
---|
111 |
|
---|
112 | if(UpLo==Upper) res.stype = 1;
|
---|
113 | if(UpLo==Lower) res.stype = -1;
|
---|
114 |
|
---|
115 | return res;
|
---|
116 | }
|
---|
117 |
|
---|
118 | /** Returns a view of the Eigen \b dense matrix \a mat as Cholmod dense matrix.
|
---|
119 | * The data are not copied but shared. */
|
---|
120 | template<typename Derived>
|
---|
121 | cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
|
---|
122 | {
|
---|
123 | EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
|
---|
124 | typedef typename Derived::Scalar Scalar;
|
---|
125 |
|
---|
126 | cholmod_dense res;
|
---|
127 | res.nrow = mat.rows();
|
---|
128 | res.ncol = mat.cols();
|
---|
129 | res.nzmax = res.nrow * res.ncol;
|
---|
130 | res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
|
---|
131 | res.x = (void*)(mat.derived().data());
|
---|
132 | res.z = 0;
|
---|
133 |
|
---|
134 | internal::cholmod_configure_matrix<Scalar>(res);
|
---|
135 |
|
---|
136 | return res;
|
---|
137 | }
|
---|
138 |
|
---|
139 | /** Returns a view of the Cholmod sparse matrix \a cm as an Eigen sparse matrix.
|
---|
140 | * The data are not copied but shared. */
|
---|
141 | template<typename Scalar, int Flags, typename Index>
|
---|
142 | MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
|
---|
143 | {
|
---|
144 | return MappedSparseMatrix<Scalar,Flags,Index>
|
---|
145 | (cm.nrow, cm.ncol, static_cast<Index*>(cm.p)[cm.ncol],
|
---|
146 | static_cast<Index*>(cm.p), static_cast<Index*>(cm.i),static_cast<Scalar*>(cm.x) );
|
---|
147 | }
|
---|
148 |
|
---|
149 | enum CholmodMode {
|
---|
150 | CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
|
---|
151 | };
|
---|
152 |
|
---|
153 |
|
---|
154 | /** \ingroup CholmodSupport_Module
|
---|
155 | * \class CholmodBase
|
---|
156 | * \brief The base class for the direct Cholesky factorization of Cholmod
|
---|
157 | * \sa class CholmodSupernodalLLT, class CholmodSimplicialLDLT, class CholmodSimplicialLLT
|
---|
158 | */
|
---|
159 | template<typename _MatrixType, int _UpLo, typename Derived>
|
---|
160 | class CholmodBase : internal::noncopyable
|
---|
161 | {
|
---|
162 | public:
|
---|
163 | typedef _MatrixType MatrixType;
|
---|
164 | enum { UpLo = _UpLo };
|
---|
165 | typedef typename MatrixType::Scalar Scalar;
|
---|
166 | typedef typename MatrixType::RealScalar RealScalar;
|
---|
167 | typedef MatrixType CholMatrixType;
|
---|
168 | typedef typename MatrixType::Index Index;
|
---|
169 |
|
---|
170 | public:
|
---|
171 |
|
---|
172 | CholmodBase()
|
---|
173 | : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
|
---|
174 | {
|
---|
175 | m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
|
---|
176 | cholmod_start(&m_cholmod);
|
---|
177 | }
|
---|
178 |
|
---|
179 | CholmodBase(const MatrixType& matrix)
|
---|
180 | : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
|
---|
181 | {
|
---|
182 | m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
|
---|
183 | cholmod_start(&m_cholmod);
|
---|
184 | compute(matrix);
|
---|
185 | }
|
---|
186 |
|
---|
187 | ~CholmodBase()
|
---|
188 | {
|
---|
189 | if(m_cholmodFactor)
|
---|
190 | cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
|
---|
191 | cholmod_finish(&m_cholmod);
|
---|
192 | }
|
---|
193 |
|
---|
194 | inline Index cols() const { return m_cholmodFactor->n; }
|
---|
195 | inline Index rows() const { return m_cholmodFactor->n; }
|
---|
196 |
|
---|
197 | Derived& derived() { return *static_cast<Derived*>(this); }
|
---|
198 | const Derived& derived() const { return *static_cast<const Derived*>(this); }
|
---|
199 |
|
---|
200 | /** \brief Reports whether previous computation was successful.
|
---|
201 | *
|
---|
202 | * \returns \c Success if computation was succesful,
|
---|
203 | * \c NumericalIssue if the matrix.appears to be negative.
|
---|
204 | */
|
---|
205 | ComputationInfo info() const
|
---|
206 | {
|
---|
207 | eigen_assert(m_isInitialized && "Decomposition is not initialized.");
|
---|
208 | return m_info;
|
---|
209 | }
|
---|
210 |
|
---|
211 | /** Computes the sparse Cholesky decomposition of \a matrix */
|
---|
212 | Derived& compute(const MatrixType& matrix)
|
---|
213 | {
|
---|
214 | analyzePattern(matrix);
|
---|
215 | factorize(matrix);
|
---|
216 | return derived();
|
---|
217 | }
|
---|
218 |
|
---|
219 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
|
---|
220 | *
|
---|
221 | * \sa compute()
|
---|
222 | */
|
---|
223 | template<typename Rhs>
|
---|
224 | inline const internal::solve_retval<CholmodBase, Rhs>
|
---|
225 | solve(const MatrixBase<Rhs>& b) const
|
---|
226 | {
|
---|
227 | eigen_assert(m_isInitialized && "LLT is not initialized.");
|
---|
228 | eigen_assert(rows()==b.rows()
|
---|
229 | && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
|
---|
230 | return internal::solve_retval<CholmodBase, Rhs>(*this, b.derived());
|
---|
231 | }
|
---|
232 |
|
---|
233 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
|
---|
234 | *
|
---|
235 | * \sa compute()
|
---|
236 | */
|
---|
237 | template<typename Rhs>
|
---|
238 | inline const internal::sparse_solve_retval<CholmodBase, Rhs>
|
---|
239 | solve(const SparseMatrixBase<Rhs>& b) const
|
---|
240 | {
|
---|
241 | eigen_assert(m_isInitialized && "LLT is not initialized.");
|
---|
242 | eigen_assert(rows()==b.rows()
|
---|
243 | && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
|
---|
244 | return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived());
|
---|
245 | }
|
---|
246 |
|
---|
247 | /** Performs a symbolic decomposition on the sparsity pattern of \a matrix.
|
---|
248 | *
|
---|
249 | * This function is particularly useful when solving for several problems having the same structure.
|
---|
250 | *
|
---|
251 | * \sa factorize()
|
---|
252 | */
|
---|
253 | void analyzePattern(const MatrixType& matrix)
|
---|
254 | {
|
---|
255 | if(m_cholmodFactor)
|
---|
256 | {
|
---|
257 | cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
|
---|
258 | m_cholmodFactor = 0;
|
---|
259 | }
|
---|
260 | cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
|
---|
261 | m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
|
---|
262 |
|
---|
263 | this->m_isInitialized = true;
|
---|
264 | this->m_info = Success;
|
---|
265 | m_analysisIsOk = true;
|
---|
266 | m_factorizationIsOk = false;
|
---|
267 | }
|
---|
268 |
|
---|
269 | /** Performs a numeric decomposition of \a matrix
|
---|
270 | *
|
---|
271 | * The given matrix must have the same sparsity pattern as the matrix on which the symbolic decomposition has been performed.
|
---|
272 | *
|
---|
273 | * \sa analyzePattern()
|
---|
274 | */
|
---|
275 | void factorize(const MatrixType& matrix)
|
---|
276 | {
|
---|
277 | eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
|
---|
278 | cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
|
---|
279 | cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
|
---|
280 |
|
---|
281 | // If the factorization failed, minor is the column at which it did. On success minor == n.
|
---|
282 | this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
|
---|
283 | m_factorizationIsOk = true;
|
---|
284 | }
|
---|
285 |
|
---|
286 | /** Returns a reference to the Cholmod's configuration structure to get a full control over the performed operations.
|
---|
287 | * See the Cholmod user guide for details. */
|
---|
288 | cholmod_common& cholmod() { return m_cholmod; }
|
---|
289 |
|
---|
290 | #ifndef EIGEN_PARSED_BY_DOXYGEN
|
---|
291 | /** \internal */
|
---|
292 | template<typename Rhs,typename Dest>
|
---|
293 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
|
---|
294 | {
|
---|
295 | eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
|
---|
296 | const Index size = m_cholmodFactor->n;
|
---|
297 | EIGEN_UNUSED_VARIABLE(size);
|
---|
298 | eigen_assert(size==b.rows());
|
---|
299 |
|
---|
300 | // note: cd stands for Cholmod Dense
|
---|
301 | Rhs& b_ref(b.const_cast_derived());
|
---|
302 | cholmod_dense b_cd = viewAsCholmod(b_ref);
|
---|
303 | cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
|
---|
304 | if(!x_cd)
|
---|
305 | {
|
---|
306 | this->m_info = NumericalIssue;
|
---|
307 | }
|
---|
308 | // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
|
---|
309 | dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
|
---|
310 | cholmod_free_dense(&x_cd, &m_cholmod);
|
---|
311 | }
|
---|
312 |
|
---|
313 | /** \internal */
|
---|
314 | template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
|
---|
315 | void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
|
---|
316 | {
|
---|
317 | eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
|
---|
318 | const Index size = m_cholmodFactor->n;
|
---|
319 | EIGEN_UNUSED_VARIABLE(size);
|
---|
320 | eigen_assert(size==b.rows());
|
---|
321 |
|
---|
322 | // note: cs stands for Cholmod Sparse
|
---|
323 | cholmod_sparse b_cs = viewAsCholmod(b);
|
---|
324 | cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
|
---|
325 | if(!x_cs)
|
---|
326 | {
|
---|
327 | this->m_info = NumericalIssue;
|
---|
328 | }
|
---|
329 | // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
|
---|
330 | dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
|
---|
331 | cholmod_free_sparse(&x_cs, &m_cholmod);
|
---|
332 | }
|
---|
333 | #endif // EIGEN_PARSED_BY_DOXYGEN
|
---|
334 |
|
---|
335 |
|
---|
336 | /** Sets the shift parameter that will be used to adjust the diagonal coefficients during the numerical factorization.
|
---|
337 | *
|
---|
338 | * During the numerical factorization, an offset term is added to the diagonal coefficients:\n
|
---|
339 | * \c d_ii = \a offset + \c d_ii
|
---|
340 | *
|
---|
341 | * The default is \a offset=0.
|
---|
342 | *
|
---|
343 | * \returns a reference to \c *this.
|
---|
344 | */
|
---|
345 | Derived& setShift(const RealScalar& offset)
|
---|
346 | {
|
---|
347 | m_shiftOffset[0] = offset;
|
---|
348 | return derived();
|
---|
349 | }
|
---|
350 |
|
---|
351 | template<typename Stream>
|
---|
352 | void dumpMemory(Stream& /*s*/)
|
---|
353 | {}
|
---|
354 |
|
---|
355 | protected:
|
---|
356 | mutable cholmod_common m_cholmod;
|
---|
357 | cholmod_factor* m_cholmodFactor;
|
---|
358 | RealScalar m_shiftOffset[2];
|
---|
359 | mutable ComputationInfo m_info;
|
---|
360 | bool m_isInitialized;
|
---|
361 | int m_factorizationIsOk;
|
---|
362 | int m_analysisIsOk;
|
---|
363 | };
|
---|
364 |
|
---|
365 | /** \ingroup CholmodSupport_Module
|
---|
366 | * \class CholmodSimplicialLLT
|
---|
367 | * \brief A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod
|
---|
368 | *
|
---|
369 | * This class allows to solve for A.X = B sparse linear problems via a simplicial LL^T Cholesky factorization
|
---|
370 | * using the Cholmod library.
|
---|
371 | * This simplicial variant is equivalent to Eigen's built-in SimplicialLLT class. Therefore, it has little practical interest.
|
---|
372 | * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
|
---|
373 | * X and B can be either dense or sparse.
|
---|
374 | *
|
---|
375 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
---|
376 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
---|
377 | * or Upper. Default is Lower.
|
---|
378 | *
|
---|
379 | * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
|
---|
380 | *
|
---|
381 | * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLLT
|
---|
382 | */
|
---|
383 | template<typename _MatrixType, int _UpLo = Lower>
|
---|
384 | class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
|
---|
385 | {
|
---|
386 | typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
|
---|
387 | using Base::m_cholmod;
|
---|
388 |
|
---|
389 | public:
|
---|
390 |
|
---|
391 | typedef _MatrixType MatrixType;
|
---|
392 |
|
---|
393 | CholmodSimplicialLLT() : Base() { init(); }
|
---|
394 |
|
---|
395 | CholmodSimplicialLLT(const MatrixType& matrix) : Base()
|
---|
396 | {
|
---|
397 | init();
|
---|
398 | Base::compute(matrix);
|
---|
399 | }
|
---|
400 |
|
---|
401 | ~CholmodSimplicialLLT() {}
|
---|
402 | protected:
|
---|
403 | void init()
|
---|
404 | {
|
---|
405 | m_cholmod.final_asis = 0;
|
---|
406 | m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
|
---|
407 | m_cholmod.final_ll = 1;
|
---|
408 | }
|
---|
409 | };
|
---|
410 |
|
---|
411 |
|
---|
412 | /** \ingroup CholmodSupport_Module
|
---|
413 | * \class CholmodSimplicialLDLT
|
---|
414 | * \brief A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod
|
---|
415 | *
|
---|
416 | * This class allows to solve for A.X = B sparse linear problems via a simplicial LDL^T Cholesky factorization
|
---|
417 | * using the Cholmod library.
|
---|
418 | * This simplicial variant is equivalent to Eigen's built-in SimplicialLDLT class. Therefore, it has little practical interest.
|
---|
419 | * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
|
---|
420 | * X and B can be either dense or sparse.
|
---|
421 | *
|
---|
422 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
---|
423 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
---|
424 | * or Upper. Default is Lower.
|
---|
425 | *
|
---|
426 | * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
|
---|
427 | *
|
---|
428 | * \sa \ref TutorialSparseDirectSolvers, class CholmodSupernodalLLT, class SimplicialLDLT
|
---|
429 | */
|
---|
430 | template<typename _MatrixType, int _UpLo = Lower>
|
---|
431 | class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
|
---|
432 | {
|
---|
433 | typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
|
---|
434 | using Base::m_cholmod;
|
---|
435 |
|
---|
436 | public:
|
---|
437 |
|
---|
438 | typedef _MatrixType MatrixType;
|
---|
439 |
|
---|
440 | CholmodSimplicialLDLT() : Base() { init(); }
|
---|
441 |
|
---|
442 | CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
|
---|
443 | {
|
---|
444 | init();
|
---|
445 | Base::compute(matrix);
|
---|
446 | }
|
---|
447 |
|
---|
448 | ~CholmodSimplicialLDLT() {}
|
---|
449 | protected:
|
---|
450 | void init()
|
---|
451 | {
|
---|
452 | m_cholmod.final_asis = 1;
|
---|
453 | m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
|
---|
454 | }
|
---|
455 | };
|
---|
456 |
|
---|
457 | /** \ingroup CholmodSupport_Module
|
---|
458 | * \class CholmodSupernodalLLT
|
---|
459 | * \brief A supernodal Cholesky (LLT) factorization and solver based on Cholmod
|
---|
460 | *
|
---|
461 | * This class allows to solve for A.X = B sparse linear problems via a supernodal LL^T Cholesky factorization
|
---|
462 | * using the Cholmod library.
|
---|
463 | * This supernodal variant performs best on dense enough problems, e.g., 3D FEM, or very high order 2D FEM.
|
---|
464 | * The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
|
---|
465 | * X and B can be either dense or sparse.
|
---|
466 | *
|
---|
467 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
---|
468 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
---|
469 | * or Upper. Default is Lower.
|
---|
470 | *
|
---|
471 | * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
|
---|
472 | *
|
---|
473 | * \sa \ref TutorialSparseDirectSolvers
|
---|
474 | */
|
---|
475 | template<typename _MatrixType, int _UpLo = Lower>
|
---|
476 | class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
|
---|
477 | {
|
---|
478 | typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
|
---|
479 | using Base::m_cholmod;
|
---|
480 |
|
---|
481 | public:
|
---|
482 |
|
---|
483 | typedef _MatrixType MatrixType;
|
---|
484 |
|
---|
485 | CholmodSupernodalLLT() : Base() { init(); }
|
---|
486 |
|
---|
487 | CholmodSupernodalLLT(const MatrixType& matrix) : Base()
|
---|
488 | {
|
---|
489 | init();
|
---|
490 | Base::compute(matrix);
|
---|
491 | }
|
---|
492 |
|
---|
493 | ~CholmodSupernodalLLT() {}
|
---|
494 | protected:
|
---|
495 | void init()
|
---|
496 | {
|
---|
497 | m_cholmod.final_asis = 1;
|
---|
498 | m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
|
---|
499 | }
|
---|
500 | };
|
---|
501 |
|
---|
502 | /** \ingroup CholmodSupport_Module
|
---|
503 | * \class CholmodDecomposition
|
---|
504 | * \brief A general Cholesky factorization and solver based on Cholmod
|
---|
505 | *
|
---|
506 | * This class allows to solve for A.X = B sparse linear problems via a LL^T or LDL^T Cholesky factorization
|
---|
507 | * using the Cholmod library. The sparse matrix A must be selfadjoint and positive definite. The vectors or matrices
|
---|
508 | * X and B can be either dense or sparse.
|
---|
509 | *
|
---|
510 | * This variant permits to change the underlying Cholesky method at runtime.
|
---|
511 | * On the other hand, it does not provide access to the result of the factorization.
|
---|
512 | * The default is to let Cholmod automatically choose between a simplicial and supernodal factorization.
|
---|
513 | *
|
---|
514 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
---|
515 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
---|
516 | * or Upper. Default is Lower.
|
---|
517 | *
|
---|
518 | * This class supports all kind of SparseMatrix<>: row or column major; upper, lower, or both; compressed or non compressed.
|
---|
519 | *
|
---|
520 | * \sa \ref TutorialSparseDirectSolvers
|
---|
521 | */
|
---|
522 | template<typename _MatrixType, int _UpLo = Lower>
|
---|
523 | class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
|
---|
524 | {
|
---|
525 | typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
|
---|
526 | using Base::m_cholmod;
|
---|
527 |
|
---|
528 | public:
|
---|
529 |
|
---|
530 | typedef _MatrixType MatrixType;
|
---|
531 |
|
---|
532 | CholmodDecomposition() : Base() { init(); }
|
---|
533 |
|
---|
534 | CholmodDecomposition(const MatrixType& matrix) : Base()
|
---|
535 | {
|
---|
536 | init();
|
---|
537 | Base::compute(matrix);
|
---|
538 | }
|
---|
539 |
|
---|
540 | ~CholmodDecomposition() {}
|
---|
541 |
|
---|
542 | void setMode(CholmodMode mode)
|
---|
543 | {
|
---|
544 | switch(mode)
|
---|
545 | {
|
---|
546 | case CholmodAuto:
|
---|
547 | m_cholmod.final_asis = 1;
|
---|
548 | m_cholmod.supernodal = CHOLMOD_AUTO;
|
---|
549 | break;
|
---|
550 | case CholmodSimplicialLLt:
|
---|
551 | m_cholmod.final_asis = 0;
|
---|
552 | m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
|
---|
553 | m_cholmod.final_ll = 1;
|
---|
554 | break;
|
---|
555 | case CholmodSupernodalLLt:
|
---|
556 | m_cholmod.final_asis = 1;
|
---|
557 | m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
|
---|
558 | break;
|
---|
559 | case CholmodLDLt:
|
---|
560 | m_cholmod.final_asis = 1;
|
---|
561 | m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
|
---|
562 | break;
|
---|
563 | default:
|
---|
564 | break;
|
---|
565 | }
|
---|
566 | }
|
---|
567 | protected:
|
---|
568 | void init()
|
---|
569 | {
|
---|
570 | m_cholmod.final_asis = 1;
|
---|
571 | m_cholmod.supernodal = CHOLMOD_AUTO;
|
---|
572 | }
|
---|
573 | };
|
---|
574 |
|
---|
575 | namespace internal {
|
---|
576 |
|
---|
577 | template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
|
---|
578 | struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
|
---|
579 | : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
|
---|
580 | {
|
---|
581 | typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
|
---|
582 | EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
|
---|
583 |
|
---|
584 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
585 | {
|
---|
586 | dec()._solve(rhs(),dst);
|
---|
587 | }
|
---|
588 | };
|
---|
589 |
|
---|
590 | template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
|
---|
591 | struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
|
---|
592 | : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
|
---|
593 | {
|
---|
594 | typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
|
---|
595 | EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
|
---|
596 |
|
---|
597 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
598 | {
|
---|
599 | dec()._solve(rhs(),dst);
|
---|
600 | }
|
---|
601 | };
|
---|
602 |
|
---|
603 | } // end namespace internal
|
---|
604 |
|
---|
605 | } // end namespace Eigen
|
---|
606 |
|
---|
607 | #endif // EIGEN_CHOLMODSUPPORT_H
|
---|