Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
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
10#ifndef EIGEN_UMFPACKSUPPORT_H
11#define EIGEN_UMFPACKSUPPORT_H
12
13// for compatibility with super old version of umfpack,
14// not sure this is really needed, but this is harmless.
15#ifndef SuiteSparse_long
16#ifdef UF_long
17#define SuiteSparse_long UF_long
18#else
19#error neither SuiteSparse_long nor UF_long are defined
20#endif
21#endif
22
23namespace Eigen {
24
25/* TODO extract L, extract U, compute det, etc... */
26
27// generic double/complex<double> wrapper functions:
28
29
30 // Defaults
31inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, int)
32{ umfpack_di_defaults(control); }
33
34inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, int)
35{ umfpack_zi_defaults(control); }
36
37inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, SuiteSparse_long)
38{ umfpack_dl_defaults(control); }
39
40inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long)
41{ umfpack_zl_defaults(control); }
42
43// Report info
44inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, int)
45{ umfpack_di_report_info(control, info);}
46
47inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, int)
48{ umfpack_zi_report_info(control, info);}
49
50inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, SuiteSparse_long)
51{ umfpack_dl_report_info(control, info);}
52
53inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, SuiteSparse_long)
54{ umfpack_zl_report_info(control, info);}
55
56// Report status
57inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, int)
58{ umfpack_di_report_status(control, status);}
59
60inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, int)
61{ umfpack_zi_report_status(control, status);}
62
63inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, SuiteSparse_long)
64{ umfpack_dl_report_status(control, status);}
65
66inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, SuiteSparse_long)
67{ umfpack_zl_report_status(control, status);}
68
69// report control
70inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, int)
71{ umfpack_di_report_control(control);}
72
73inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, int)
74{ umfpack_zi_report_control(control);}
75
76inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, SuiteSparse_long)
77{ umfpack_dl_report_control(control);}
78
79inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long)
80{ umfpack_zl_report_control(control);}
81
82// Free numeric
83inline void umfpack_free_numeric(void **Numeric, double, int)
84{ umfpack_di_free_numeric(Numeric); *Numeric = 0; }
85
86inline void umfpack_free_numeric(void **Numeric, std::complex<double>, int)
87{ umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
88
89inline void umfpack_free_numeric(void **Numeric, double, SuiteSparse_long)
90{ umfpack_dl_free_numeric(Numeric); *Numeric = 0; }
91
92inline void umfpack_free_numeric(void **Numeric, std::complex<double>, SuiteSparse_long)
93{ umfpack_zl_free_numeric(Numeric); *Numeric = 0; }
94
95// Free symbolic
96inline void umfpack_free_symbolic(void **Symbolic, double, int)
97{ umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
98
99inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, int)
100{ umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
101
102inline void umfpack_free_symbolic(void **Symbolic, double, SuiteSparse_long)
103{ umfpack_dl_free_symbolic(Symbolic); *Symbolic = 0; }
104
105inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, SuiteSparse_long)
106{ umfpack_zl_free_symbolic(Symbolic); *Symbolic = 0; }
107
108// Symbolic
109inline int umfpack_symbolic(int n_row,int n_col,
110 const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
111 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
112{
113 return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
114}
115
116inline int umfpack_symbolic(int n_row,int n_col,
117 const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
118 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
119{
120 return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
121}
122inline SuiteSparse_long umfpack_symbolic( SuiteSparse_long n_row,SuiteSparse_long n_col,
123 const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[], void **Symbolic,
124 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
125{
126 return umfpack_dl_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
127}
128
129inline SuiteSparse_long umfpack_symbolic( SuiteSparse_long n_row,SuiteSparse_long n_col,
130 const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[], void **Symbolic,
131 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
132{
133 return umfpack_zl_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
134}
135
136// Numeric
137inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
138 void *Symbolic, void **Numeric,
139 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
140{
141 return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
142}
143
144inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
145 void *Symbolic, void **Numeric,
146 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
147{
148 return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
149}
150inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[],
151 void *Symbolic, void **Numeric,
152 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
153{
154 return umfpack_dl_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
155}
156
157inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[],
158 void *Symbolic, void **Numeric,
159 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
160{
161 return umfpack_zl_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
162}
163
164// solve
165inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
166 double X[], const double B[], void *Numeric,
167 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
168{
169 return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
170}
171
172inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
173 std::complex<double> X[], const std::complex<double> B[], void *Numeric,
174 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
175{
176 return umfpack_zi_solve(sys,Ap,Ai,&numext::real_ref(Ax[0]),0,&numext::real_ref(X[0]),0,&numext::real_ref(B[0]),0,Numeric,Control,Info);
177}
178
179inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[],
180 double X[], const double B[], void *Numeric,
181 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
182{
183 return umfpack_dl_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
184}
185
186inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[],
187 std::complex<double> X[], const std::complex<double> B[], void *Numeric,
188 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
189{
190 return umfpack_zl_solve(sys,Ap,Ai,&numext::real_ref(Ax[0]),0,&numext::real_ref(X[0]),0,&numext::real_ref(B[0]),0,Numeric,Control,Info);
191}
192
193// Get Lunz
194inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
195{
196 return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
197}
198
199inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
200{
201 return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
202}
203
204inline SuiteSparse_long umfpack_get_lunz( SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row, SuiteSparse_long *n_col,
205 SuiteSparse_long *nz_udiag, void *Numeric, double)
206{
207 return umfpack_dl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
208}
209
210inline SuiteSparse_long umfpack_get_lunz( SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row, SuiteSparse_long *n_col,
211 SuiteSparse_long *nz_udiag, void *Numeric, std::complex<double>)
212{
213 return umfpack_zl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
214}
215
216// Get Numeric
217inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
218 int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
219{
220 return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
221}
222
223inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
224 int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
225{
226 double& lx0_real = numext::real_ref(Lx[0]);
227 double& ux0_real = numext::real_ref(Ux[0]);
228 double& dx0_real = numext::real_ref(Dx[0]);
229 return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
230 Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
231}
232inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], double Lx[], SuiteSparse_long Up[], SuiteSparse_long Ui[], double Ux[],
233 SuiteSparse_long P[], SuiteSparse_long Q[], double Dx[], SuiteSparse_long *do_recip, double Rs[], void *Numeric)
234{
235 return umfpack_dl_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
236}
237
238inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], std::complex<double> Lx[], SuiteSparse_long Up[], SuiteSparse_long Ui[], std::complex<double> Ux[],
239 SuiteSparse_long P[], SuiteSparse_long Q[], std::complex<double> Dx[], SuiteSparse_long *do_recip, double Rs[], void *Numeric)
240{
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_zl_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
245 Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
246}
247
248// Get Determinant
249inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int)
250{
251 return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
252}
253
254inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int)
255{
256 double& mx_real = numext::real_ref(*Mx);
257 return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
258}
259
260inline SuiteSparse_long umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], SuiteSparse_long)
261{
262 return umfpack_dl_get_determinant(Mx,Ex,NumericHandle,User_Info);
263}
264
265inline SuiteSparse_long umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], SuiteSparse_long)
266{
267 double& mx_real = numext::real_ref(*Mx);
268 return umfpack_zl_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
269}
270
271
287template<typename _MatrixType>
288class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
289{
290 protected:
292 using Base::m_isInitialized;
293 public:
294 using Base::_solve_impl;
295 typedef _MatrixType MatrixType;
296 typedef typename MatrixType::Scalar Scalar;
297 typedef typename MatrixType::RealScalar RealScalar;
298 typedef typename MatrixType::StorageIndex StorageIndex;
305 enum {
306 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
307 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
308 };
309
310 public:
311
314
315 UmfPackLU()
316 : m_dummy(0,0), mp_matrix(m_dummy)
317 {
318 init();
319 }
320
321 template<typename InputMatrixType>
322 explicit UmfPackLU(const InputMatrixType& matrix)
323 : mp_matrix(matrix)
324 {
325 init();
326 compute(matrix);
327 }
328
329 ~UmfPackLU()
330 {
331 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(), StorageIndex());
332 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(), StorageIndex());
333 }
334
335 inline Index rows() const { return mp_matrix.rows(); }
336 inline Index cols() const { return mp_matrix.cols(); }
337
344 {
345 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
346 return m_info;
347 }
348
349 inline const LUMatrixType& matrixL() const
350 {
351 if (m_extractedDataAreDirty) extractData();
352 return m_l;
353 }
354
355 inline const LUMatrixType& matrixU() const
356 {
357 if (m_extractedDataAreDirty) extractData();
358 return m_u;
359 }
360
361 inline const IntColVectorType& permutationP() const
362 {
363 if (m_extractedDataAreDirty) extractData();
364 return m_p;
365 }
366
367 inline const IntRowVectorType& permutationQ() const
368 {
369 if (m_extractedDataAreDirty) extractData();
370 return m_q;
371 }
372
377 template<typename InputMatrixType>
378 void compute(const InputMatrixType& matrix)
379 {
380 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
381 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
382 grab(matrix.derived());
383 analyzePattern_impl();
384 factorize_impl();
385 }
386
393 template<typename InputMatrixType>
394 void analyzePattern(const InputMatrixType& matrix)
395 {
396 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
397 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
398
399 grab(matrix.derived());
400
401 analyzePattern_impl();
402 }
403
409 inline int umfpackFactorizeReturncode() const
410 {
411 eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()");
412 return m_fact_errorCode;
413 }
414
421 inline const UmfpackControl& umfpackControl() const
422 {
423 return m_control;
424 }
425
433 {
434 return m_control;
435 }
436
443 template<typename InputMatrixType>
444 void factorize(const InputMatrixType& matrix)
445 {
446 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
447 if(m_numeric)
448 umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
449
450 grab(matrix.derived());
451
452 factorize_impl();
453 }
454
460 {
461 umfpack_report_control(m_control.data(), Scalar(),StorageIndex());
462 }
463
469 {
470 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
471 umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar(),StorageIndex());
472 }
473
479 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
480 umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar(),StorageIndex());
481 }
482
484 template<typename BDerived,typename XDerived>
485 bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
486
487 Scalar determinant() const;
488
489 void extractData() const;
490
491 protected:
492
493 void init()
494 {
495 m_info = InvalidInput;
496 m_isInitialized = false;
497 m_numeric = 0;
498 m_symbolic = 0;
499 m_extractedDataAreDirty = true;
500
501 umfpack_defaults(m_control.data(), Scalar(),StorageIndex());
502 }
503
504 void analyzePattern_impl()
505 {
506 m_fact_errorCode = umfpack_symbolic(internal::convert_index<StorageIndex>(mp_matrix.rows()),
507 internal::convert_index<StorageIndex>(mp_matrix.cols()),
508 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
509 &m_symbolic, m_control.data(), m_umfpackInfo.data());
510
511 m_isInitialized = true;
512 m_info = m_fact_errorCode ? InvalidInput : Success;
513 m_analysisIsOk = true;
514 m_factorizationIsOk = false;
515 m_extractedDataAreDirty = true;
516 }
517
518 void factorize_impl()
519 {
520
521 m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
522 m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data());
523
524 m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
525 m_factorizationIsOk = true;
526 m_extractedDataAreDirty = true;
527 }
528
529 template<typename MatrixDerived>
530 void grab(const EigenBase<MatrixDerived> &A)
531 {
532 mp_matrix.~UmfpackMatrixRef();
533 ::new (&mp_matrix) UmfpackMatrixRef(A.derived());
534 }
535
536 void grab(const UmfpackMatrixRef &A)
537 {
538 if(&(A.derived()) != &mp_matrix)
539 {
540 mp_matrix.~UmfpackMatrixRef();
541 ::new (&mp_matrix) UmfpackMatrixRef(A);
542 }
543 }
544
545 // cached data to reduce reallocation, etc.
546 mutable LUMatrixType m_l;
547 StorageIndex m_fact_errorCode;
548 UmfpackControl m_control;
549 mutable UmfpackInfo m_umfpackInfo;
550
551 mutable LUMatrixType m_u;
552 mutable IntColVectorType m_p;
553 mutable IntRowVectorType m_q;
554
555 UmfpackMatrixType m_dummy;
556 UmfpackMatrixRef mp_matrix;
557
558 void* m_numeric;
559 void* m_symbolic;
560
561 mutable ComputationInfo m_info;
562 int m_factorizationIsOk;
563 int m_analysisIsOk;
564 mutable bool m_extractedDataAreDirty;
565
566 private:
567 UmfPackLU(const UmfPackLU& ) { }
568};
569
570
571template<typename MatrixType>
572void UmfPackLU<MatrixType>::extractData() const
573{
574 if (m_extractedDataAreDirty)
575 {
576 // get size of the data
577 StorageIndex lnz, unz, rows, cols, nz_udiag;
578 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
579
580 // allocate data
581 m_l.resize(rows,(std::min)(rows,cols));
582 m_l.resizeNonZeros(lnz);
583
584 m_u.resize((std::min)(rows,cols),cols);
585 m_u.resizeNonZeros(unz);
586
587 m_p.resize(rows);
588 m_q.resize(cols);
589
590 // extract
591 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
592 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
593 m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
594
595 m_extractedDataAreDirty = false;
596 }
597}
598
599template<typename MatrixType>
600typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const
601{
602 Scalar det;
603 umfpack_get_determinant(&det, 0, m_numeric, 0, StorageIndex());
604 return det;
605}
606
607template<typename MatrixType>
608template<typename BDerived,typename XDerived>
609bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
610{
611 Index rhsCols = b.cols();
612 eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
613 eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
614 eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
615
616 Scalar* x_ptr = 0;
617 Matrix<Scalar,Dynamic,1> x_tmp;
618 if(x.innerStride()!=1)
619 {
620 x_tmp.resize(x.rows());
621 x_ptr = x_tmp.data();
622 }
623 for (int j=0; j<rhsCols; ++j)
624 {
625 if(x.innerStride()==1)
626 x_ptr = &x.col(j).coeffRef(0);
627 StorageIndex errorCode = umfpack_solve(UMFPACK_A,
628 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
629 x_ptr, &b.const_cast_derived().col(j).coeffRef(0),
630 m_numeric, m_control.data(), m_umfpackInfo.data());
631 if(x.innerStride()!=1)
632 x.col(j) = x_tmp;
633 if (errorCode!=0)
634 return false;
635 }
636
637 return true;
638}
639
640} // end namespace Eigen
641
642#endif // EIGEN_UMFPACKSUPPORT_H
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:47
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
const Scalar * data() const
Definition: PlainObjectBase.h:247
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:283
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
A sparse LU factorization and solver based on UmfPack.
Definition: UmfPackSupport.h:289
void compute(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:378
void factorize(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:444
void printUmfpackControl()
Definition: UmfPackSupport.h:459
const UmfpackControl & umfpackControl() const
Definition: UmfPackSupport.h:421
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: UmfPackSupport.h:343
UmfpackControl & umfpackControl()
Definition: UmfPackSupport.h:432
int umfpackFactorizeReturncode() const
Definition: UmfPackSupport.h:409
void printUmfpackInfo()
Definition: UmfPackSupport.h:468
void analyzePattern(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:394
void printUmfpackStatus()
Definition: UmfPackSupport.h:478
ComputationInfo
Definition: Constants.h:440
@ NumericalIssue
Definition: Constants.h:444
@ InvalidInput
Definition: Constants.h:449
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74