Eigen  5.0.1
 
Loading...
Searching...
No Matches
UmfPackSupport.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2011 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// SPDX-License-Identifier: MPL-2.0
10
11#ifndef EIGEN_UMFPACKSUPPORT_H
12#define EIGEN_UMFPACKSUPPORT_H
13
14// for compatibility with super old version of umfpack,
15// This may not be strictly needed, but it is harmless.
16#ifndef SuiteSparse_long
17#ifdef UF_long
18#define SuiteSparse_long UF_long
19#else
20#error neither SuiteSparse_long nor UF_long are defined
21#endif
22#endif
23
24// IWYU pragma: private
25#include "./InternalHeaderCheck.h"
26
27namespace Eigen {
28
29/* TODO extract L, extract U, compute det, etc... */
30
31// generic double/complex<double> wrapper functions:
32
33// Defaults
34inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, int) { umfpack_di_defaults(control); }
35
36inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, int) {
37 umfpack_zi_defaults(control);
38}
39
40inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, SuiteSparse_long) {
41 umfpack_dl_defaults(control);
42}
43
44inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long) {
45 umfpack_zl_defaults(control);
46}
47
48// Report info
49inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, int) {
50 umfpack_di_report_info(control, info);
51}
52
53inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, int) {
54 umfpack_zi_report_info(control, info);
55}
56
57inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, SuiteSparse_long) {
58 umfpack_dl_report_info(control, info);
59}
60
61inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>,
62 SuiteSparse_long) {
63 umfpack_zl_report_info(control, info);
64}
65
66// Report status
67inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, int) {
68 umfpack_di_report_status(control, status);
69}
70
71inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, int) {
72 umfpack_zi_report_status(control, status);
73}
74
75inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, SuiteSparse_long) {
76 umfpack_dl_report_status(control, status);
77}
78
79inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, SuiteSparse_long) {
80 umfpack_zl_report_status(control, status);
81}
82
83// report control
84inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, int) { umfpack_di_report_control(control); }
85
86inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, int) {
87 umfpack_zi_report_control(control);
88}
89
90inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, SuiteSparse_long) {
91 umfpack_dl_report_control(control);
92}
93
94inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long) {
95 umfpack_zl_report_control(control);
96}
97
98// Free numeric
99inline void umfpack_free_numeric(void **Numeric, double, int) {
100 umfpack_di_free_numeric(Numeric);
101 *Numeric = 0;
102}
103
104inline void umfpack_free_numeric(void **Numeric, std::complex<double>, int) {
105 umfpack_zi_free_numeric(Numeric);
106 *Numeric = 0;
107}
108
109inline void umfpack_free_numeric(void **Numeric, double, SuiteSparse_long) {
110 umfpack_dl_free_numeric(Numeric);
111 *Numeric = 0;
112}
113
114inline void umfpack_free_numeric(void **Numeric, std::complex<double>, SuiteSparse_long) {
115 umfpack_zl_free_numeric(Numeric);
116 *Numeric = 0;
117}
118
119// Free symbolic
120inline void umfpack_free_symbolic(void **Symbolic, double, int) {
121 umfpack_di_free_symbolic(Symbolic);
122 *Symbolic = 0;
123}
124
125inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, int) {
126 umfpack_zi_free_symbolic(Symbolic);
127 *Symbolic = 0;
128}
129
130inline void umfpack_free_symbolic(void **Symbolic, double, SuiteSparse_long) {
131 umfpack_dl_free_symbolic(Symbolic);
132 *Symbolic = 0;
133}
134
135inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, SuiteSparse_long) {
136 umfpack_zl_free_symbolic(Symbolic);
137 *Symbolic = 0;
138}
139
140// Symbolic
141inline int umfpack_symbolic(int n_row, int n_col, const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
142 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
143 return umfpack_di_symbolic(n_row, n_col, Ap, Ai, Ax, Symbolic, Control, Info);
144}
145
146inline int umfpack_symbolic(int n_row, int n_col, const int Ap[], const int Ai[], const std::complex<double> Ax[],
147 void **Symbolic, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
148 return umfpack_zi_symbolic(n_row, n_col, Ap, Ai, &numext::real_ref(Ax[0]), 0, Symbolic, Control, Info);
149}
150inline SuiteSparse_long umfpack_symbolic(SuiteSparse_long n_row, SuiteSparse_long n_col, const SuiteSparse_long Ap[],
151 const SuiteSparse_long Ai[], const double Ax[], void **Symbolic,
152 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
153 return umfpack_dl_symbolic(n_row, n_col, Ap, Ai, Ax, Symbolic, Control, Info);
154}
155
156inline SuiteSparse_long umfpack_symbolic(SuiteSparse_long n_row, SuiteSparse_long n_col, const SuiteSparse_long Ap[],
157 const SuiteSparse_long Ai[], const std::complex<double> Ax[], void **Symbolic,
158 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
159 return umfpack_zl_symbolic(n_row, n_col, Ap, Ai, &numext::real_ref(Ax[0]), 0, Symbolic, Control, Info);
160}
161
162// Numeric
163inline int umfpack_numeric(const int Ap[], const int Ai[], const double Ax[], void *Symbolic, void **Numeric,
164 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
165 return umfpack_di_numeric(Ap, Ai, Ax, Symbolic, Numeric, Control, Info);
166}
167
168inline int umfpack_numeric(const int Ap[], const int Ai[], const std::complex<double> Ax[], void *Symbolic,
169 void **Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
170 return umfpack_zi_numeric(Ap, Ai, &numext::real_ref(Ax[0]), 0, Symbolic, Numeric, Control, Info);
171}
172inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[],
173 void *Symbolic, void **Numeric, const double Control[UMFPACK_CONTROL],
174 double Info[UMFPACK_INFO]) {
175 return umfpack_dl_numeric(Ap, Ai, Ax, Symbolic, Numeric, Control, Info);
176}
177
178inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[],
179 const std::complex<double> Ax[], void *Symbolic, void **Numeric,
180 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
181 return umfpack_zl_numeric(Ap, Ai, &numext::real_ref(Ax[0]), 0, Symbolic, Numeric, Control, Info);
182}
183
184// solve
185inline int umfpack_solve(int sys, const int Ap[], const int Ai[], const double Ax[], double X[], const double B[],
186 void *Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
187 return umfpack_di_solve(sys, Ap, Ai, Ax, X, B, Numeric, Control, Info);
188}
189
190inline int umfpack_solve(int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
191 std::complex<double> X[], const std::complex<double> B[], void *Numeric,
192 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
193 return umfpack_zi_solve(sys, Ap, Ai, &numext::real_ref(Ax[0]), 0, &numext::real_ref(X[0]), 0, &numext::real_ref(B[0]),
194 0, Numeric, Control, Info);
195}
196
197inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[],
198 const double Ax[], double X[], const double B[], void *Numeric,
199 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
200 return umfpack_dl_solve(sys, Ap, Ai, Ax, X, B, Numeric, Control, Info);
201}
202
203inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[],
204 const std::complex<double> Ax[], std::complex<double> X[],
205 const std::complex<double> B[], void *Numeric,
206 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO]) {
207 return umfpack_zl_solve(sys, Ap, Ai, &numext::real_ref(Ax[0]), 0, &numext::real_ref(X[0]), 0, &numext::real_ref(B[0]),
208 0, Numeric, Control, Info);
209}
210
211// Get Lunz
212inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double) {
213 return umfpack_di_get_lunz(lnz, unz, n_row, n_col, nz_udiag, Numeric);
214}
215
216inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric,
217 std::complex<double>) {
218 return umfpack_zi_get_lunz(lnz, unz, n_row, n_col, nz_udiag, Numeric);
219}
220
221inline SuiteSparse_long umfpack_get_lunz(SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row,
222 SuiteSparse_long *n_col, SuiteSparse_long *nz_udiag, void *Numeric, double) {
223 return umfpack_dl_get_lunz(lnz, unz, n_row, n_col, nz_udiag, Numeric);
224}
225
226inline SuiteSparse_long umfpack_get_lunz(SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row,
227 SuiteSparse_long *n_col, SuiteSparse_long *nz_udiag, void *Numeric,
228 std::complex<double>) {
229 return umfpack_zl_get_lunz(lnz, unz, n_row, n_col, nz_udiag, Numeric);
230}
231
232// Get Numeric
233inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[], int P[], int Q[],
234 double Dx[], int *do_recip, double Rs[], void *Numeric) {
235 return umfpack_di_get_numeric(Lp, Lj, Lx, Up, Ui, Ux, P, Q, Dx, do_recip, Rs, Numeric);
236}
237
238inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[],
239 std::complex<double> Ux[], int P[], int Q[], std::complex<double> Dx[], int *do_recip,
240 double Rs[], void *Numeric) {
241 double &lx0_real = numext::real_ref(Lx[0]);
242 double &ux0_real = numext::real_ref(Ux[0]);
243 double &dx0_real = numext::real_ref(Dx[0]);
244 return umfpack_zi_get_numeric(Lp, Lj, Lx ? &lx0_real : 0, 0, Up, Ui, Ux ? &ux0_real : 0, 0, P, Q, Dx ? &dx0_real : 0,
245 0, do_recip, Rs, Numeric);
246}
247inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], double Lx[],
248 SuiteSparse_long Up[], SuiteSparse_long Ui[], double Ux[],
249 SuiteSparse_long P[], SuiteSparse_long Q[], double Dx[],
250 SuiteSparse_long *do_recip, double Rs[], void *Numeric) {
251 return umfpack_dl_get_numeric(Lp, Lj, Lx, Up, Ui, Ux, P, Q, Dx, do_recip, Rs, Numeric);
252}
253
254inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], std::complex<double> Lx[],
255 SuiteSparse_long Up[], SuiteSparse_long Ui[], std::complex<double> Ux[],
256 SuiteSparse_long P[], SuiteSparse_long Q[], std::complex<double> Dx[],
257 SuiteSparse_long *do_recip, double Rs[], void *Numeric) {
258 double &lx0_real = numext::real_ref(Lx[0]);
259 double &ux0_real = numext::real_ref(Ux[0]);
260 double &dx0_real = numext::real_ref(Dx[0]);
261 return umfpack_zl_get_numeric(Lp, Lj, Lx ? &lx0_real : 0, 0, Up, Ui, Ux ? &ux0_real : 0, 0, P, Q, Dx ? &dx0_real : 0,
262 0, do_recip, Rs, Numeric);
263}
264
265// Get Determinant
266inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info[UMFPACK_INFO], int) {
267 return umfpack_di_get_determinant(Mx, Ex, NumericHandle, User_Info);
268}
269
270inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle,
271 double User_Info[UMFPACK_INFO], int) {
272 double &mx_real = numext::real_ref(*Mx);
273 return umfpack_zi_get_determinant(&mx_real, 0, Ex, NumericHandle, User_Info);
274}
275
276inline SuiteSparse_long umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle,
277 double User_Info[UMFPACK_INFO], SuiteSparse_long) {
278 return umfpack_dl_get_determinant(Mx, Ex, NumericHandle, User_Info);
279}
280
281inline SuiteSparse_long umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle,
282 double User_Info[UMFPACK_INFO], SuiteSparse_long) {
283 double &mx_real = numext::real_ref(*Mx);
284 return umfpack_zl_get_determinant(&mx_real, 0, Ex, NumericHandle, User_Info);
285}
286
302template <typename MatrixType_>
303class UmfPackLU : public SparseSolverBase<UmfPackLU<MatrixType_> > {
304 protected:
305 typedef SparseSolverBase<UmfPackLU<MatrixType_> > Base;
306 using Base::m_isInitialized;
307
308 public:
309 using Base::_solve_impl;
310 typedef MatrixType_ MatrixType;
311 typedef typename MatrixType::Scalar Scalar;
312 typedef typename MatrixType::RealScalar RealScalar;
313 typedef typename MatrixType::StorageIndex StorageIndex;
314 typedef Matrix<Scalar, Dynamic, 1> Vector;
315 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
316 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
317 typedef SparseMatrix<Scalar> LUMatrixType;
318 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> UmfpackMatrixType;
319 typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef;
320 enum { ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime };
321
322 public:
323 typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl;
324 typedef Array<double, UMFPACK_INFO, 1> UmfpackInfo;
325
326 UmfPackLU() : m_dummy(0, 0), mp_matrix(m_dummy) { init(); }
327
328 template <typename InputMatrixType>
329 explicit UmfPackLU(const InputMatrixType &matrix) : mp_matrix(matrix) {
330 init();
331 compute(matrix);
332 }
333
334 ~UmfPackLU() {
335 if (m_symbolic) umfpack_free_symbolic(&m_symbolic, Scalar(), StorageIndex());
336 if (m_numeric) umfpack_free_numeric(&m_numeric, Scalar(), StorageIndex());
337 }
338
339 inline Index rows() const { return mp_matrix.rows(); }
340 inline Index cols() const { return mp_matrix.cols(); }
341
347 ComputationInfo info() const {
348 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
349 return m_info;
350 }
351
352 inline const LUMatrixType &matrixL() const {
353 if (m_extractedDataAreDirty) extractData();
354 return m_l;
355 }
356
357 inline const LUMatrixType &matrixU() const {
358 if (m_extractedDataAreDirty) extractData();
359 return m_u;
360 }
361
362 inline const IntColVectorType &permutationP() const {
363 if (m_extractedDataAreDirty) extractData();
364 return m_p;
365 }
366
367 inline const IntRowVectorType &permutationQ() const {
368 if (m_extractedDataAreDirty) extractData();
369 return m_q;
370 }
371
376 template <typename InputMatrixType>
377 void compute(const InputMatrixType &matrix) {
378 if (m_symbolic) umfpack_free_symbolic(&m_symbolic, Scalar(), StorageIndex());
379 if (m_numeric) umfpack_free_numeric(&m_numeric, Scalar(), StorageIndex());
380 grab(matrix.derived());
381 analyzePattern_impl();
382 factorize_impl();
383 }
384
391 template <typename InputMatrixType>
392 void analyzePattern(const InputMatrixType &matrix) {
393 if (m_symbolic) umfpack_free_symbolic(&m_symbolic, Scalar(), StorageIndex());
394 if (m_numeric) umfpack_free_numeric(&m_numeric, Scalar(), StorageIndex());
395
396 grab(matrix.derived());
397
398 analyzePattern_impl();
399 }
400
406 inline int umfpackFactorizeReturncode() const {
407 eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()");
408 return m_fact_errorCode;
409 }
410
417 inline const UmfpackControl &umfpackControl() const { return m_control; }
418
425 inline UmfpackControl &umfpackControl() { return m_control; }
426
433 template <typename InputMatrixType>
434 void factorize(const InputMatrixType &matrix) {
435 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
436 if (m_numeric) umfpack_free_numeric(&m_numeric, Scalar(), StorageIndex());
437
438 grab(matrix.derived());
439
440 factorize_impl();
441 }
442
447 void printUmfpackControl() { umfpack_report_control(m_control.data(), Scalar(), StorageIndex()); }
448
453 void printUmfpackInfo() {
454 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
455 umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar(), StorageIndex());
456 }
457
463 void printUmfpackStatus() {
464 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
465 umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar(), StorageIndex());
466 }
467
469 template <typename BDerived, typename XDerived>
470 bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
471
472 Scalar determinant() const;
473
474 void extractData() const;
475
476 protected:
477 void init() {
478 m_info = InvalidInput;
479 m_isInitialized = false;
480 m_numeric = 0;
481 m_symbolic = 0;
482 m_extractedDataAreDirty = true;
483
484 umfpack_defaults(m_control.data(), Scalar(), StorageIndex());
485 }
486
487 void analyzePattern_impl() {
488 m_fact_errorCode = umfpack_symbolic(internal::convert_index<StorageIndex>(mp_matrix.rows()),
489 internal::convert_index<StorageIndex>(mp_matrix.cols()),
490 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
491 &m_symbolic, m_control.data(), m_umfpackInfo.data());
492
493 m_isInitialized = true;
494 m_info = m_fact_errorCode ? InvalidInput : Success;
495 m_analysisIsOk = true;
496 m_factorizationIsOk = false;
497 m_extractedDataAreDirty = true;
498 }
499
500 void factorize_impl() {
501 m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
502 m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data());
503
504 m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
505 m_factorizationIsOk = true;
506 m_extractedDataAreDirty = true;
507 }
508
509 template <typename MatrixDerived>
510 void grab(const EigenBase<MatrixDerived> &A) {
511 internal::destroy_at(&mp_matrix);
512 internal::construct_at(&mp_matrix, A.derived());
513 }
514
515 void grab(const UmfpackMatrixRef &A) {
516 if (&(A.derived()) != &mp_matrix) {
517 internal::destroy_at(&mp_matrix);
518 internal::construct_at(&mp_matrix, A);
519 }
520 }
521
522 // cached data to reduce reallocation, etc.
523 mutable LUMatrixType m_l;
524 StorageIndex m_fact_errorCode;
525 UmfpackControl m_control;
526 mutable UmfpackInfo m_umfpackInfo;
527
528 mutable LUMatrixType m_u;
529 mutable IntColVectorType m_p;
530 mutable IntRowVectorType m_q;
531
532 UmfpackMatrixType m_dummy;
533 UmfpackMatrixRef mp_matrix;
534
535 void *m_numeric;
536 void *m_symbolic;
537
538 mutable ComputationInfo m_info;
539 int m_factorizationIsOk;
540 int m_analysisIsOk;
541 mutable bool m_extractedDataAreDirty;
542
543 private:
544 UmfPackLU(const UmfPackLU &) {}
545};
546
547template <typename MatrixType>
548void UmfPackLU<MatrixType>::extractData() const {
549 if (m_extractedDataAreDirty) {
550 // get size of the data
551 StorageIndex lnz, unz, rows, cols, nz_udiag;
552 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
553
554 // allocate data
555 m_l.resize(rows, (std::min)(rows, cols));
556 m_l.resizeNonZeros(lnz);
557
558 m_u.resize((std::min)(rows, cols), cols);
559 m_u.resizeNonZeros(unz);
560
561 m_p.resize(rows);
562 m_q.resize(cols);
563
564 // extract
565 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(), m_u.outerIndexPtr(),
566 m_u.innerIndexPtr(), m_u.valuePtr(), m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
567
568 m_extractedDataAreDirty = false;
569 }
570}
571
572template <typename MatrixType>
573typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const {
574 Scalar det;
575 umfpack_get_determinant(&det, 0, m_numeric, 0, StorageIndex());
576 return det;
577}
578
579template <typename MatrixType>
580template <typename BDerived, typename XDerived>
581bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const {
582 Index rhsCols = b.cols();
583 eigen_assert((BDerived::Flags & RowMajorBit) == 0 && "UmfPackLU backend does not support non col-major rhs yet");
584 eigen_assert((XDerived::Flags & RowMajorBit) == 0 && "UmfPackLU backend does not support non col-major result yet");
585 eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
586
587 Scalar *x_ptr = 0;
588 Matrix<Scalar, Dynamic, 1> x_tmp;
589 if (x.innerStride() != 1) {
590 x_tmp.resize(x.rows());
591 x_ptr = x_tmp.data();
592 }
593 for (int j = 0; j < rhsCols; ++j) {
594 if (x.innerStride() == 1) x_ptr = &x.col(j).coeffRef(0);
595 StorageIndex errorCode =
596 umfpack_solve(UMFPACK_A, mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(), x_ptr,
597 &b.const_cast_derived().col(j).coeffRef(0), m_numeric, m_control.data(), m_umfpackInfo.data());
598 if (x.innerStride() != 1) x.col(j) = x_tmp;
599 if (errorCode != 0) return false;
600 }
601
602 return true;
603}
604
605} // end namespace Eigen
606
607#endif // EIGEN_UMFPACKSUPPORT_H
Eigen::Array
General-purpose arrays with easy API for coefficient-wise operations.
Definition Array.h:55
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:53
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:188
Eigen::PlainObjectBase::data
constexpr const Scalar * data() const
Definition PlainObjectBase.h:261
Eigen::Ref
A matrix or vector expression mapping an existing expression.
Definition Ref.h:262
Eigen::SparseMatrix
A versatile sparse matrix representation.
Definition SparseMatrix.h:122
Eigen::UmfPackLU::printUmfpackControl
void printUmfpackControl()
Definition UmfPackSupport.h:447
Eigen::UmfPackLU::compute
void compute(const InputMatrixType &matrix)
Definition UmfPackSupport.h:377
Eigen::UmfPackLU::factorize
void factorize(const InputMatrixType &matrix)
Definition UmfPackSupport.h:434
Eigen::UmfPackLU::umfpackFactorizeReturncode
int umfpackFactorizeReturncode() const
Definition UmfPackSupport.h:406
Eigen::UmfPackLU::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition UmfPackSupport.h:347
Eigen::UmfPackLU::umfpackControl
UmfpackControl & umfpackControl()
Definition UmfPackSupport.h:425
Eigen::UmfPackLU::printUmfpackInfo
void printUmfpackInfo()
Definition UmfPackSupport.h:453
Eigen::UmfPackLU::umfpackControl
const UmfpackControl & umfpackControl() const
Definition UmfPackSupport.h:417
Eigen::UmfPackLU::analyzePattern
void analyzePattern(const InputMatrixType &matrix)
Definition UmfPackSupport.h:392
Eigen::UmfPackLU::printUmfpackStatus
void printUmfpackStatus()
Definition UmfPackSupport.h:463
Eigen::ComputationInfo
ComputationInfo
Definition Constants.h:455
Eigen::NumericalIssue
@ NumericalIssue
Definition Constants.h:459
Eigen::InvalidInput
@ InvalidInput
Definition Constants.h:464
Eigen::Success
@ Success
Definition Constants.h:457
Eigen::RowMajorBit
constexpr unsigned int RowMajorBit
Definition Constants.h:71