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 inline void umfpack_free_numeric(void **Numeric, double)
20 { umfpack_di_free_numeric(Numeric); *Numeric = 0; }
21 
22 inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
23 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
24 
25 inline void umfpack_free_symbolic(void **Symbolic, double)
26 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
27 
28 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
29 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
30 
31 inline int umfpack_symbolic(int n_row,int n_col,
32  const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
33  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
34 {
35  return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
36 }
37 
38 inline int umfpack_symbolic(int n_row,int n_col,
39  const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
40  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
41 {
42  return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
43 }
44 
45 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
46  void *Symbolic, void **Numeric,
47  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
48 {
49  return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
50 }
51 
52 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
53  void *Symbolic, void **Numeric,
54  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
55 {
56  return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
57 }
58 
59 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
60  double X[], const double B[], void *Numeric,
61  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
62 {
63  return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
64 }
65 
66 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
67  std::complex<double> X[], const std::complex<double> B[], void *Numeric,
68  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
69 {
70  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);
71 }
72 
73 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
74 {
75  return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
76 }
77 
78 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
79 {
80  return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
81 }
82 
83 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
84  int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
85 {
86  return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
87 }
88 
89 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
90  int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
91 {
92  double& lx0_real = numext::real_ref(Lx[0]);
93  double& ux0_real = numext::real_ref(Ux[0]);
94  double& dx0_real = numext::real_ref(Dx[0]);
95  return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
96  Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
97 }
98 
99 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
100 {
101  return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
102 }
103 
104 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
105 {
106  double& mx_real = numext::real_ref(*Mx);
107  return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
108 }
109 
123 template<typename _MatrixType>
125 {
126  public:
127  typedef _MatrixType MatrixType;
128  typedef typename MatrixType::Scalar Scalar;
129  typedef typename MatrixType::RealScalar RealScalar;
130  typedef typename MatrixType::Index Index;
136 
137  public:
138 
139  UmfPackLU() { init(); }
140 
141  UmfPackLU(const MatrixType& matrix)
142  {
143  init();
144  compute(matrix);
145  }
146 
148  {
151  }
152 
153  inline Index rows() const { return m_copyMatrix.rows(); }
154  inline Index cols() const { return m_copyMatrix.cols(); }
155 
162  {
163  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
164  return m_info;
165  }
166 
167  inline const LUMatrixType& matrixL() const
168  {
170  return m_l;
171  }
172 
173  inline const LUMatrixType& matrixU() const
174  {
176  return m_u;
177  }
178 
179  inline const IntColVectorType& permutationP() const
180  {
182  return m_p;
183  }
184 
185  inline const IntRowVectorType& permutationQ() const
186  {
188  return m_q;
189  }
190 
195  void compute(const MatrixType& matrix)
196  {
197  analyzePattern(matrix);
198  factorize(matrix);
199  }
200 
205  template<typename Rhs>
207  {
208  eigen_assert(m_isInitialized && "UmfPackLU is not initialized.");
209  eigen_assert(rows()==b.rows()
210  && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
211  return internal::solve_retval<UmfPackLU, Rhs>(*this, b.derived());
212  }
213 
218  template<typename Rhs>
220  {
221  eigen_assert(m_isInitialized && "UmfPackLU is not initialized.");
222  eigen_assert(rows()==b.rows()
223  && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
225  }
226 
233  void analyzePattern(const MatrixType& matrix)
234  {
235  if(m_symbolic)
237  if(m_numeric)
239 
240  grapInput(matrix);
241 
242  int errorCode = 0;
243  errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
244  &m_symbolic, 0, 0);
245 
246  m_isInitialized = true;
247  m_info = errorCode ? InvalidInput : Success;
248  m_analysisIsOk = true;
249  m_factorizationIsOk = false;
250  }
251 
258  void factorize(const MatrixType& matrix)
259  {
260  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
261  if(m_numeric)
263 
264  grapInput(matrix);
265 
266  int errorCode;
268  m_symbolic, &m_numeric, 0, 0);
269 
270  m_info = errorCode ? NumericalIssue : Success;
271  m_factorizationIsOk = true;
272  }
273 
274  #ifndef EIGEN_PARSED_BY_DOXYGEN
275 
276  template<typename BDerived,typename XDerived>
277  bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
278  #endif
279 
280  Scalar determinant() const;
281 
282  void extractData() const;
283 
284  protected:
285 
286 
287  void init()
288  {
290  m_isInitialized = false;
291  m_numeric = 0;
292  m_symbolic = 0;
293  m_outerIndexPtr = 0;
294  m_innerIndexPtr = 0;
295  m_valuePtr = 0;
296  }
297 
298  void grapInput(const MatrixType& mat)
299  {
300  m_copyMatrix.resize(mat.rows(), mat.cols());
301  if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() )
302  {
303  // non supported input -> copy
304  m_copyMatrix = mat;
308  }
309  else
310  {
311  m_outerIndexPtr = mat.outerIndexPtr();
312  m_innerIndexPtr = mat.innerIndexPtr();
313  m_valuePtr = mat.valuePtr();
314  }
315  }
316 
317  // cached data to reduce reallocation, etc.
318  mutable LUMatrixType m_l;
319  mutable LUMatrixType m_u;
320  mutable IntColVectorType m_p;
321  mutable IntRowVectorType m_q;
322 
323  UmfpackMatrixType m_copyMatrix;
324  const Scalar* m_valuePtr;
325  const int* m_outerIndexPtr;
326  const int* m_innerIndexPtr;
327  void* m_numeric;
328  void* m_symbolic;
329 
335 
336  private:
338 };
339 
340 
341 template<typename MatrixType>
343 {
345  {
346  // get size of the data
347  int lnz, unz, rows, cols, nz_udiag;
348  umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
349 
350  // allocate data
351  m_l.resize(rows,(std::min)(rows,cols));
352  m_l.resizeNonZeros(lnz);
353 
354  m_u.resize((std::min)(rows,cols),cols);
355  m_u.resizeNonZeros(unz);
356 
357  m_p.resize(rows);
358  m_q.resize(cols);
359 
360  // extract
363  m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
364 
365  m_extractedDataAreDirty = false;
366  }
367 }
368 
369 template<typename MatrixType>
371 {
372  Scalar det;
373  umfpack_get_determinant(&det, 0, m_numeric, 0);
374  return det;
375 }
376 
377 template<typename MatrixType>
378 template<typename BDerived,typename XDerived>
380 {
381  const int rhsCols = b.cols();
382  eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
383  eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
384  eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
385 
386  int errorCode;
387  for (int j=0; j<rhsCols; ++j)
388  {
389  errorCode = umfpack_solve(UMFPACK_A,
391  &x.col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0);
392  if (errorCode!=0)
393  return false;
394  }
395 
396  return true;
397 }
398 
399 
400 namespace internal {
401 
402 template<typename _MatrixType, typename Rhs>
403 struct solve_retval<UmfPackLU<_MatrixType>, Rhs>
404  : solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
405 {
407  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
408 
409  template<typename Dest> void evalTo(Dest& dst) const
410  {
411  dec()._solve(rhs(),dst);
412  }
413 };
414 
415 template<typename _MatrixType, typename Rhs>
416 struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs>
417  : sparse_solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
418 {
421 
422  template<typename Dest> void evalTo(Dest& dst) const
423  {
424  this->defaultEvalTo(dst);
425  }
426 };
427 
428 } // end namespace internal
429 
430 } // end namespace Eigen
431 
432 #endif // EIGEN_UMFPACKSUPPORT_H
Matrix< int, MatrixType::RowsAtCompileTime, 1 > IntColVectorType
MatrixType::Index Index
int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
ColXpr col(Index i)
Definition: DenseBase.h:709
Index cols() const
Definition: SparseMatrix.h:121
A sparse LU factorization and solver based on UmfPack.
void compute(const MatrixType &matrix)
void extractData() const
IntRowVectorType m_q
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])
void resizeNonZeros(Index size)
Definition: SparseMatrix.h:619
const internal::solve_retval< UmfPackLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Matrix< Scalar, Dynamic, 1 > Vector
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
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 unsigned int RowMajorBit
MatrixType::RealScalar RealScalar
const Index * outerIndexPtr() const
Definition: SparseMatrix.h:149
Index rows() const
const IntRowVectorType & permutationQ() const
SparseMatrix< Scalar > LUMatrixType
void umfpack_free_numeric(void **Numeric, double)
IntColVectorType m_p
Index cols() const
Base class of any sparse matrices or sparse expressions.
ComputationInfo m_info
_MatrixType MatrixType
void grapInput(const MatrixType &mat)
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
UmfPackLU(UmfPackLU &)
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:596
const LUMatrixType & matrixL() const
const Scalar * m_valuePtr
const int * m_outerIndexPtr
#define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: SparseSolve.h:71
ComputationInfo info() const
Reports whether previous computation was successful.
const LUMatrixType & matrixU() const
void rhs(const real_t *x, real_t *f)
EIGEN_STRONG_INLINE const Scalar * data() const
Scalar determinant() const
UmfpackMatrixType m_copyMatrix
Matrix< int, 1, MatrixType::ColsAtCompileTime > IntRowVectorType
const Derived & derived() const
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
void factorize(const MatrixType &matrix)
const int * m_innerIndexPtr
const Index * innerIndexPtr() const
Definition: SparseMatrix.h:140
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
internal::add_const_on_value_type< EIGEN_MATHFUNC_RETVAL(real_ref, Scalar) >::type real_ref(const Scalar &x)
const internal::sparse_solve_retval< UmfPackLU, Rhs > solve(const SparseMatrixBase< Rhs > &b) 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])
#define eigen_assert(x)
const Scalar * valuePtr() const
Definition: SparseMatrix.h:131
bool _solve(const MatrixBase< BDerived > &b, MatrixBase< XDerived > &x) const
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Index rows() const
Definition: SparseMatrix.h:119
void umfpack_free_symbolic(void **Symbolic, double)
SparseMatrix< Scalar, ColMajor, int > UmfpackMatrixType
UmfPackLU(const MatrixType &matrix)
void analyzePattern(const MatrixType &matrix)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:15