10 #ifndef EIGEN_UMFPACKSUPPORT_H 11 #define EIGEN_UMFPACKSUPPORT_H 19 inline void umfpack_free_numeric(
void **Numeric,
double)
20 { umfpack_di_free_numeric(Numeric); *Numeric = 0; }
22 inline void umfpack_free_numeric(
void **Numeric, std::complex<double>)
23 { umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
25 inline void umfpack_free_symbolic(
void **Symbolic,
double)
26 { umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
28 inline void umfpack_free_symbolic(
void **Symbolic, std::complex<double>)
29 { umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
31 inline int umfpack_symbolic(
int n_row,
int n_col,
32 const int Ap[],
const int Ai[],
const double Ax[],
void **Symbolic,
33 const double Control [UMFPACK_CONTROL],
double Info [UMFPACK_INFO])
35 return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
38 inline int umfpack_symbolic(
int n_row,
int n_col,
39 const int Ap[],
const int Ai[],
const std::complex<double> Ax[],
void **Symbolic,
40 const double Control [UMFPACK_CONTROL],
double Info [UMFPACK_INFO])
42 return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
45 inline int umfpack_numeric(
const int Ap[],
const int Ai[],
const double Ax[],
46 void *Symbolic,
void **Numeric,
47 const double Control[UMFPACK_CONTROL],
double Info [UMFPACK_INFO])
49 return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
52 inline int umfpack_numeric(
const int Ap[],
const int Ai[],
const std::complex<double> Ax[],
53 void *Symbolic,
void **Numeric,
54 const double Control[UMFPACK_CONTROL],
double Info [UMFPACK_INFO])
56 return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
59 inline int umfpack_solve(
int sys,
const int Ap[],
const int Ai[],
const double Ax[],
60 double X[],
const double B[],
void *Numeric,
61 const double Control[UMFPACK_CONTROL],
double Info[UMFPACK_INFO])
63 return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
66 inline int umfpack_solve(
int sys,
const int Ap[],
const int Ai[],
const std::complex<double> Ax[],
67 std::complex<double> X[],
const std::complex<double> B[],
void *Numeric,
68 const double Control[UMFPACK_CONTROL],
double Info[UMFPACK_INFO])
70 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);
73 inline int umfpack_get_lunz(
int *lnz,
int *unz,
int *n_row,
int *n_col,
int *nz_udiag,
void *Numeric,
double)
75 return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
78 inline int umfpack_get_lunz(
int *lnz,
int *unz,
int *n_row,
int *n_col,
int *nz_udiag,
void *Numeric, std::complex<double>)
80 return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
83 inline int umfpack_get_numeric(
int Lp[],
int Lj[],
double Lx[],
int Up[],
int Ui[],
double Ux[],
84 int P[],
int Q[],
double Dx[],
int *do_recip,
double Rs[],
void *Numeric)
86 return umfpack_di_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
89 inline int umfpack_get_numeric(
int Lp[],
int Lj[], std::complex<double> Lx[],
int Up[],
int Ui[], std::complex<double> Ux[],
90 int P[],
int Q[], std::complex<double> Dx[],
int *do_recip,
double Rs[],
void *Numeric)
92 double& lx0_real = numext::real_ref(Lx[0]);
93 double& ux0_real = numext::real_ref(Ux[0]);
94 double& dx0_real = numext::real_ref(Dx[0]);
95 return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
96 Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
99 inline int umfpack_get_determinant(
double *Mx,
double *Ex,
void *NumericHandle,
double User_Info [UMFPACK_INFO])
101 return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
104 inline int umfpack_get_determinant(std::complex<double> *Mx,
double *Ex,
void *NumericHandle,
double User_Info [UMFPACK_INFO])
106 double& mx_real = numext::real_ref(*Mx);
107 return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
112 template<
typename Scalar,
int Options,
typename StorageIndex>
115 template<
typename Scalar,
int Options,
typename StorageIndex>
133 template<
typename _MatrixType>
137 typedef _MatrixType MatrixType;
139 typedef typename MatrixType::RealScalar RealScalar;
140 typedef typename MatrixType::Index Index;
159 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,
Scalar());
160 if(m_numeric) umfpack_free_numeric(&m_numeric,
Scalar());
163 inline Index rows()
const {
return m_copyMatrix.rows(); }
164 inline Index cols()
const {
return m_copyMatrix.cols(); }
173 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
177 inline const LUMatrixType& matrixL()
const 179 if (m_extractedDataAreDirty) extractData();
183 inline const LUMatrixType& matrixU()
const 185 if (m_extractedDataAreDirty) extractData();
189 inline const IntColVectorType& permutationP()
const 191 if (m_extractedDataAreDirty) extractData();
195 inline const IntRowVectorType& permutationQ()
const 197 if (m_extractedDataAreDirty) extractData();
205 template<
typename InputMatrixType>
208 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,
Scalar());
209 if(m_numeric) umfpack_free_numeric(&m_numeric,
Scalar());
210 grapInput(matrix.derived());
211 analyzePattern_impl();
219 template<
typename Rhs>
222 eigen_assert(m_isInitialized &&
"UmfPackLU is not initialized.");
223 eigen_assert(rows()==b.rows()
224 &&
"UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
232 template<
typename Rhs>
235 eigen_assert(m_isInitialized &&
"UmfPackLU is not initialized.");
236 eigen_assert(rows()==b.
rows()
237 &&
"UmfPackLU::solve(): invalid number of rows of the right hand side matrix b");
247 template<
typename InputMatrixType>
250 if(m_symbolic) umfpack_free_symbolic(&m_symbolic,
Scalar());
251 if(m_numeric) umfpack_free_numeric(&m_numeric,
Scalar());
253 grapInput(matrix.derived());
255 analyzePattern_impl();
264 template<
typename InputMatrixType>
267 eigen_assert(m_analysisIsOk &&
"UmfPackLU: you must first call analyzePattern()");
269 umfpack_free_numeric(&m_numeric,
Scalar());
271 grapInput(matrix.derived());
276 #ifndef EIGEN_PARSED_BY_DOXYGEN 278 template<
typename BDerived,
typename XDerived>
282 Scalar determinant()
const;
284 void extractData()
const;
291 m_isInitialized =
false;
297 m_extractedDataAreDirty =
true;
300 template<
typename InputMatrixType>
303 m_copyMatrix.resize(mat.rows(), mat.cols());
304 if( ((MatrixType::Flags&
RowMajorBit)==RowMajorBit) ||
sizeof(
typename MatrixType::Index)!=
sizeof(
int) || !mat.isCompressed() )
308 m_outerIndexPtr = m_copyMatrix.outerIndexPtr();
309 m_innerIndexPtr = m_copyMatrix.innerIndexPtr();
310 m_valuePtr = m_copyMatrix.valuePtr();
314 m_outerIndexPtr = mat.outerIndexPtr();
315 m_innerIndexPtr = mat.innerIndexPtr();
316 m_valuePtr = mat.valuePtr();
320 template<
typename InputMatrixType>
324 m_outerIndexPtr = m_copyMatrix.outerIndexPtr();
325 m_innerIndexPtr = m_copyMatrix.innerIndexPtr();
326 m_valuePtr = m_copyMatrix.valuePtr();
329 template<
typename InputMatrixType>
330 void grapInput(
const InputMatrixType& mat)
335 void analyzePattern_impl()
338 errorCode = umfpack_symbolic(m_copyMatrix.rows(), m_copyMatrix.cols(), m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
341 m_isInitialized =
true;
343 m_analysisIsOk =
true;
344 m_factorizationIsOk =
false;
345 m_extractedDataAreDirty =
true;
348 void factorize_impl()
351 errorCode = umfpack_numeric(m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
352 m_symbolic, &m_numeric, 0, 0);
355 m_factorizationIsOk =
true;
356 m_extractedDataAreDirty =
true;
360 mutable LUMatrixType m_l;
361 mutable LUMatrixType m_u;
362 mutable IntColVectorType m_p;
363 mutable IntRowVectorType m_q;
365 UmfpackMatrixType m_copyMatrix;
366 const Scalar* m_valuePtr;
367 const int* m_outerIndexPtr;
368 const int* m_innerIndexPtr;
373 bool m_isInitialized;
374 int m_factorizationIsOk;
376 mutable bool m_extractedDataAreDirty;
383 template<
typename MatrixType>
386 if (m_extractedDataAreDirty)
389 int lnz, unz, rows, cols, nz_udiag;
390 umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric,
Scalar());
393 m_l.resize(rows,(std::min)(rows,cols));
394 m_l.resizeNonZeros(lnz);
396 m_u.resize((std::min)(rows,cols),cols);
397 m_u.resizeNonZeros(unz);
403 umfpack_get_numeric(m_l.outerIndexPtr(), m_l.innerIndexPtr(), m_l.valuePtr(),
404 m_u.outerIndexPtr(), m_u.innerIndexPtr(), m_u.valuePtr(),
405 m_p.data(), m_q.data(), 0, 0, 0, m_numeric);
407 m_extractedDataAreDirty =
false;
411 template<
typename MatrixType>
415 umfpack_get_determinant(&det, 0, m_numeric, 0);
419 template<
typename MatrixType>
420 template<
typename BDerived,
typename XDerived>
423 const int rhsCols = b.cols();
424 eigen_assert((BDerived::Flags&
RowMajorBit)==0 &&
"UmfPackLU backend does not support non col-major rhs yet");
425 eigen_assert((XDerived::Flags&RowMajorBit)==0 &&
"UmfPackLU backend does not support non col-major result yet");
426 eigen_assert(b.derived().data() != x.derived().data() &&
" Umfpack does not support inplace solve");
429 for (
int j=0; j<rhsCols; ++j)
431 errorCode = umfpack_solve(UMFPACK_A,
432 m_outerIndexPtr, m_innerIndexPtr, m_valuePtr,
433 &x.
col(j).coeffRef(0), &b.const_cast_derived().
col(j).coeffRef(0), m_numeric, 0, 0);
444 template<
typename _MatrixType,
typename Rhs>
449 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
451 template<
typename Dest>
void evalTo(Dest& dst)
const 453 dec()._solve(rhs(),dst);
457 template<
typename _MatrixType,
typename Rhs>
462 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
464 template<
typename Dest>
void evalTo(Dest& dst)
const 466 this->defaultEvalTo(dst);
474 #endif // EIGEN_UMFPACKSUPPORT_H Eigen::Matrix< Scalar, n, 1 > Vector
A vector of length n.
Definition: FlexibleKalmanBase.h:98
Definition: ForwardDeclarations.h:124
ColXpr col(Index i)
Definition: DenseBase.h:733
A sparse LU factorization and solver based on UmfPack.
Definition: UmfPackSupport.h:134
void factorize(const InputMatrixType &matrix)
Performs a numeric decomposition of matrix.
Definition: UmfPackSupport.h:265
A versatible sparse matrix representation.
Definition: SparseMatrix.h:85
void compute(const InputMatrixType &matrix)
Computes the sparse Cholesky decomposition of matrix Note that the matrix should be column-major...
Definition: UmfPackSupport.h:206
void analyzePattern(const InputMatrixType &matrix)
Performs a symbolic decomposition on the sparcity of matrix.
Definition: UmfPackSupport.h:248
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: TestIMU_Common.h:87
Definition: SparseSolve.h:18
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: UmfPackSupport.h:171
const unsigned int RowMajorBit
for a matrix, this means that the storage order is row-major.
Definition: Constants.h:53
Index rows() const
Definition: SparseMatrixBase.h:159
Base class of any sparse matrices or sparse expressions.
Definition: ForwardDeclarations.h:239
const internal::sparse_solve_retval< UmfPackLU, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
Definition: UmfPackSupport.h:233
Definition: UmfPackSupport.h:111
The provided data did not satisfy the prerequisites.
Definition: Constants.h:378
The inputs are invalid, or the algorithm has been improperly called.
Definition: Constants.h:383
Computation was successful.
Definition: Constants.h:376
Definition: BandTriangularSolver.h:13
Definition: SparseSolve.h:17
Definition: ForwardDeclarations.h:125
const internal::solve_retval< UmfPackLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: UmfPackSupport.h:220
ComputationInfo
Enum for reporting the status of a computation.
Definition: Constants.h:374
Sparse matrix.
Definition: MappedSparseMatrix.h:31
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Definition: osvr_print_tree.cpp:52
double Scalar
Common scalar type.
Definition: FlexibleKalmanBase.h:48