Eigen  3.2.10
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 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 
110 namespace internal {
111  template<typename T> struct umfpack_helper_is_sparse_plain : false_type {};
112  template<typename Scalar, int Options, typename StorageIndex>
113  struct umfpack_helper_is_sparse_plain<SparseMatrix<Scalar,Options,StorageIndex> >
114  : true_type {};
115  template<typename Scalar, int Options, typename StorageIndex>
116  struct umfpack_helper_is_sparse_plain<MappedSparseMatrix<Scalar,Options,StorageIndex> >
117  : true_type {};
118 }
119 
133 template<typename _MatrixType>
134 class UmfPackLU : internal::noncopyable
135 {
136  public:
137  typedef _MatrixType MatrixType;
138  typedef typename MatrixType::Scalar Scalar;
139  typedef typename MatrixType::RealScalar RealScalar;
140  typedef typename MatrixType::Index Index;
146 
147  public:
148 
149  UmfPackLU() { init(); }
150 
151  template<typename InputMatrixType>
152  UmfPackLU(const InputMatrixType& matrix)
153  {
154  init();
155  compute(matrix);
156  }
157 
158  ~UmfPackLU()
159  {
160  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
161  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
162  }
163 
164  inline Index rows() const { return m_copyMatrix.rows(); }
165  inline Index cols() const { return m_copyMatrix.cols(); }
166 
173  {
174  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
175  return m_info;
176  }
177 
178  inline const LUMatrixType& matrixL() const
179  {
180  if (m_extractedDataAreDirty) extractData();
181  return m_l;
182  }
183 
184  inline const LUMatrixType& matrixU() const
185  {
186  if (m_extractedDataAreDirty) extractData();
187  return m_u;
188  }
189 
190  inline const IntColVectorType& permutationP() const
191  {
192  if (m_extractedDataAreDirty) extractData();
193  return m_p;
194  }
195 
196  inline const IntRowVectorType& permutationQ() const
197  {
198  if (m_extractedDataAreDirty) extractData();
199  return m_q;
200  }
201 
206  template<typename InputMatrixType>
207  void compute(const InputMatrixType& matrix)
208  {
209  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
210  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
211  grapInput(matrix.derived());
212  analyzePattern_impl();
213  factorize_impl();
214  }
215 
220  template<typename Rhs>
221  inline const internal::solve_retval<UmfPackLU, Rhs> solve(const MatrixBase<Rhs>& b) const
222  {
223  eigen_assert(m_isInitialized && "UmfPackLU is not initialized.");
224  eigen_assert(rows()==b.rows()
225  && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
226  return internal::solve_retval<UmfPackLU, Rhs>(*this, b.derived());
227  }
228 
233  template<typename Rhs>
234  inline const internal::sparse_solve_retval<UmfPackLU, Rhs> solve(const SparseMatrixBase<Rhs>& b) const
235  {
236  eigen_assert(m_isInitialized && "UmfPackLU is not initialized.");
237  eigen_assert(rows()==b.rows()
238  && "UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
239  return internal::sparse_solve_retval<UmfPackLU, Rhs>(*this, b.derived());
240  }
241 
248  template<typename InputMatrixType>
249  void analyzePattern(const InputMatrixType& matrix)
250  {
251  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
252  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
253 
254  grapInput(matrix.derived());
255 
256  analyzePattern_impl();
257  }
258 
265  template<typename InputMatrixType>
266  void factorize(const InputMatrixType& matrix)
267  {
268  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
269  if(m_numeric)
270  umfpack_free_numeric(&m_numeric,Scalar());
271 
272  grapInput(matrix.derived());
273 
274  factorize_impl();
275  }
276 
277  #ifndef EIGEN_PARSED_BY_DOXYGEN
278 
279  template<typename BDerived,typename XDerived>
280  bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
281  #endif
282 
283  Scalar determinant() const;
284 
285  void extractData() const;
286 
287  protected:
288 
289  void init()
290  {
291  m_info = InvalidInput;
292  m_isInitialized = false;
293  m_numeric = 0;
294  m_symbolic = 0;
295  m_outerIndexPtr = 0;
296  m_innerIndexPtr = 0;
297  m_valuePtr = 0;
298  m_extractedDataAreDirty = true;
299  }
300 
301  template<typename InputMatrixType>
302  void grapInput_impl(const InputMatrixType& mat, internal::true_type)
303  {
304  m_copyMatrix.resize(mat.rows(), mat.cols());
305  if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() )
306  {
307  // non supported input -> copy
308  m_copyMatrix = mat;
309  m_outerIndexPtr = m_copyMatrix.outerIndexPtr();
310  m_innerIndexPtr = m_copyMatrix.innerIndexPtr();
311  m_valuePtr = m_copyMatrix.valuePtr();
312  }
313  else
314  {
315  m_outerIndexPtr = mat.outerIndexPtr();
316  m_innerIndexPtr = mat.innerIndexPtr();
317  m_valuePtr = mat.valuePtr();
318  }
319  }
320 
321  template<typename InputMatrixType>
322  void grapInput_impl(const InputMatrixType& mat, internal::false_type)
323  {
324  m_copyMatrix = mat;
325  m_outerIndexPtr = m_copyMatrix.outerIndexPtr();
326  m_innerIndexPtr = m_copyMatrix.innerIndexPtr();
327  m_valuePtr = m_copyMatrix.valuePtr();
328  }
329 
330  template<typename InputMatrixType>
331  void grapInput(const InputMatrixType& mat)
332  {
333  grapInput_impl(mat, internal::umfpack_helper_is_sparse_plain<InputMatrixType>());
334  }
335 
336  void analyzePattern_impl()
337  {
338  int errorCode = 0;
339  errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
340  &m_symbolic, 0, 0);
341 
342  m_isInitialized = true;
343  m_info = errorCode ? InvalidInput : Success;
344  m_analysisIsOk = true;
345  m_factorizationIsOk = false;
346  m_extractedDataAreDirty = true;
347  }
348 
349  void factorize_impl()
350  {
351  int errorCode;
352  errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
353  m_symbolic, &m_numeric, 0, 0);
354 
355  m_info = errorCode ? NumericalIssue : Success;
356  m_factorizationIsOk = true;
357  m_extractedDataAreDirty = true;
358  }
359 
360  // cached data to reduce reallocation, etc.
361  mutable LUMatrixType m_l;
362  mutable LUMatrixType m_u;
363  mutable IntColVectorType m_p;
364  mutable IntRowVectorType m_q;
365 
366  UmfpackMatrixType m_copyMatrix;
367  const Scalar* m_valuePtr;
368  const int* m_outerIndexPtr;
369  const int* m_innerIndexPtr;
370  void* m_numeric;
371  void* m_symbolic;
372 
373  mutable ComputationInfo m_info;
374  bool m_isInitialized;
375  int m_factorizationIsOk;
376  int m_analysisIsOk;
377  mutable bool m_extractedDataAreDirty;
378 
379  private:
380  UmfPackLU(UmfPackLU& ) { }
381 };
382 
383 
384 template<typename MatrixType>
386 {
387  if (m_extractedDataAreDirty)
388  {
389  // get size of the data
390  int lnz, unz, rows, cols, nz_udiag;
391  umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
392 
393  // allocate data
394  m_l.resize(rows,(std::min)(rows,cols));
395  m_l.resizeNonZeros(lnz);
396 
397  m_u.resize((std::min)(rows,cols),cols);
398  m_u.resizeNonZeros(unz);
399 
400  m_p.resize(rows);
401  m_q.resize(cols);
402 
403  // extract
404  umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
405  m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
406  m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
407 
408  m_extractedDataAreDirty = false;
409  }
410 }
411 
412 template<typename MatrixType>
413 typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const
414 {
415  Scalar det;
416  umfpack_get_determinant(&det, 0, m_numeric, 0);
417  return det;
418 }
419 
420 template<typename MatrixType>
421 template<typename BDerived,typename XDerived>
423 {
424  const int rhsCols = b.cols();
425  eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
426  eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
427  eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
428 
429  int errorCode;
430  for (int j=0; j<rhsCols; ++j)
431  {
432  errorCode = umfpack_solve(UMFPACK_A,
433  m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
434  &x.col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0);
435  if (errorCode!=0)
436  return false;
437  }
438 
439  return true;
440 }
441 
442 
443 namespace internal {
444 
445 template<typename _MatrixType, typename Rhs>
446 struct solve_retval<UmfPackLU<_MatrixType>, Rhs>
447  : solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
448 {
449  typedef UmfPackLU<_MatrixType> Dec;
450  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
451 
452  template<typename Dest> void evalTo(Dest& dst) const
453  {
454  dec()._solve(rhs(),dst);
455  }
456 };
457 
458 template<typename _MatrixType, typename Rhs>
459 struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs>
460  : sparse_solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
461 {
462  typedef UmfPackLU<_MatrixType> Dec;
463  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
464 
465  template<typename Dest> void evalTo(Dest& dst) const
466  {
467  this->defaultEvalTo(dst);
468  }
469 };
470 
471 } // end namespace internal
472 
473 } // end namespace Eigen
474 
475 #endif // EIGEN_UMFPACKSUPPORT_H
A sparse LU factorization and solver based on UmfPack.
Definition: UmfPackSupport.h:134
void factorize(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:266
Definition: Constants.h:378
Definition: LDLT.h:16
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: UmfPackSupport.h:172
ColXpr col(Index i)
Definition: DenseBase.h:733
Index rows() const
Definition: SparseMatrixBase.h:160
void compute(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:207
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:239
Derived & derived()
Definition: EigenBase.h:34
Definition: Constants.h:383
const internal::solve_retval< UmfPackLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: UmfPackSupport.h:221
Definition: Eigen_Colamd.h:50
const internal::sparse_solve_retval< UmfPackLU, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
Definition: UmfPackSupport.h:234
Definition: Constants.h:376
const unsigned int RowMajorBit
Definition: Constants.h:53
void analyzePattern(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:249
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48