UmfPackSupport.h
Go to the documentation of this file.
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 namespace Eigen {
14 
15 /* TODO extract L, extract U, compute det, etc... */
16 
17 // generic double/complex<double> wrapper functions:
18 
19 
20 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double)
21 { umfpack_di_defaults(control); }
22 
23 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>)
24 { umfpack_zi_defaults(control); }
25 
26 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double)
27 { umfpack_di_report_info(control, info);}
28 
29 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>)
30 { umfpack_zi_report_info(control, info);}
31 
32 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double)
33 { umfpack_di_report_status(control, status);}
34 
35 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>)
36 { umfpack_zi_report_status(control, status);}
37 
38 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double)
39 { umfpack_di_report_control(control);}
40 
41 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>)
42 { umfpack_zi_report_control(control);}
43 
44 inline void umfpack_free_numeric(void **Numeric, double)
45 { umfpack_di_free_numeric(Numeric); *Numeric = 0; }
46 
47 inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
48 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
49 
50 inline void umfpack_free_symbolic(void **Symbolic, double)
51 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
52 
53 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
54 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
55 
56 inline 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 
63 inline 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 
70 inline 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 
77 inline 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 
84 inline 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 
91 inline 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 
98 inline 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 
103 inline 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 
108 inline 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 
114 inline 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 
124 inline 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 
129 inline 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 
151 template<typename _MatrixType>
152 class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
153 {
154  protected:
156  using Base::m_isInitialized;
157  public:
158  using Base::_solve_impl;
159  typedef _MatrixType MatrixType;
160  typedef typename MatrixType::Scalar Scalar;
162  typedef typename MatrixType::StorageIndex StorageIndex;
169  enum {
170  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
171  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
172  };
173 
174  public:
175 
178 
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 
194  {
197  }
198 
199  inline Index rows() const { return mp_matrix.rows(); }
200  inline Index cols() const { return mp_matrix.cols(); }
201 
208  {
209  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
210  return m_info;
211  }
212 
213  inline const LUMatrixType& matrixL() const
214  {
216  return m_l;
217  }
218 
219  inline const LUMatrixType& matrixU() const
220  {
222  return m_u;
223  }
224 
225  inline const IntColVectorType& permutationP() const
226  {
228  return m_p;
229  }
230 
231  inline const IntRowVectorType& permutationQ() const
232  {
234  return m_q;
235  }
236 
241  template<typename InputMatrixType>
242  void compute(const InputMatrixType& matrix)
243  {
246  grab(matrix.derived());
248  factorize_impl();
249  }
250 
257  template<typename InputMatrixType>
258  void analyzePattern(const InputMatrixType& matrix)
259  {
262 
263  grab(matrix.derived());
264 
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)
313 
314  grab(matrix.derived());
315 
316  factorize_impl();
317  }
318 
324  {
326  }
327 
333  {
334  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
336  }
337 
343  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
345  }
346 
348  template<typename BDerived,typename XDerived>
350 
351  Scalar determinant() const;
352 
353  void extractData() const;
354 
355  protected:
356 
357  void init()
358  {
360  m_isInitialized = false;
361  m_numeric = 0;
362  m_symbolic = 0;
364 
366  }
367 
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(),
374 
375  m_isInitialized = true;
377  m_analysisIsOk = true;
378  m_factorizationIsOk = false;
380  }
381 
383  {
384 
385  m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
387 
388  m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
389  m_factorizationIsOk = true;
391  }
392 
393  template<typename MatrixDerived>
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;
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 
429 
430  private:
431  UmfPackLU(const UmfPackLU& ) { }
432 };
433 
434 
435 template<typename MatrixType>
437 {
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
457  m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
458 
459  m_extractedDataAreDirty = false;
460  }
461 }
462 
463 template<typename MatrixType>
465 {
466  Scalar det;
467  umfpack_get_determinant(&det, 0, m_numeric, 0);
468  return det;
469 }
470 
471 template<typename MatrixType>
472 template<typename BDerived,typename XDerived>
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;
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_DEVICE_FUNC ColXpr col(Index i)
Definition: DenseBase.h:839
void _solve_impl(const SparseMatrixBase< Rhs > &b, SparseMatrixBase< Dest > &dest) const
Matrix< int, MatrixType::RowsAtCompileTime, 1 > IntColVectorType
SCALAR Scalar
Definition: bench_gemm.cpp:33
int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
A sparse LU factorization and solver based on UmfPack.
EIGEN_DEVICE_FUNC internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar &x)
void factorize(const InputMatrixType &matrix)
Scalar * b
Definition: benchVecAdd.cpp:17
void umfpackReportStatus()
bool _solve_impl(const MatrixBase< BDerived > &b, MatrixBase< XDerived > &x) const
void extractData() const
void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double)
#define min(a, b)
Definition: datatypes.h:19
void compute(const InputMatrixType &matrix)
void analyzePattern(const InputMatrixType &matrix)
void umfpack_defaults(double control[UMFPACK_CONTROL], double)
IntRowVectorType m_q
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:621
A base class for sparse solvers.
int umfpack_numeric(const int Ap[], const int Ai[], const double Ax[], void *Symbolic, void **Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
Matrix< Scalar, Dynamic, 1 > Vector
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
const IntColVectorType & permutationP() const
int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info[UMFPACK_INFO])
LUMatrixType m_u
int umfpack_symbolic(int n_row, int n_col, const int Ap[], const int Ai[], const double Ax[], void **Symbolic, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
const Scalar * valuePtr() const
Definition: SparseMatrix.h:148
const unsigned int RowMajorBit
Definition: Constants.h:61
void resizeNonZeros(Index size)
Definition: SparseMatrix.h:644
MatrixType::RealScalar RealScalar
UmfpackMatrixType m_dummy
Index rows() const
const IntRowVectorType & permutationQ() const
SparseMatrix< Scalar > LUMatrixType
void umfpack_free_numeric(void **Numeric, double)
else if n * info
void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double)
IntColVectorType m_p
void grab(const UmfpackMatrixRef &A)
MatrixType::StorageIndex StorageIndex
Index cols() const
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
ComputationInfo m_info
UmfpackControl m_control
_MatrixType MatrixType
UmfpackInfo m_umfpackInfo
void umfpackReportControl()
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
#define eigen_assert(x)
Definition: Macros.h:579
UmfpackMatrixRef mp_matrix
const LUMatrixType & matrixL() const
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:34
void grab(const EigenBase< MatrixDerived > &A)
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:166
ComputationInfo info() const
Reports whether previous computation was successful.
const LUMatrixType & matrixU() const
Ref< const UmfpackMatrixType, StandardCompressedFormat > UmfpackMatrixRef
Scalar determinant() const
Array< double, UMFPACK_CONTROL, 1 > UmfpackControl
void analyzePattern_impl()
Matrix< int, 1, MatrixType::ColsAtCompileTime > IntRowVectorType
The quaternion class used to represent 3D orientations and rotations.
LUMatrixType m_l
int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[], int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
MatrixType::Scalar Scalar
const UmfpackControl & umfpackControl() const
UmfpackControl & umfpackControl()
int umfpackFactorizeReturncode() const
int umfpack_solve(int sys, const int Ap[], const int Ai[], const double Ax[], double X[], const double B[], void *Numeric, const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:157
Map< Matrix< T, Dynamic, Dynamic, ColMajor >, 0, OuterStride<> > matrix(T *data, int rows, int cols, int stride)
UmfPackLU(const UmfPackLU &)
#define X
Definition: icosphere.cpp:20
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Array< double, UMFPACK_INFO, 1 > UmfpackInfo
ComputationInfo
Definition: Constants.h:430
UmfPackLU(const InputMatrixType &matrix)
EIGEN_DEVICE_FUNC Derived & derived()
Definition: EigenBase.h:45
SparseSolverBase< UmfPackLU< _MatrixType > > Base
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
void umfpack_free_symbolic(void **Symbolic, double)
SparseMatrix< Scalar, ColMajor, int > UmfpackMatrixType
std::ptrdiff_t j
void umfpack_report_control(double control[UMFPACK_CONTROL], double)


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:51:21