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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_UMFPACKSUPPORT_H
26 #define EIGEN_UMFPACKSUPPORT_H
27 
28 namespace Eigen {
29 
30 /* TODO extract L, extract U, compute det, etc... */
31 
32 // generic double/complex<double> wrapper functions:
33 
34 inline void umfpack_free_numeric(void **Numeric, double)
35 { umfpack_di_free_numeric(Numeric); *Numeric = 0; }
36 
37 inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
38 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
39 
40 inline void umfpack_free_symbolic(void **Symbolic, double)
41 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
42 
43 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
44 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
45 
46 inline int umfpack_symbolic(int n_row,int n_col,
47  const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
48  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
49 {
50  return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
51 }
52 
53 inline int umfpack_symbolic(int n_row,int n_col,
54  const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
55  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
56 {
57  return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&internal::real_ref(Ax[0]),0,Symbolic,Control,Info);
58 }
59 
60 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
61  void *Symbolic, void **Numeric,
62  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
63 {
64  return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
65 }
66 
67 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
68  void *Symbolic, void **Numeric,
69  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
70 {
71  return umfpack_zi_numeric(Ap,Ai,&internal::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
72 }
73 
74 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
75  double X[], const double B[], void *Numeric,
76  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
77 {
78  return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
79 }
80 
81 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
82  std::complex<double> X[], const std::complex<double> B[], void *Numeric,
83  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
84 {
85  return umfpack_zi_solve(sys,Ap,Ai,&internal::real_ref(Ax[0]),0,&internal::real_ref(X[0]),0,&internal::real_ref(B[0]),0,Numeric,Control,Info);
86 }
87 
88 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
89 {
90  return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
91 }
92 
93 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
94 {
95  return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
96 }
97 
98 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
99  int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
100 {
101  return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
102 }
103 
104 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
105  int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
106 {
107  double& lx0_real = internal::real_ref(Lx[0]);
108  double& ux0_real = internal::real_ref(Ux[0]);
109  double& dx0_real = internal::real_ref(Dx[0]);
110  return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
111  Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
112 }
113 
114 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
115 {
116  return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
117 }
118 
119 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
120 {
121  double& mx_real = internal::real_ref(*Mx);
122  return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
123 }
124 
138 template<typename _MatrixType>
139 class UmfPackLU : internal::noncopyable
140 {
141  public:
142  typedef _MatrixType MatrixType;
143  typedef typename MatrixType::Scalar Scalar;
144  typedef typename MatrixType::RealScalar RealScalar;
145  typedef typename MatrixType::Index Index;
151 
152  public:
153 
154  UmfPackLU() { init(); }
155 
156  UmfPackLU(const MatrixType& matrix)
157  {
158  init();
159  compute(matrix);
160  }
161 
163  {
166  }
167 
168  inline Index rows() const { return m_copyMatrix.rows(); }
169  inline Index cols() const { return m_copyMatrix.cols(); }
170 
177  {
178  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
179  return m_info;
180  }
181 
182  inline const LUMatrixType& matrixL() const
183  {
185  return m_l;
186  }
187 
188  inline const LUMatrixType& matrixU() const
189  {
191  return m_u;
192  }
193 
194  inline const IntColVectorType& permutationP() const
195  {
197  return m_p;
198  }
199 
200  inline const IntRowVectorType& permutationQ() const
201  {
203  return m_q;
204  }
205 
210  void compute(const MatrixType& matrix)
211  {
212  analyzePattern(matrix);
213  factorize(matrix);
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  void analyzePattern(const MatrixType& matrix)
249  {
250  if(m_symbolic)
252  if(m_numeric)
254 
255  grapInput(matrix);
256 
257  int errorCode = 0;
258  errorCode = umfpack_symbolic(matrix.rows(), matrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
259  &m_symbolic, 0, 0);
260 
261  m_isInitialized = true;
262  m_info = errorCode ? InvalidInput : Success;
263  m_analysisIsOk = true;
264  m_factorizationIsOk = false;
265  }
266 
273  void factorize(const MatrixType& matrix)
274  {
275  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
276  if(m_numeric)
278 
279  grapInput(matrix);
280 
281  int errorCode;
283  m_symbolic, &m_numeric, 0, 0);
284 
285  m_info = errorCode ? NumericalIssue : Success;
286  m_factorizationIsOk = true;
287  }
288 
289  #ifndef EIGEN_PARSED_BY_DOXYGEN
290 
291  template<typename BDerived,typename XDerived>
292  bool _solve(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
293  #endif
294 
295  Scalar determinant() const;
296 
297  void extractData() const;
298 
299  protected:
300 
301 
302  void init()
303  {
305  m_isInitialized = false;
306  m_numeric = 0;
307  m_symbolic = 0;
308  m_outerIndexPtr = 0;
309  m_innerIndexPtr = 0;
310  m_valuePtr = 0;
311  }
312 
313  void grapInput(const MatrixType& mat)
314  {
315  m_copyMatrix.resize(mat.rows(), mat.cols());
316  if( ((MatrixType::Flags&RowMajorBit)==RowMajorBit) || sizeof(typename MatrixType::Index)!=sizeof(int) || !mat.isCompressed() )
317  {
318  // non supported input -> copy
319  m_copyMatrix = mat;
323  }
324  else
325  {
326  m_outerIndexPtr = mat.outerIndexPtr();
327  m_innerIndexPtr = mat.innerIndexPtr();
328  m_valuePtr = mat.valuePtr();
329  }
330  }
331 
332  // cached data to reduce reallocation, etc.
333  mutable LUMatrixType m_l;
334  mutable LUMatrixType m_u;
337 
340  const int* m_outerIndexPtr;
341  const int* m_innerIndexPtr;
342  void* m_numeric;
343  void* m_symbolic;
344 
350 
351  private:
352  UmfPackLU(UmfPackLU& ) { }
353 };
354 
355 
356 template<typename MatrixType>
358 {
359  if (m_extractedDataAreDirty)
360  {
361  // get size of the data
362  int lnz, unz, rows, cols, nz_udiag;
363  umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
364 
365  // allocate data
366  m_l.resize(rows,(std::min)(rows,cols));
367  m_l.resizeNonZeros(lnz);
368 
369  m_u.resize((std::min)(rows,cols),cols);
370  m_u.resizeNonZeros(unz);
371 
372  m_p.resize(rows);
373  m_q.resize(cols);
374 
375  // extract
376  umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
377  m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
378  m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
379 
380  m_extractedDataAreDirty = false;
381  }
382 }
383 
384 template<typename MatrixType>
386 {
387  Scalar det;
388  umfpack_get_determinant(&det, 0, m_numeric, 0);
389  return det;
390 }
391 
392 template<typename MatrixType>
393 template<typename BDerived,typename XDerived>
395 {
396  const int rhsCols = b.cols();
397  eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
398  eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
399 
400  int errorCode;
401  for (int j=0; j<rhsCols; ++j)
402  {
403  errorCode = umfpack_solve(UMFPACK_A,
404  m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
405  &x.col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0);
406  if (errorCode!=0)
407  return false;
408  }
409 
410  return true;
411 }
412 
413 
414 namespace internal {
415 
416 template<typename _MatrixType, typename Rhs>
417 struct solve_retval<UmfPackLU<_MatrixType>, Rhs>
418  : solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
419 {
420  typedef UmfPackLU<_MatrixType> Dec;
421  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
422 
423  template<typename Dest> void evalTo(Dest& dst) const
424  {
425  dec()._solve(rhs(),dst);
426  }
427 };
428 
429 template<typename _MatrixType, typename Rhs>
430 struct sparse_solve_retval<UmfPackLU<_MatrixType>, Rhs>
431  : sparse_solve_retval_base<UmfPackLU<_MatrixType>, Rhs>
432 {
433  typedef UmfPackLU<_MatrixType> Dec;
435 
436  template<typename Dest> void evalTo(Dest& dst) const
437  {
438  dec()._solve(rhs(),dst);
439  }
440 };
441 
442 } // end namespace internal
443 
444 } // end namespace Eigen
445 
446 #endif // EIGEN_UMFPACKSUPPORT_H