10#ifndef EIGEN_UMFPACKSUPPORT_H
11#define EIGEN_UMFPACKSUPPORT_H
15#ifndef SuiteSparse_long
17#define SuiteSparse_long UF_long
19#error neither SuiteSparse_long nor UF_long are defined
31inline void umfpack_defaults(
double control[UMFPACK_CONTROL],
double,
int)
32{ umfpack_di_defaults(control); }
34inline void umfpack_defaults(
double control[UMFPACK_CONTROL], std::complex<double>,
int)
35{ umfpack_zi_defaults(control); }
37inline void umfpack_defaults(
double control[UMFPACK_CONTROL],
double, SuiteSparse_long)
38{ umfpack_dl_defaults(control); }
40inline void umfpack_defaults(
double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long)
41{ umfpack_zl_defaults(control); }
44inline void umfpack_report_info(
double control[UMFPACK_CONTROL],
double info[UMFPACK_INFO],
double,
int)
45{ umfpack_di_report_info(control, info);}
47inline void umfpack_report_info(
double control[UMFPACK_CONTROL],
double info[UMFPACK_INFO], std::complex<double>,
int)
48{ umfpack_zi_report_info(control, info);}
50inline void umfpack_report_info(
double control[UMFPACK_CONTROL],
double info[UMFPACK_INFO],
double, SuiteSparse_long)
51{ umfpack_dl_report_info(control, info);}
53inline void umfpack_report_info(
double control[UMFPACK_CONTROL],
double info[UMFPACK_INFO], std::complex<double>, SuiteSparse_long)
54{ umfpack_zl_report_info(control, info);}
57inline void umfpack_report_status(
double control[UMFPACK_CONTROL],
int status,
double,
int)
58{ umfpack_di_report_status(control, status);}
60inline void umfpack_report_status(
double control[UMFPACK_CONTROL],
int status, std::complex<double>,
int)
61{ umfpack_zi_report_status(control, status);}
63inline void umfpack_report_status(
double control[UMFPACK_CONTROL],
int status,
double, SuiteSparse_long)
64{ umfpack_dl_report_status(control, status);}
66inline void umfpack_report_status(
double control[UMFPACK_CONTROL],
int status, std::complex<double>, SuiteSparse_long)
67{ umfpack_zl_report_status(control, status);}
70inline void umfpack_report_control(
double control[UMFPACK_CONTROL],
double,
int)
71{ umfpack_di_report_control(control);}
73inline void umfpack_report_control(
double control[UMFPACK_CONTROL], std::complex<double>,
int)
74{ umfpack_zi_report_control(control);}
76inline void umfpack_report_control(
double control[UMFPACK_CONTROL],
double, SuiteSparse_long)
77{ umfpack_dl_report_control(control);}
79inline void umfpack_report_control(
double control[UMFPACK_CONTROL], std::complex<double>, SuiteSparse_long)
80{ umfpack_zl_report_control(control);}
83inline void umfpack_free_numeric(
void **Numeric,
double,
int)
84{ umfpack_di_free_numeric(Numeric); *Numeric = 0; }
86inline void umfpack_free_numeric(
void **Numeric, std::complex<double>,
int)
87{ umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
89inline void umfpack_free_numeric(
void **Numeric,
double, SuiteSparse_long)
90{ umfpack_dl_free_numeric(Numeric); *Numeric = 0; }
92inline void umfpack_free_numeric(
void **Numeric, std::complex<double>, SuiteSparse_long)
93{ umfpack_zl_free_numeric(Numeric); *Numeric = 0; }
96inline void umfpack_free_symbolic(
void **Symbolic,
double,
int)
97{ umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
99inline void umfpack_free_symbolic(
void **Symbolic, std::complex<double>,
int)
100{ umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
102inline void umfpack_free_symbolic(
void **Symbolic,
double, SuiteSparse_long)
103{ umfpack_dl_free_symbolic(Symbolic); *Symbolic = 0; }
105inline void umfpack_free_symbolic(
void **Symbolic, std::complex<double>, SuiteSparse_long)
106{ umfpack_zl_free_symbolic(Symbolic); *Symbolic = 0; }
109inline 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])
113 return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
116inline 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])
120 return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
122inline 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])
126 return umfpack_dl_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
129inline 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])
133 return umfpack_zl_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
137inline 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])
141 return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
144inline 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])
148 return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
150inline 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])
154 return umfpack_dl_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
157inline 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])
161 return umfpack_zl_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
165inline 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])
169 return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
172inline 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])
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);
179inline 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])
183 return umfpack_dl_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
186inline 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])
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);
194inline int umfpack_get_lunz(
int *lnz,
int *unz,
int *n_row,
int *n_col,
int *nz_udiag,
void *Numeric,
double)
196 return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
199inline int umfpack_get_lunz(
int *lnz,
int *unz,
int *n_row,
int *n_col,
int *nz_udiag,
void *Numeric, std::complex<double>)
201 return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
204inline 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)
207 return umfpack_dl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
210inline 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>)
213 return umfpack_zl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
217inline 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)
220 return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
223inline 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)
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);
232inline 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)
235 return umfpack_dl_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
238inline 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)
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);
249inline int umfpack_get_determinant(
double *Mx,
double *Ex,
void *NumericHandle,
double User_Info [UMFPACK_INFO],
int)
251 return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
254inline int umfpack_get_determinant(std::complex<double> *Mx,
double *Ex,
void *NumericHandle,
double User_Info [UMFPACK_INFO],
int)
256 double& mx_real = numext::real_ref(*Mx);
257 return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
260inline SuiteSparse_long umfpack_get_determinant(
double *Mx,
double *Ex,
void *NumericHandle,
double User_Info [UMFPACK_INFO], SuiteSparse_long)
262 return umfpack_dl_get_determinant(Mx,Ex,NumericHandle,User_Info);
265inline SuiteSparse_long umfpack_get_determinant(std::complex<double> *Mx,
double *Ex,
void *NumericHandle,
double User_Info [UMFPACK_INFO], SuiteSparse_long)
267 double& mx_real = numext::real_ref(*Mx);
268 return umfpack_zl_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
287template<
typename _MatrixType>
292 using Base::m_isInitialized;
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;
306 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
307 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
316 : m_dummy(0,0), mp_matrix(m_dummy)
321 template<
typename InputMatrixType>
322 explicit UmfPackLU(
const InputMatrixType& matrix)
331 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(), StorageIndex());
332 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(), StorageIndex());
335 inline Index rows()
const {
return mp_matrix.rows(); }
336 inline Index cols()
const {
return mp_matrix.cols(); }
345 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
349 inline const LUMatrixType& matrixL()
const
351 if (m_extractedDataAreDirty) extractData();
355 inline const LUMatrixType& matrixU()
const
357 if (m_extractedDataAreDirty) extractData();
361 inline const IntColVectorType& permutationP()
const
363 if (m_extractedDataAreDirty) extractData();
367 inline const IntRowVectorType& permutationQ()
const
369 if (m_extractedDataAreDirty) extractData();
377 template<
typename InputMatrixType>
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();
393 template<
typename InputMatrixType>
396 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
397 if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
399 grab(matrix.derived());
401 analyzePattern_impl();
411 eigen_assert(m_numeric &&
"UmfPackLU: you must first call factorize()");
412 return m_fact_errorCode;
443 template<
typename InputMatrixType>
446 eigen_assert(m_analysisIsOk &&
"UmfPackLU: you must first call analyzePattern()");
448 umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
450 grab(matrix.derived());
461 umfpack_report_control(m_control.
data(), Scalar(),StorageIndex());
470 eigen_assert(m_analysisIsOk &&
"UmfPackLU: you must first call analyzePattern()");
471 umfpack_report_info(m_control.
data(), m_umfpackInfo.
data(), Scalar(),StorageIndex());
479 eigen_assert(m_analysisIsOk &&
"UmfPackLU: you must first call analyzePattern()");
480 umfpack_report_status(m_control.
data(), m_fact_errorCode, Scalar(),StorageIndex());
484 template<
typename BDerived,
typename XDerived>
487 Scalar determinant()
const;
489 void extractData()
const;
496 m_isInitialized =
false;
499 m_extractedDataAreDirty =
true;
501 umfpack_defaults(m_control.
data(), Scalar(),StorageIndex());
504 void analyzePattern_impl()
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());
511 m_isInitialized =
true;
513 m_analysisIsOk =
true;
514 m_factorizationIsOk =
false;
515 m_extractedDataAreDirty =
true;
518 void factorize_impl()
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());
525 m_factorizationIsOk =
true;
526 m_extractedDataAreDirty =
true;
529 template<
typename MatrixDerived>
530 void grab(
const EigenBase<MatrixDerived> &A)
532 mp_matrix.~UmfpackMatrixRef();
533 ::new (&mp_matrix) UmfpackMatrixRef(A.derived());
536 void grab(
const UmfpackMatrixRef &A)
538 if(&(A.derived()) != &mp_matrix)
540 mp_matrix.~UmfpackMatrixRef();
541 ::new (&mp_matrix) UmfpackMatrixRef(A);
546 mutable LUMatrixType m_l;
547 StorageIndex m_fact_errorCode;
548 UmfpackControl m_control;
549 mutable UmfpackInfo m_umfpackInfo;
551 mutable LUMatrixType m_u;
552 mutable IntColVectorType m_p;
553 mutable IntRowVectorType m_q;
555 UmfpackMatrixType m_dummy;
556 UmfpackMatrixRef mp_matrix;
562 int m_factorizationIsOk;
564 mutable bool m_extractedDataAreDirty;
567 UmfPackLU(
const UmfPackLU& ) { }
571template<
typename MatrixType>
572void UmfPackLU<MatrixType>::extractData()
const
574 if (m_extractedDataAreDirty)
577 StorageIndex lnz, unz, rows, cols, nz_udiag;
578 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
581 m_l.resize(rows,(std::min)(rows,cols));
582 m_l.resizeNonZeros(lnz);
584 m_u.resize((std::min)(rows,cols),cols);
585 m_u.resizeNonZeros(unz);
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);
595 m_extractedDataAreDirty =
false;
599template<
typename MatrixType>
600typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant()
const
603 umfpack_get_determinant(&det, 0, m_numeric, 0, StorageIndex());
607template<
typename MatrixType>
608template<
typename BDerived,
typename XDerived>
609bool UmfPackLU<MatrixType>::_solve_impl(
const MatrixBase<BDerived> &b, MatrixBase<XDerived> &x)
const
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");
617 Matrix<Scalar,Dynamic,1> x_tmp;
618 if(x.innerStride()!=1)
620 x_tmp.resize(x.rows());
621 x_ptr = x_tmp.data();
623 for (
int j=0; j<rhsCols; ++j)
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)
General-purpose arrays with easy API for coefficient-wise operations.
Definition: Array.h:47
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
const Scalar * data() const
Definition: PlainObjectBase.h:247
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:283
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
A sparse LU factorization and solver based on UmfPack.
Definition: UmfPackSupport.h:289
void compute(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:378
void factorize(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:444
void printUmfpackControl()
Definition: UmfPackSupport.h:459
const UmfpackControl & umfpackControl() const
Definition: UmfPackSupport.h:421
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: UmfPackSupport.h:343
UmfpackControl & umfpackControl()
Definition: UmfPackSupport.h:432
int umfpackFactorizeReturncode() const
Definition: UmfPackSupport.h:409
void printUmfpackInfo()
Definition: UmfPackSupport.h:468
void analyzePattern(const InputMatrixType &matrix)
Definition: UmfPackSupport.h:394
void printUmfpackStatus()
Definition: UmfPackSupport.h:478
ComputationInfo
Definition: Constants.h:440
@ NumericalIssue
Definition: Constants.h:444
@ InvalidInput
Definition: Constants.h:449
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74