Eigen  3.3.9
 
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
13namespace Eigen {
14
15/* TODO extract L, extract U, compute det, etc... */
16
17// generic double/complex<double> wrapper functions:
18
19
20inline void umfpack_defaults(double control[UMFPACK_CONTROL], double)
21{ umfpack_di_defaults(control); }
22
23inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>)
24{ umfpack_zi_defaults(control); }
25
26inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double)
27{ umfpack_di_report_info(control, info);}
28
29inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>)
30{ umfpack_zi_report_info(control, info);}
31
32inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double)
33{ umfpack_di_report_status(control, status);}
34
35inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>)
36{ umfpack_zi_report_status(control, status);}
37
38inline void umfpack_report_control(double control[UMFPACK_CONTROL], double)
39{ umfpack_di_report_control(control);}
40
41inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>)
42{ umfpack_zi_report_control(control);}
43
44inline void umfpack_free_numeric(void **Numeric, double)
45{ umfpack_di_free_numeric(Numeric); *Numeric = 0; }
46
47inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
48{ umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
49
50inline void umfpack_free_symbolic(void **Symbolic, double)
51{ umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
52
53inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
54{ umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
55
56inline int umfpack_symbolic(int n_row,int n_col,
57 const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
58 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
59{
60 return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
61}
62
63inline int umfpack_symbolic(int n_row,int n_col,
64 const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
65 const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
66{
67 return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
68}
69
70inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
71 void *Symbolic, void **Numeric,
72 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
73{
74 return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
75}
76
77inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
78 void *Symbolic, void **Numeric,
79 const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
80{
81 return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
82}
83
84inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
85 double X[], const double B[], void *Numeric,
86 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
87{
88 return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
89}
90
91inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
92 std::complex<double> X[], const std::complex<double> B[], void *Numeric,
93 const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
94{
95 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);
96}
97
98inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
99{
100 return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
101}
102
103inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
104{
105 return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
106}
107
108inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
109 int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
110{
111 return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
112}
113
114inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
115 int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
116{
117 double& lx0_real = numext::real_ref(Lx[0]);
118 double& ux0_real = numext::real_ref(Ux[0]);
119 double& dx0_real = numext::real_ref(Dx[0]);
120 return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
121 Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
122}
123
124inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
125{
126 return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
127}
128
129inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
130{
131 double& mx_real = numext::real_ref(*Mx);
132 return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
133}
134
135
151template<typename _MatrixType>
152class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
153{
154 protected:
155 typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base;
156 using Base::m_isInitialized;
157 public:
158 using Base::_solve_impl;
159 typedef _MatrixType MatrixType;
160 typedef typename MatrixType::Scalar Scalar;
161 typedef typename MatrixType::RealScalar RealScalar;
162 typedef typename MatrixType::StorageIndex StorageIndex;
163 typedef Matrix<Scalar,Dynamic,1> Vector;
164 typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
165 typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
166 typedef SparseMatrix<Scalar> LUMatrixType;
167 typedef SparseMatrix<Scalar,ColMajor,int> UmfpackMatrixType;
168 typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef;
169 enum {
170 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
171 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
172 };
173
174 public:
175
176 typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl;
177 typedef Array<double, UMFPACK_INFO, 1> UmfpackInfo;
178
179 UmfPackLU()
180 : m_dummy(0,0), mp_matrix(m_dummy)
181 {
182 init();
183 }
184
185 template<typename InputMatrixType>
186 explicit UmfPackLU(const InputMatrixType& matrix)
187 : mp_matrix(matrix)
188 {
189 init();
190 compute(matrix);
191 }
192
193 ~UmfPackLU()
194 {
195 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
196 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
197 }
198
199 inline Index rows() const { return mp_matrix.rows(); }
200 inline Index cols() const { return mp_matrix.cols(); }
201
207 ComputationInfo info() const
208 {
209 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
210 return m_info;
211 }
212
213 inline const LUMatrixType& matrixL() const
214 {
215 if (m_extractedDataAreDirty) extractData();
216 return m_l;
217 }
218
219 inline const LUMatrixType& matrixU() const
220 {
221 if (m_extractedDataAreDirty) extractData();
222 return m_u;
223 }
224
225 inline const IntColVectorType& permutationP() const
226 {
227 if (m_extractedDataAreDirty) extractData();
228 return m_p;
229 }
230
231 inline const IntRowVectorType& permutationQ() const
232 {
233 if (m_extractedDataAreDirty) extractData();
234 return m_q;
235 }
236
241 template<typename InputMatrixType>
242 void compute(const InputMatrixType& matrix)
243 {
244 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
245 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
246 grab(matrix.derived());
247 analyzePattern_impl();
248 factorize_impl();
249 }
250
257 template<typename InputMatrixType>
258 void analyzePattern(const InputMatrixType& matrix)
259 {
260 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
261 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
262
263 grab(matrix.derived());
264
265 analyzePattern_impl();
266 }
267
273 inline int umfpackFactorizeReturncode() const
274 {
275 eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()");
276 return m_fact_errorCode;
277 }
278
285 inline const UmfpackControl& umfpackControl() const
286 {
287 return m_control;
288 }
289
296 inline UmfpackControl& umfpackControl()
297 {
298 return m_control;
299 }
300
307 template<typename InputMatrixType>
308 void factorize(const InputMatrixType& matrix)
309 {
310 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
311 if(m_numeric)
312 umfpack_free_numeric(&m_numeric,Scalar());
313
314 grab(matrix.derived());
315
316 factorize_impl();
317 }
318
323 void umfpackReportControl()
324 {
325 umfpack_report_control(m_control.data(), Scalar());
326 }
327
332 void umfpackReportInfo()
333 {
334 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
335 umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar());
336 }
337
342 void umfpackReportStatus() {
343 eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
344 umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar());
345 }
346
348 template<typename BDerived,typename XDerived>
349 bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
350
351 Scalar determinant() const;
352
353 void extractData() const;
354
355 protected:
356
357 void init()
358 {
359 m_info = InvalidInput;
360 m_isInitialized = false;
361 m_numeric = 0;
362 m_symbolic = 0;
363 m_extractedDataAreDirty = true;
364
365 umfpack_defaults(m_control.data(), Scalar());
366 }
367
368 void analyzePattern_impl()
369 {
370 m_fact_errorCode = umfpack_symbolic(internal::convert_index<int>(mp_matrix.rows()),
371 internal::convert_index<int>(mp_matrix.cols()),
372 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
373 &m_symbolic, m_control.data(), m_umfpackInfo.data());
374
375 m_isInitialized = true;
376 m_info = m_fact_errorCode ? InvalidInput : Success;
377 m_analysisIsOk = true;
378 m_factorizationIsOk = false;
379 m_extractedDataAreDirty = true;
380 }
381
382 void factorize_impl()
383 {
384
385 m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
386 m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data());
387
388 m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
389 m_factorizationIsOk = true;
390 m_extractedDataAreDirty = true;
391 }
392
393 template<typename MatrixDerived>
394 void grab(const EigenBase<MatrixDerived> &A)
395 {
396 mp_matrix.~UmfpackMatrixRef();
397 ::new (&mp_matrix) UmfpackMatrixRef(A.derived());
398 }
399
400 void grab(const UmfpackMatrixRef &A)
401 {
402 if(&(A.derived()) != &mp_matrix)
403 {
404 mp_matrix.~UmfpackMatrixRef();
405 ::new (&mp_matrix) UmfpackMatrixRef(A);
406 }
407 }
408
409 // cached data to reduce reallocation, etc.
410 mutable LUMatrixType m_l;
411 int m_fact_errorCode;
412 UmfpackControl m_control;
413 mutable UmfpackInfo m_umfpackInfo;
414
415 mutable LUMatrixType m_u;
416 mutable IntColVectorType m_p;
417 mutable IntRowVectorType m_q;
418
419 UmfpackMatrixType m_dummy;
420 UmfpackMatrixRef mp_matrix;
421
422 void* m_numeric;
423 void* m_symbolic;
424
425 mutable ComputationInfo m_info;
426 int m_factorizationIsOk;
427 int m_analysisIsOk;
428 mutable bool m_extractedDataAreDirty;
429
430 private:
431 UmfPackLU(const UmfPackLU& ) { }
432};
433
434
435template<typename MatrixType>
436void UmfPackLU<MatrixType>::extractData() const
437{
438 if (m_extractedDataAreDirty)
439 {
440 // get size of the data
441 int lnz, unz, rows, cols, nz_udiag;
442 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
443
444 // allocate data
445 m_l.resize(rows,(std::min)(rows,cols));
446 m_l.resizeNonZeros(lnz);
447
448 m_u.resize((std::min)(rows,cols),cols);
449 m_u.resizeNonZeros(unz);
450
451 m_p.resize(rows);
452 m_q.resize(cols);
453
454 // extract
455 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
456 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
457 m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
458
459 m_extractedDataAreDirty = false;
460 }
461}
462
463template<typename MatrixType>
464typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const
465{
466 Scalar det;
467 umfpack_get_determinant(&det, 0, m_numeric, 0);
468 return det;
469}
470
471template<typename MatrixType>
472template<typename BDerived,typename XDerived>
473bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
474{
475 Index rhsCols = b.cols();
476 eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
477 eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
478 eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
479
480 int errorCode;
481 Scalar* x_ptr = 0;
482 Matrix<Scalar,Dynamic,1> x_tmp;
483 if(x.innerStride()!=1)
484 {
485 x_tmp.resize(x.rows());
486 x_ptr = x_tmp.data();
487 }
488 for (int j=0; j<rhsCols; ++j)
489 {
490 if(x.innerStride()==1)
491 x_ptr = &x.col(j).coeffRef(0);
492 errorCode = umfpack_solve(UMFPACK_A,
493 mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
494 x_ptr, &b.const_cast_derived().col(j).coeffRef(0), m_numeric, m_control.data(), m_umfpackInfo.data());
495 if(x.innerStride()!=1)
496 x.col(j) = x_tmp;
497 if (errorCode!=0)
498 return false;
499 }
500
501 return true;
502}
503
504} // end namespace Eigen
505
506#endif // EIGEN_UMFPACKSUPPORT_H
Eigen::Array
General-purpose arrays with easy API for coefficient-wise operations.
Definition Array.h:47
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:50
Eigen::Matrix
The matrix class, also used for vectors and row-vectors.
Definition Matrix.h:180
Eigen::PlainObjectBase::data
const Scalar * data() const
Definition PlainObjectBase.h:255
Eigen::Ref
A matrix or vector expression mapping an existing expression.
Definition Ref.h:195
Eigen::SparseMatrix
A versatible sparse matrix representation.
Definition SparseMatrix.h:98
Eigen::SparseSolverBase< UmfPackLU< _MatrixType > >::SparseSolverBase
SparseSolverBase()
Definition SparseSolverBase.h:72
Eigen::UmfPackLU::compute
void compute(const InputMatrixType &matrix)
Definition UmfPackSupport.h:242
Eigen::UmfPackLU::factorize
void factorize(const InputMatrixType &matrix)
Definition UmfPackSupport.h:308
Eigen::UmfPackLU::umfpackControl
const UmfpackControl & umfpackControl() const
Definition UmfPackSupport.h:285
Eigen::UmfPackLU::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition UmfPackSupport.h:207
Eigen::UmfPackLU::umfpackControl
UmfpackControl & umfpackControl()
Definition UmfPackSupport.h:296
Eigen::UmfPackLU::umfpackFactorizeReturncode
int umfpackFactorizeReturncode() const
Definition UmfPackSupport.h:273
Eigen::UmfPackLU::umfpackReportStatus
void umfpackReportStatus()
Definition UmfPackSupport.h:342
Eigen::UmfPackLU::umfpackReportControl
void umfpackReportControl()
Definition UmfPackSupport.h:323
Eigen::UmfPackLU::analyzePattern
void analyzePattern(const InputMatrixType &matrix)
Definition UmfPackSupport.h:258
Eigen::UmfPackLU::umfpackReportInfo
void umfpackReportInfo()
Definition UmfPackSupport.h:332
Eigen::ComputationInfo
ComputationInfo
Definition Constants.h:430
Eigen::NumericalIssue
@ NumericalIssue
Definition Constants.h:434
Eigen::InvalidInput
@ InvalidInput
Definition Constants.h:439
Eigen::Success
@ Success
Definition Constants.h:432
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition Constants.h:61
Eigen
Namespace containing all symbols from the Eigen library.
Definition A05_PortingFrom2To3.dox:1
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition Meta.h:65