Eigen  3.2.92
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 
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_free_numeric(void **Numeric, double)
27 { umfpack_di_free_numeric(Numeric); *Numeric = 0; }
28 
29 inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
30 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
31 
32 inline void umfpack_free_symbolic(void **Symbolic, double)
33 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
34 
35 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
36 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
37 
38 inline int umfpack_symbolic(int n_row,int n_col,
39  const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
40  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
41 {
42  return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
43 }
44 
45 inline int umfpack_symbolic(int n_row,int n_col,
46  const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
47  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
48 {
49  return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
50 }
51 
52 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
53  void *Symbolic, void **Numeric,
54  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
55 {
56  return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
57 }
58 
59 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
60  void *Symbolic, void **Numeric,
61  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
62 {
63  return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
64 }
65 
66 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
67  double X[], const double B[], void *Numeric,
68  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
69 {
70  return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
71 }
72 
73 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
74  std::complex<double> X[], const std::complex<double> B[], void *Numeric,
75  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
76 {
77  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);
78 }
79 
80 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
81 {
82  return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
83 }
84 
85 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
86 {
87  return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
88 }
89 
90 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
91  int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
92 {
93  return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
94 }
95 
96 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
97  int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
98 {
99  double& lx0_real = numext::real_ref(Lx[0]);
100  double& ux0_real = numext::real_ref(Ux[0]);
101  double& dx0_real = numext::real_ref(Dx[0]);
102  return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
103  Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
104 }
105 
106 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
107 {
108  return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
109 }
110 
111 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
112 {
113  double& mx_real = numext::real_ref(*Mx);
114  return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
115 }
116 
117 
131 template<typename _MatrixType>
132 class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
133 {
134  protected:
136  using Base::m_isInitialized;
137  public:
138  using Base::_solve_impl;
139  typedef _MatrixType MatrixType;
140  typedef typename MatrixType::Scalar Scalar;
141  typedef typename MatrixType::RealScalar RealScalar;
142  typedef typename MatrixType::StorageIndex StorageIndex;
149  enum {
150  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
151  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
152  };
153 
154  public:
155 
157 
158  UmfPackLU()
159  : m_dummy(0,0), mp_matrix(m_dummy)
160  {
161  init();
162  }
163 
164  template<typename InputMatrixType>
165  explicit UmfPackLU(const InputMatrixType& matrix)
166  : mp_matrix(matrix)
167  {
168  init();
169  compute(matrix);
170  }
171 
172  ~UmfPackLU()
173  {
174  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
175  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
176  }
177 
178  inline Index rows() const { return mp_matrix.rows(); }
179  inline Index cols() const { return mp_matrix.cols(); }
180 
187  {
188  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
189  return m_info;
190  }
191 
192  inline const LUMatrixType& matrixL() const
193  {
194  if (m_extractedDataAreDirty) extractData();
195  return m_l;
196  }
197 
198  inline const LUMatrixType& matrixU() const
199  {
200  if (m_extractedDataAreDirty) extractData();
201  return m_u;
202  }
203 
204  inline const IntColVectorType& permutationP() const
205  {
206  if (m_extractedDataAreDirty) extractData();
207  return m_p;
208  }
209 
210  inline const IntRowVectorType& permutationQ() const
211  {
212  if (m_extractedDataAreDirty) extractData();
213  return m_q;
214  }
215 
220  template<typename InputMatrixType>
221  void compute(const InputMatrixType& matrix)
222  {
223  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
224  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
225  grab(matrix.derived());
226  analyzePattern_impl();
227  factorize_impl();
228  }
229 
236  template<typename InputMatrixType>
237  void analyzePattern(const InputMatrixType& matrix)
238  {
239  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
240  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
241 
242  grab(matrix.derived());
243 
244  analyzePattern_impl();
245  }
246 
252  inline int umfpackFactorizeReturncode() const
253  {
254  eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()");
255  return m_fact_errorCode;
256  }
257 
264  inline const UmfpackControl& umfpackControl() const
265  {
266  return m_control;
267  }
268 
275  inline UmfpackControl& umfpackControl()
276  {
277  return m_control;
278  }
279 
286  template<typename InputMatrixType>
287  void factorize(const InputMatrixType& matrix)
288  {
289  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
290  if(m_numeric)
291  umfpack_free_numeric(&m_numeric,Scalar());
292 
293  grab(matrix.derived());
294 
295  factorize_impl();
296  }
297 
299  template<typename BDerived,typename XDerived>
300  bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
301 
302  Scalar determinant() const;
303 
304  void extractData() const;
305 
306  protected:
307 
308  void init()
309  {
310  m_info = InvalidInput;
311  m_isInitialized = false;
312  m_numeric = 0;
313  m_symbolic = 0;
314  m_extractedDataAreDirty = true;
315  }
316 
317  void analyzePattern_impl()
318  {
319  umfpack_defaults(m_control.data(), Scalar());
320  int errorCode = 0;
321  errorCode = umfpack_symbolic(internal::convert_index<int>(mp_matrix.rows()),
322  internal::convert_index<int>(mp_matrix.cols()),
323  mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
324  &m_symbolic, m_control.data(), 0);
325 
326  m_isInitialized = true;
327  m_info = errorCode ? InvalidInput : Success;
328  m_analysisIsOk = true;
329  m_factorizationIsOk = false;
330  m_extractedDataAreDirty = true;
331  }
332 
333  void factorize_impl()
334  {
335  m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
336  m_symbolic, &m_numeric, m_control.data(), 0);
337 
338  m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
339  m_factorizationIsOk = true;
340  m_extractedDataAreDirty = true;
341  }
342 
343  template<typename MatrixDerived>
344  void grab(const EigenBase<MatrixDerived> &A)
345  {
346  mp_matrix.~UmfpackMatrixRef();
347  ::new (&mp_matrix) UmfpackMatrixRef(A.derived());
348  }
349 
350  void grab(const UmfpackMatrixRef &A)
351  {
352  if(&(A.derived()) != &mp_matrix)
353  {
354  mp_matrix.~UmfpackMatrixRef();
355  ::new (&mp_matrix) UmfpackMatrixRef(A);
356  }
357  }
358 
359  // cached data to reduce reallocation, etc.
360  mutable LUMatrixType m_l;
361  int m_fact_errorCode;
362  UmfpackControl m_control;
363 
364  mutable LUMatrixType m_u;
365  mutable IntColVectorType m_p;
366  mutable IntRowVectorType m_q;
367 
368  UmfpackMatrixType m_dummy;
369  UmfpackMatrixRef mp_matrix;
370 
371  void* m_numeric;
372  void* m_symbolic;
373 
374  mutable ComputationInfo m_info;
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>
385 void UmfPackLU<MatrixType>::extractData() const
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>
422 bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
423 {
424  Index 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  Scalar* x_ptr = 0;
431  Matrix<Scalar,Dynamic,1> x_tmp;
432  if(x.innerStride()!=1)
433  {
434  x_tmp.resize(x.rows());
435  x_ptr = x_tmp.data();
436  }
437  for (int j=0; j<rhsCols; ++j)
438  {
439  if(x.innerStride()==1)
440  x_ptr = &x.col(j).coeffRef(0);
441  errorCode = umfpack_solve(UMFPACK_A,
442  mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
443  x_ptr, &b.const_cast_derived().col(j).coeffRef(0), m_numeric, m_control.data(), 0);
444  if(x.innerStride()!=1)
445  x.col(j) = x_tmp;
446  if (errorCode!=0)
447  return false;
448  }
449 
450  return true;
451 }
452 
453 } // end namespace Eigen
454 
455 #endif // EIGEN_UMFPACKSUPPORT_H
A sparse LU factorization and solver based on UmfPack.
Definition: UmfPackSupport.h:132
void factorize(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:287
void compute(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:221
void analyzePattern(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:237
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Definition: LDLT.h:16
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:161
const unsigned int RowMajorBit
Definition: Constants.h:61
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
Definition: Constants.h:434
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:611
Definition: Constants.h:439
Definition: Constants.h:432
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: UmfPackSupport.h:186
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:152
const UmfpackControl & umfpackControl() const
Definition: UmfPackSupport.h:264
UmfpackControl & umfpackControl()
Definition: UmfPackSupport.h:275
int umfpackFactorizeReturncode() const
Definition: UmfPackSupport.h:252
const Scalar * data() const
Definition: PlainObjectBase.h:228
const Scalar * valuePtr() const
Definition: SparseMatrix.h:143
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48