• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
umfpack_defaults(double control[UMFPACK_CONTROL],double)20 inline void umfpack_defaults(double control[UMFPACK_CONTROL], double)
21 { umfpack_di_defaults(control); }
22 
umfpack_defaults(double control[UMFPACK_CONTROL],std::complex<double>)23 inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>)
24 { umfpack_zi_defaults(control); }
25 
umfpack_report_info(double control[UMFPACK_CONTROL],double info[UMFPACK_INFO],double)26 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double)
27 { umfpack_di_report_info(control, info);}
28 
umfpack_report_info(double control[UMFPACK_CONTROL],double info[UMFPACK_INFO],std::complex<double>)29 inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>)
30 { umfpack_zi_report_info(control, info);}
31 
umfpack_report_status(double control[UMFPACK_CONTROL],int status,double)32 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double)
33 { umfpack_di_report_status(control, status);}
34 
umfpack_report_status(double control[UMFPACK_CONTROL],int status,std::complex<double>)35 inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>)
36 { umfpack_zi_report_status(control, status);}
37 
umfpack_report_control(double control[UMFPACK_CONTROL],double)38 inline void umfpack_report_control(double control[UMFPACK_CONTROL], double)
39 { umfpack_di_report_control(control);}
40 
umfpack_report_control(double control[UMFPACK_CONTROL],std::complex<double>)41 inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>)
42 { umfpack_zi_report_control(control);}
43 
umfpack_free_numeric(void ** Numeric,double)44 inline void umfpack_free_numeric(void **Numeric, double)
45 { umfpack_di_free_numeric(Numeric); *Numeric = 0; }
46 
umfpack_free_numeric(void ** Numeric,std::complex<double>)47 inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
48 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
49 
umfpack_free_symbolic(void ** Symbolic,double)50 inline void umfpack_free_symbolic(void **Symbolic, double)
51 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
52 
umfpack_free_symbolic(void ** Symbolic,std::complex<double>)53 inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
54 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
55 
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])56 inline int umfpack_symbolic(int n_row,int n_col,
57                             const int Ap[], const int Ai[], const double Ax[], void **Symbolic,
58                             const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
59 {
60   return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
61 }
62 
umfpack_symbolic(int n_row,int n_col,const int Ap[],const int Ai[],const std::complex<double> Ax[],void ** Symbolic,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])63 inline int umfpack_symbolic(int n_row,int n_col,
64                             const int Ap[], const int Ai[], const std::complex<double> Ax[], void **Symbolic,
65                             const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
66 {
67   return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
68 }
69 
umfpack_numeric(const int Ap[],const int Ai[],const double Ax[],void * Symbolic,void ** Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])70 inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
71                             void *Symbolic, void **Numeric,
72                             const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
73 {
74   return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
75 }
76 
umfpack_numeric(const int Ap[],const int Ai[],const std::complex<double> Ax[],void * Symbolic,void ** Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])77 inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<double> Ax[],
78                             void *Symbolic, void **Numeric,
79                             const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
80 {
81   return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
82 }
83 
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])84 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const double Ax[],
85                           double X[], const double B[], void *Numeric,
86                           const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
87 {
88   return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
89 }
90 
umfpack_solve(int sys,const int Ap[],const int Ai[],const std::complex<double> Ax[],std::complex<double> X[],const std::complex<double> B[],void * Numeric,const double Control[UMFPACK_CONTROL],double Info[UMFPACK_INFO])91 inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::complex<double> Ax[],
92                           std::complex<double> X[], const std::complex<double> B[], void *Numeric,
93                           const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
94 {
95   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);
96 }
97 
umfpack_get_lunz(int * lnz,int * unz,int * n_row,int * n_col,int * nz_udiag,void * Numeric,double)98 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
99 {
100   return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
101 }
102 
umfpack_get_lunz(int * lnz,int * unz,int * n_row,int * n_col,int * nz_udiag,void * Numeric,std::complex<double>)103 inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, std::complex<double>)
104 {
105   return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
106 }
107 
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)108 inline int umfpack_get_numeric(int Lp[], int Lj[], double Lx[], int Up[], int Ui[], double Ux[],
109                                int P[], int Q[], double Dx[], int *do_recip, double Rs[], void *Numeric)
110 {
111   return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
112 }
113 
umfpack_get_numeric(int Lp[],int Lj[],std::complex<double> Lx[],int Up[],int Ui[],std::complex<double> Ux[],int P[],int Q[],std::complex<double> Dx[],int * do_recip,double Rs[],void * Numeric)114 inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], int Up[], int Ui[], std::complex<double> Ux[],
115                                int P[], int Q[], std::complex<double> Dx[], int *do_recip, double Rs[], void *Numeric)
116 {
117   double& lx0_real = numext::real_ref(Lx[0]);
118   double& ux0_real = numext::real_ref(Ux[0]);
119   double& dx0_real = numext::real_ref(Dx[0]);
120   return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
121                                 Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
122 }
123 
umfpack_get_determinant(double * Mx,double * Ex,void * NumericHandle,double User_Info[UMFPACK_INFO])124 inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
125 {
126   return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
127 }
128 
umfpack_get_determinant(std::complex<double> * Mx,double * Ex,void * NumericHandle,double User_Info[UMFPACK_INFO])129 inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
130 {
131   double& mx_real = numext::real_ref(*Mx);
132   return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
133 }
134 
135 
136 /** \ingroup UmfPackSupport_Module
137   * \brief A sparse LU factorization and solver based on UmfPack
138   *
139   * This class allows to solve for A.X = B sparse linear problems via a LU factorization
140   * using the UmfPack library. The sparse matrix A must be squared and full rank.
141   * The vectors or matrices X and B can be either dense or sparse.
142   *
143   * \warning The input matrix A should be in a \b compressed and \b column-major form.
144   * Otherwise an expensive copy will be made. You can call the inexpensive makeCompressed() to get a compressed matrix.
145   * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
146   *
147   * \implsparsesolverconcept
148   *
149   * \sa \ref TutorialSparseSolverConcept, class SparseLU
150   */
151 template<typename _MatrixType>
152 class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
153 {
154   protected:
155     typedef SparseSolverBase<UmfPackLU<_MatrixType> > Base;
156     using Base::m_isInitialized;
157   public:
158     using Base::_solve_impl;
159     typedef _MatrixType MatrixType;
160     typedef typename MatrixType::Scalar Scalar;
161     typedef typename MatrixType::RealScalar RealScalar;
162     typedef typename MatrixType::StorageIndex StorageIndex;
163     typedef Matrix<Scalar,Dynamic,1> Vector;
164     typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
165     typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
166     typedef SparseMatrix<Scalar> LUMatrixType;
167     typedef SparseMatrix<Scalar,ColMajor,int> UmfpackMatrixType;
168     typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef;
169     enum {
170       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
171       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
172     };
173 
174   public:
175 
176     typedef Array<double, UMFPACK_CONTROL, 1> UmfpackControl;
177     typedef Array<double, UMFPACK_INFO, 1> UmfpackInfo;
178 
UmfPackLU()179     UmfPackLU()
180       : m_dummy(0,0), mp_matrix(m_dummy)
181     {
182       init();
183     }
184 
185     template<typename InputMatrixType>
UmfPackLU(const InputMatrixType & matrix)186     explicit UmfPackLU(const InputMatrixType& matrix)
187       : mp_matrix(matrix)
188     {
189       init();
190       compute(matrix);
191     }
192 
~UmfPackLU()193     ~UmfPackLU()
194     {
195       if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
196       if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar());
197     }
198 
rows()199     inline Index rows() const { return mp_matrix.rows(); }
cols()200     inline Index cols() const { return mp_matrix.cols(); }
201 
202     /** \brief Reports whether previous computation was successful.
203       *
204       * \returns \c Success if computation was succesful,
205       *          \c NumericalIssue if the matrix.appears to be negative.
206       */
info()207     ComputationInfo info() const
208     {
209       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
210       return m_info;
211     }
212 
matrixL()213     inline const LUMatrixType& matrixL() const
214     {
215       if (m_extractedDataAreDirty) extractData();
216       return m_l;
217     }
218 
matrixU()219     inline const LUMatrixType& matrixU() const
220     {
221       if (m_extractedDataAreDirty) extractData();
222       return m_u;
223     }
224 
permutationP()225     inline const IntColVectorType& permutationP() const
226     {
227       if (m_extractedDataAreDirty) extractData();
228       return m_p;
229     }
230 
permutationQ()231     inline const IntRowVectorType& permutationQ() const
232     {
233       if (m_extractedDataAreDirty) extractData();
234       return m_q;
235     }
236 
237     /** Computes the sparse Cholesky decomposition of \a matrix
238      *  Note that the matrix should be column-major, and in compressed format for best performance.
239      *  \sa SparseMatrix::makeCompressed().
240      */
241     template<typename InputMatrixType>
compute(const InputMatrixType & matrix)242     void compute(const InputMatrixType& matrix)
243     {
244       if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
245       if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar());
246       grab(matrix.derived());
247       analyzePattern_impl();
248       factorize_impl();
249     }
250 
251     /** Performs a symbolic decomposition on the sparcity of \a matrix.
252       *
253       * This function is particularly useful when solving for several problems having the same structure.
254       *
255       * \sa factorize(), compute()
256       */
257     template<typename InputMatrixType>
analyzePattern(const InputMatrixType & matrix)258     void analyzePattern(const InputMatrixType& matrix)
259     {
260       if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
261       if(m_numeric)  umfpack_free_numeric(&m_numeric,Scalar());
262 
263       grab(matrix.derived());
264 
265       analyzePattern_impl();
266     }
267 
268     /** Provides the return status code returned by UmfPack during the numeric
269       * factorization.
270       *
271       * \sa factorize(), compute()
272       */
umfpackFactorizeReturncode()273     inline int umfpackFactorizeReturncode() const
274     {
275       eigen_assert(m_numeric && "UmfPackLU: you must first call factorize()");
276       return m_fact_errorCode;
277     }
278 
279     /** Provides access to the control settings array used by UmfPack.
280       *
281       * If this array contains NaN's, the default values are used.
282       *
283       * See UMFPACK documentation for details.
284       */
umfpackControl()285     inline const UmfpackControl& umfpackControl() const
286     {
287       return m_control;
288     }
289 
290     /** Provides access to the control settings array used by UmfPack.
291       *
292       * If this array contains NaN's, the default values are used.
293       *
294       * See UMFPACK documentation for details.
295       */
umfpackControl()296     inline UmfpackControl& umfpackControl()
297     {
298       return m_control;
299     }
300 
301     /** Performs a numeric decomposition of \a matrix
302       *
303       * The given matrix must has the same sparcity than the matrix on which the pattern anylysis has been performed.
304       *
305       * \sa analyzePattern(), compute()
306       */
307     template<typename InputMatrixType>
factorize(const InputMatrixType & matrix)308     void factorize(const InputMatrixType& matrix)
309     {
310       eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
311       if(m_numeric)
312         umfpack_free_numeric(&m_numeric,Scalar());
313 
314       grab(matrix.derived());
315 
316       factorize_impl();
317     }
318 
319     /** Prints the current UmfPack control settings.
320       *
321       * \sa umfpackControl()
322       */
umfpackReportControl()323     void umfpackReportControl()
324     {
325       umfpack_report_control(m_control.data(), Scalar());
326     }
327 
328     /** Prints statistics collected by UmfPack.
329       *
330       * \sa analyzePattern(), compute()
331       */
umfpackReportInfo()332     void umfpackReportInfo()
333     {
334       eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
335       umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar());
336     }
337 
338     /** Prints the status of the previous factorization operation performed by UmfPack (symbolic or numerical factorization).
339       *
340       * \sa analyzePattern(), compute()
341       */
umfpackReportStatus()342     void umfpackReportStatus() {
343       eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
344       umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar());
345     }
346 
347     /** \internal */
348     template<typename BDerived,typename XDerived>
349     bool _solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const;
350 
351     Scalar determinant() const;
352 
353     void extractData() const;
354 
355   protected:
356 
init()357     void init()
358     {
359       m_info                  = InvalidInput;
360       m_isInitialized         = false;
361       m_numeric               = 0;
362       m_symbolic              = 0;
363       m_extractedDataAreDirty = true;
364 
365       umfpack_defaults(m_control.data(), Scalar());
366     }
367 
analyzePattern_impl()368     void analyzePattern_impl()
369     {
370       m_fact_errorCode = umfpack_symbolic(internal::convert_index<int>(mp_matrix.rows()),
371                                           internal::convert_index<int>(mp_matrix.cols()),
372                                           mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
373                                           &m_symbolic, m_control.data(), m_umfpackInfo.data());
374 
375       m_isInitialized = true;
376       m_info = m_fact_errorCode ? InvalidInput : Success;
377       m_analysisIsOk = true;
378       m_factorizationIsOk = false;
379       m_extractedDataAreDirty = true;
380     }
381 
factorize_impl()382     void factorize_impl()
383     {
384 
385       m_fact_errorCode = umfpack_numeric(mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
386                                          m_symbolic, &m_numeric, m_control.data(), m_umfpackInfo.data());
387 
388       m_info = m_fact_errorCode == UMFPACK_OK ? Success : NumericalIssue;
389       m_factorizationIsOk = true;
390       m_extractedDataAreDirty = true;
391     }
392 
393     template<typename MatrixDerived>
grab(const EigenBase<MatrixDerived> & A)394     void grab(const EigenBase<MatrixDerived> &A)
395     {
396       mp_matrix.~UmfpackMatrixRef();
397       ::new (&mp_matrix) UmfpackMatrixRef(A.derived());
398     }
399 
grab(const UmfpackMatrixRef & A)400     void grab(const UmfpackMatrixRef &A)
401     {
402       if(&(A.derived()) != &mp_matrix)
403       {
404         mp_matrix.~UmfpackMatrixRef();
405         ::new (&mp_matrix) UmfpackMatrixRef(A);
406       }
407     }
408 
409     // cached data to reduce reallocation, etc.
410     mutable LUMatrixType m_l;
411     int m_fact_errorCode;
412     UmfpackControl m_control;
413     mutable UmfpackInfo m_umfpackInfo;
414 
415     mutable LUMatrixType m_u;
416     mutable IntColVectorType m_p;
417     mutable IntRowVectorType m_q;
418 
419     UmfpackMatrixType m_dummy;
420     UmfpackMatrixRef mp_matrix;
421 
422     void* m_numeric;
423     void* m_symbolic;
424 
425     mutable ComputationInfo m_info;
426     int m_factorizationIsOk;
427     int m_analysisIsOk;
428     mutable bool m_extractedDataAreDirty;
429 
430   private:
UmfPackLU(const UmfPackLU &)431     UmfPackLU(const UmfPackLU& ) { }
432 };
433 
434 
435 template<typename MatrixType>
extractData()436 void UmfPackLU<MatrixType>::extractData() const
437 {
438   if (m_extractedDataAreDirty)
439   {
440     // get size of the data
441     int lnz, unz, rows, cols, nz_udiag;
442     umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
443 
444     // allocate data
445     m_l.resize(rows,(std::min)(rows,cols));
446     m_l.resizeNonZeros(lnz);
447 
448     m_u.resize((std::min)(rows,cols),cols);
449     m_u.resizeNonZeros(unz);
450 
451     m_p.resize(rows);
452     m_q.resize(cols);
453 
454     // extract
455     umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
456                         m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
457                         m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
458 
459     m_extractedDataAreDirty = false;
460   }
461 }
462 
463 template<typename MatrixType>
determinant()464 typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const
465 {
466   Scalar det;
467   umfpack_get_determinant(&det, 0, m_numeric, 0);
468   return det;
469 }
470 
471 template<typename MatrixType>
472 template<typename BDerived,typename XDerived>
_solve_impl(const MatrixBase<BDerived> & b,MatrixBase<XDerived> & x)473 bool UmfPackLU<MatrixType>::_solve_impl(const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x) const
474 {
475   Index rhsCols = b.cols();
476   eigen_assert((BDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major rhs yet");
477   eigen_assert((XDerived::Flags&RowMajorBit)==0 && "UmfPackLU backend does not support non col-major result yet");
478   eigen_assert(b.derived().data() != x.derived().data() && " Umfpack does not support inplace solve");
479 
480   int errorCode;
481   Scalar* x_ptr = 0;
482   Matrix<Scalar,Dynamic,1> x_tmp;
483   if(x.innerStride()!=1)
484   {
485     x_tmp.resize(x.rows());
486     x_ptr = x_tmp.data();
487   }
488   for (int j=0; j<rhsCols; ++j)
489   {
490     if(x.innerStride()==1)
491       x_ptr = &x.col(j).coeffRef(0);
492     errorCode = umfpack_solve(UMFPACK_A,
493         mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
494         x_ptr, &b.const_cast_derived().col(j).coeffRef(0), m_numeric, m_control.data(), m_umfpackInfo.data());
495     if(x.innerStride()!=1)
496       x.col(j) = x_tmp;
497     if (errorCode!=0)
498       return false;
499   }
500 
501   return true;
502 }
503 
504 } // end namespace Eigen
505 
506 #endif // EIGEN_UMFPACKSUPPORT_H
507