Path Tracer
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 // for compatibility with super old version of umfpack,
14 // not sure this is really needed, but this is harmless.
15 #ifndef SuiteSparse_long
16 #ifdef UF_long
17 #define SuiteSparse_long UF_long
18 #else
19 #error neither SuiteSparse_long nor UF_long are defined
20 #endif
21 #endif
22 
23 namespace Eigen {
24 
25 /* TODO extract L, extract U, compute det, etc... */
26 
27 // generic double/complex<double> wrapper functions:
28 
29 
30  // Defaults
31 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, int)
32 { umfpack_di_defaults(control); }
33 
34 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, int)
35 { umfpack_zi_defaults(control); }
36 
37 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, SuiteSparse_long)
38 { umfpack_dl_defaults(control); }
39 
40 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long)
41 { umfpack_zl_defaults(control); }
42 
43 // Report info
44 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, int)
45 { umfpack_di_report_info(control, info);}
46 
47 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, int)
48 { umfpack_zi_report_info(control, info);}
49 
50 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, SuiteSparse_long)
51 { umfpack_dl_report_info(control, info);}
52 
53 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, SuiteSparse_long)
54 { umfpack_zl_report_info(control, info);}
55 
56 // Report status
57 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, int)
58 { umfpack_di_report_status(control, status);}
59 
60 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, int)
61 { umfpack_zi_report_status(control, status);}
62 
63 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, SuiteSparse_long)
64 { umfpack_dl_report_status(control, status);}
65 
66 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, SuiteSparse_long)
67 { umfpack_zl_report_status(control, status);}
68 
69 // report control
70 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, int)
71 { umfpack_di_report_control(control);}
72 
73 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, int)
74 { umfpack_zi_report_control(control);}
75 
76 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, SuiteSparse_long)
77 { umfpack_dl_report_control(control);}
78 
79 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long)
80 { umfpack_zl_report_control(control);}
81 
82 // Free numeric
83 inline void umfpack_free_numeric(void **Numeric, double, int)
84 { umfpack_di_free_numeric(Numeric); *Numeric = 0; }
85 
86 inline void umfpack_free_numeric(void **Numeric, std::complex<double>, int)
87 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
88 
89 inline void umfpack_free_numeric(void **Numeric, double, SuiteSparse_long)
90 { umfpack_dl_free_numeric(Numeric); *Numeric = 0; }
91 
92 inline void umfpack_free_numeric(void **Numeric, std::complex<double>, SuiteSparse_long)
93 { umfpack_zl_free_numeric(Numeric); *Numeric = 0; }
94 
95 // Free symbolic
96 inline void umfpack_free_symbolic(void **Symbolic, double, int)
97 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
98 
99 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, int)
100 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
101 
102 inline void umfpack_free_symbolic(void **Symbolic, double, SuiteSparse_long)
103 { umfpack_dl_free_symbolic(Symbolic); *Symbolic = 0; }
104 
105 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, SuiteSparse_long)
106 { umfpack_zl_free_symbolic(Symbolic); *Symbolic = 0; }
107 
108 // Symbolic
109 inline int umfpack_symbolic(int n_row,int n_col,
110  const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
111  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
112 {
113  return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
114 }
115 
116 inline int umfpack_symbolic(int n_row,int n_col,
117  const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
118  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
119 {
120  return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
121 }
122 inline SuiteSparse_long umfpack_symbolic( SuiteSparse_long n_row,SuiteSparse_long n_col,
123  const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[], void **Symbolic,
124  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
125 {
126  return umfpack_dl_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
127 }
128 
129 inline SuiteSparse_long umfpack_symbolic( SuiteSparse_long n_row,SuiteSparse_long n_col,
130  const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[], void **Symbolic,
131  const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
132 {
133  return umfpack_zl_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
134 }
135 
136 // Numeric
137 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
138  void *Symbolic, void **Numeric,
139  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
140 {
141  return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
142 }
143 
144 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
145  void *Symbolic, void **Numeric,
146  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
147 {
148  return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
149 }
150 inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[],
151  void *Symbolic, void **Numeric,
152  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
153 {
154  return umfpack_dl_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
155 }
156 
157 inline SuiteSparse_long umfpack_numeric(const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[],
158  void *Symbolic, void **Numeric,
159  const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
160 {
161  return umfpack_zl_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
162 }
163 
164 // solve
165 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
166  double X[], const double B[], void *Numeric,
167  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
168 {
169  return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
170 }
171 
172 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
173  std::complex<double> X[], const std::complex<double> B[], void *Numeric,
174  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
175 {
176  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);
177 }
178 
179 inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const double Ax[],
180  double X[], const double B[], void *Numeric,
181  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
182 {
183  return umfpack_dl_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
184 }
185 
186 inline SuiteSparse_long umfpack_solve(int sys, const SuiteSparse_long Ap[], const SuiteSparse_long Ai[], const std::complex<double> Ax[],
187  std::complex<double> X[], const std::complex<double> B[], void *Numeric,
188  const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
189 {
190  return umfpack_zl_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);
191 }
192 
193 // Get Lunz
194 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
195 {
196  return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
197 }
198 
199 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
200 {
201  return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
202 }
203 
204 inline SuiteSparse_long umfpack_get_lunz( SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row, SuiteSparse_long *n_col,
205  SuiteSparse_long *nz_udiag, void *Numeric, double)
206 {
207  return umfpack_dl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
208 }
209 
210 inline SuiteSparse_long umfpack_get_lunz( SuiteSparse_long *lnz, SuiteSparse_long *unz, SuiteSparse_long *n_row, SuiteSparse_long *n_col,
211  SuiteSparse_long *nz_udiag, void *Numeric, std::complex<double>)
212 {
213  return umfpack_zl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
214 }
215 
216 // Get Numeric
217 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
218  int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
219 {
220  return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
221 }
222 
223 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
224  int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
225 {
226  double& lx0_real = numext::real_ref(Lx[0]);
227  double& ux0_real = numext::real_ref(Ux[0]);
228  double& dx0_real = numext::real_ref(Dx[0]);
229  return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
230  Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
231 }
232 inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], double Lx[], SuiteSparse_long Up[], SuiteSparse_long Ui[], double Ux[],
233  SuiteSparse_long P[], SuiteSparse_long Q[], double Dx[], SuiteSparse_long *do_recip, double Rs[], void *Numeric)
234 {
235  return umfpack_dl_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
236 }
237 
238 inline SuiteSparse_long umfpack_get_numeric(SuiteSparse_long Lp[], SuiteSparse_long Lj[], std::complex<double> Lx[], SuiteSparse_long Up[], SuiteSparse_long Ui[], std::complex<double> Ux[],
239  SuiteSparse_long P[], SuiteSparse_long Q[], std::complex<double> Dx[], SuiteSparse_long *do_recip, double Rs[], void *Numeric)
240 {
241  double& lx0_real = numext::real_ref(Lx[0]);
242  double& ux0_real = numext::real_ref(Ux[0]);
243  double& dx0_real = numext::real_ref(Dx[0]);
244  return umfpack_zl_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
245  Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
246 }
247 
248 // Get Determinant
249 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int)
250 {
251  return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
252 }
253 
254 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int)
255 {
256  double& mx_real = numext::real_ref(*Mx);
257  return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
258 }
259 
260 inline SuiteSparse_long umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], SuiteSparse_long)
261 {
262  return umfpack_dl_get_determinant(Mx,Ex,NumericHandle,User_Info);
263 }
264 
265 inline SuiteSparse_long umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], SuiteSparse_long)
266 {
267  double& mx_real = numext::real_ref(*Mx);
268  return umfpack_zl_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
269 }
270 
271 
287 template<typename _MatrixType>
288 class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
289 {
290  protected:
291  typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base;
292  using Base::m_isInitialized;
293  public:
294  using Base::_solve_impl;
295  typedef _MatrixType MatrixType;
296  typedef typename MatrixType::Scalar Scalar;
297  typedef typename MatrixType::RealScalar RealScalar;
298  typedef typename MatrixType::StorageIndex StorageIndex;
299  typedef Matrix<Scalar,Dynamic,1> Vector;
300  typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
301  typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
302  typedef SparseMatrix<Scalar> LUMatrixType;
303  typedef SparseMatrix<Scalar,ColMajor,StorageIndex> UmfpackMatrixType;
304  typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef;
305  enum {
306  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
307  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
308  };
309 
310  public:
311 
312  typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl;
313  typedef Array<double, UMFPACK_INFO, 1> UmfpackInfo;
314 
315  UmfPackLU()
316  : m_dummy(0,0), mp_matrix(m_dummy)
317  {
318  init();
319  }
320 
321  template<typename InputMatrixType>
322  explicit UmfPackLU(const InputMatrixType& matrix)
323  : mp_matrix(matrix)
324  {
325  init();
326  compute(matrix);
327  }
328 
329  ~UmfPackLU()
330  {
331  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(), StorageIndex());
332  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(), StorageIndex());
333  }
334 
335  inline Index rows() const { return mp_matrix.rows(); }
336  inline Index cols() const { return mp_matrix.cols(); }
337 
343  ComputationInfo info() const
344  {
345  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
346  return m_info;
347  }
348 
349  inline const LUMatrixType& matrixL() const
350  {
351  if (m_extractedDataAreDirty) extractData();
352  return m_l;
353  }
354 
355  inline const LUMatrixType& matrixU() const
356  {
357  if (m_extractedDataAreDirty) extractData();
358  return m_u;
359  }
360 
361  inline const IntColVectorType& permutationP() const
362  {
363  if (m_extractedDataAreDirty) extractData();
364  return m_p;
365  }
366 
367  inline const IntRowVectorType& permutationQ() const
368  {
369  if (m_extractedDataAreDirty) extractData();
370  return m_q;
371  }
372 
377  template<typename InputMatrixType>
378  void compute(const InputMatrixType& matrix)
379  {
380  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
381  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
382  grab(matrix.derived());
383  analyzePattern_impl();
384  factorize_impl();
385  }
386 
393  template<typename InputMatrixType>
394  void analyzePattern(const InputMatrixType& matrix)
395  {
396  if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
397  if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
398 
399  grab(matrix.derived());
400 
401  analyzePattern_impl();
402  }
403 
409  inline int umfpackFactorizeReturncode() const
410  {
411  eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()");
412  return m_fact_errorCode;
413  }
414 
421  inline const UmfpackControl& umfpackControl() const
422  {
423  return m_control;
424  }
425 
432  inline UmfpackControl& umfpackControl()
433  {
434  return m_control;
435  }
436 
443  template<typename InputMatrixType>
444  void factorize(const InputMatrixType& matrix)
445  {
446  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
447  if(m_numeric)
448  umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
449 
450  grab(matrix.derived());
451 
452  factorize_impl();
453  }
454 
459  void printUmfpackControl()
460  {
461  umfpack_report_control(m_control.data(), Scalar(),StorageIndex());
462  }
463 
468  void printUmfpackInfo()
469  {
470  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
471  umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar(),StorageIndex());
472  }
473 
478  void printUmfpackStatus() {
479  eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
480  umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar(),StorageIndex());
481  }
482 
484  template<typename BDerived,typename XDerived>
485  bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
486 
487  Scalar determinant() const;
488 
489  void extractData() const;
490 
491  protected:
492 
493  void init()
494  {
495  m_info = InvalidInput;
496  m_isInitialized = false;
497  m_numeric = 0;
498  m_symbolic = 0;
499  m_extractedDataAreDirty = true;
500 
501  umfpack_defaults(m_control.data(), Scalar(),StorageIndex());
502  }
503 
504  void analyzePattern_impl()
505  {
506  m_fact_errorCode = umfpack_symbolic(internal::convert_index<StorageIndex>(mp_matrix.rows()),
507  internal::convert_index<StorageIndex>(mp_matrix.cols()),
508  mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
509  &m_symbolic, m_control.data(), m_umfpackInfo.data());
510 
511  m_isInitialized = true;
512  m_info = m_fact_errorCode ? InvalidInput : Success;
513  m_analysisIsOk = true;
514  m_factorizationIsOk = false;
515  m_extractedDataAreDirty = true;
516  }
517 
518  void factorize_impl()
519  {
520 
521  m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
522  m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data());
523 
524  m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
525  m_factorizationIsOk = true;
526  m_extractedDataAreDirty = true;
527  }
528 
529  template<typename MatrixDerived>
530  void grab(const EigenBase<MatrixDerived> &A)
531  {
532  mp_matrix.~UmfpackMatrixRef();
533  ::new (&mp_matrix) UmfpackMatrixRef(A.derived());
534  }
535 
536  void grab(const UmfpackMatrixRef &A)
537  {
538  if(&(A.derived()) != &mp_matrix)
539  {
540  mp_matrix.~UmfpackMatrixRef();
541  ::new (&mp_matrix) UmfpackMatrixRef(A);
542  }
543  }
544 
545  // cached data to reduce reallocation, etc.
546  mutable LUMatrixType m_l;
547  StorageIndex m_fact_errorCode;
548  UmfpackControl m_control;
549  mutable UmfpackInfo m_umfpackInfo;
550 
551  mutable LUMatrixType m_u;
552  mutable IntColVectorType m_p;
553  mutable IntRowVectorType m_q;
554 
555  UmfpackMatrixType m_dummy;
556  UmfpackMatrixRef mp_matrix;
557 
558  void* m_numeric;
559  void* m_symbolic;
560 
561  mutable ComputationInfo m_info;
562  int m_factorizationIsOk;
563  int m_analysisIsOk;
564  mutable bool m_extractedDataAreDirty;
565 
566  private:
567  UmfPackLU(const UmfPackLU& ) { }
568 };
569 
570 
571 template<typename MatrixType>
572 void UmfPackLU<MatrixType>::extractData() const
573 {
574  if (m_extractedDataAreDirty)
575  {
576  // get size of the data
577  StorageIndex lnz, unz, rows, cols, nz_udiag;
578  umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
579 
580  // allocate data
581  m_l.resize(rows,(std::min)(rows,cols));
582  m_l.resizeNonZeros(lnz);
583 
584  m_u.resize((std::min)(rows,cols),cols);
585  m_u.resizeNonZeros(unz);
586 
587  m_p.resize(rows);
588  m_q.resize(cols);
589 
590  // extract
591  umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
592  m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
593  m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
594 
595  m_extractedDataAreDirty = false;
596  }
597 }
598 
599 template<typename MatrixType>
600 typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const
601 {
602  Scalar det;
603  umfpack_get_determinant(&det, 0, m_numeric, 0, StorageIndex());
604  return det;
605 }
606 
607 template<typename MatrixType>
608 template<typename BDerived,typename XDerived>
609 bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
610 {
611  Index rhsCols = b.cols();
612  eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
613  eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
614  eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
615 
616  Scalar* x_ptr = 0;
617  Matrix<Scalar,Dynamic,1> x_tmp;
618  if(x.innerStride()!=1)
619  {
620  x_tmp.resize(x.rows());
621  x_ptr = x_tmp.data();
622  }
623  for (int j=0; j<rhsCols; ++j)
624  {
625  if(x.innerStride()==1)
626  x_ptr = &x.col(j).coeffRef(0);
627  StorageIndex errorCode = umfpack_solve(UMFPACK_A,
628  mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
629  x_ptr, &b.const_cast_derived().col(j).coeffRef(0),
630  m_numeric, m_control.data(), m_umfpackInfo.data());
631  if(x.innerStride()!=1)
632  x.col(j) = x_tmp;
633  if (errorCode!=0)
634  return false;
635  }
636 
637  return true;
638 }
639 
640 } // end namespace Eigen
641 
642 #endif // EIGEN_UMFPACKSUPPORT_H
Eigen::NumericalIssue
@ NumericalIssue
Definition: Constants.h:443
Eigen
Namespace containing all symbols from the Eigen library.
Definition: LDLT.h:16
Eigen::SparseMatrix< Scalar >
Eigen::UmfPackLU::umfpackControl
const UmfpackControl & umfpackControl() const
Definition: UmfPackSupport.h:421
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition: Constants.h:65
Eigen::Success
@ Success
Definition: Constants.h:441
Eigen::Array< double, UMFPACK_CONTROL, 1 >
Eigen::UmfPackLU::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: UmfPackSupport.h:343
Eigen::UmfPackLU::umfpackControl
UmfpackControl & umfpackControl()
Definition: UmfPackSupport.h:432
Eigen::UmfPackLU::compute
void compute(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:378
Eigen::UmfPackLU::printUmfpackInfo
void printUmfpackInfo()
Definition: UmfPackSupport.h:468
Eigen::UmfPackLU::printUmfpackStatus
void printUmfpackStatus()
Definition: UmfPackSupport.h:478
Eigen::UmfPackLU
A sparse LU factorization and solver based on UmfPack.
Definition: UmfPackSupport.h:289
Eigen::UmfPackLU::printUmfpackControl
void printUmfpackControl()
Definition: UmfPackSupport.h:459
Eigen::SparseSolverBase
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
Eigen::UmfPackLU::analyzePattern
void analyzePattern(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:394
Eigen::Ref< const UmfpackMatrixType, StandardCompressedFormat >
Eigen::UmfPackLU::factorize
void factorize(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:444
Eigen::Matrix< Scalar, Dynamic, 1 >
Eigen::InvalidInput
@ InvalidInput
Definition: Constants.h:448
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
Eigen::UmfPackLU::umfpackFactorizeReturncode
int umfpackFactorizeReturncode() const
Definition: UmfPackSupport.h:409
Eigen::ComputationInfo
ComputationInfo
Definition: Constants.h:439
Eigen::PlainObjectBase::data
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
Definition: PlainObjectBase.h:255
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:42
Generated by   doxygen 1.8.19