mirror of
https://gitlab.com/libeigen/eigen.git
synced 2024-12-21 07:19:46 +08:00
Add interface to umfpack_*l_* functions
This commit is contained in:
vhuber
parent
624df50945
commit
baf9a5a776
@ -12,47 +12,92 @@
|
||||
|
||||
namespace Eigen {
|
||||
|
||||
#define UF_long __int64
|
||||
|
||||
/* TODO extract L, extract U, compute det, etc... */
|
||||
|
||||
// generic double/complex<double> wrapper functions:
|
||||
|
||||
|
||||
inline void umfpack_defaults(double control[UMFPACK_CONTROL], double)
|
||||
// Defaults
|
||||
inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, int)
|
||||
{ umfpack_di_defaults(control); }
|
||||
|
||||
inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>)
|
||||
inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, int)
|
||||
{ umfpack_zi_defaults(control); }
|
||||
|
||||
inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double)
|
||||
inline void umfpack_defaults(double control[UMFPACK_CONTROL], double, UF_long)
|
||||
{ umfpack_dl_defaults(control); }
|
||||
|
||||
inline void umfpack_defaults(double control[UMFPACK_CONTROL], std::complex<double>, UF_long)
|
||||
{ umfpack_zl_defaults(control); }
|
||||
|
||||
// Report info
|
||||
inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, int)
|
||||
{ umfpack_di_report_info(control, info);}
|
||||
|
||||
inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>)
|
||||
inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, int)
|
||||
{ umfpack_zi_report_info(control, info);}
|
||||
|
||||
inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double)
|
||||
inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], double, UF_long)
|
||||
{ umfpack_dl_report_info(control, info);}
|
||||
|
||||
inline void umfpack_report_info(double control[UMFPACK_CONTROL], double info[UMFPACK_INFO], std::complex<double>, UF_long)
|
||||
{ umfpack_zl_report_info(control, info);}
|
||||
|
||||
// Report status
|
||||
inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, int)
|
||||
{ umfpack_di_report_status(control, status);}
|
||||
|
||||
inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>)
|
||||
inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, int)
|
||||
{ umfpack_zi_report_status(control, status);}
|
||||
|
||||
inline void umfpack_report_control(double control[UMFPACK_CONTROL], double)
|
||||
inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, double, UF_long)
|
||||
{ umfpack_dl_report_status(control, status);}
|
||||
|
||||
inline void umfpack_report_status(double control[UMFPACK_CONTROL], int status, std::complex<double>, UF_long)
|
||||
{ umfpack_zl_report_status(control, status);}
|
||||
|
||||
// report control
|
||||
inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, int)
|
||||
{ umfpack_di_report_control(control);}
|
||||
|
||||
inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>)
|
||||
inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, int)
|
||||
{ umfpack_zi_report_control(control);}
|
||||
|
||||
inline void umfpack_free_numeric(void **Numeric, double)
|
||||
inline void umfpack_report_control(double control[UMFPACK_CONTROL], double, UF_long)
|
||||
{ umfpack_dl_report_control(control);}
|
||||
|
||||
inline void umfpack_report_control(double control[UMFPACK_CONTROL], std::complex<double>, UF_long)
|
||||
{ umfpack_zl_report_control(control);}
|
||||
|
||||
// Free numeric
|
||||
inline void umfpack_free_numeric(void **Numeric, double, int)
|
||||
{ umfpack_di_free_numeric(Numeric); *Numeric = 0; }
|
||||
|
||||
inline void umfpack_free_numeric(void **Numeric, std::complex<double>)
|
||||
inline void umfpack_free_numeric(void **Numeric, std::complex<double>, int)
|
||||
{ umfpack_zi_free_numeric(Numeric); *Numeric = 0; }
|
||||
|
||||
inline void umfpack_free_symbolic(void **Symbolic, double)
|
||||
inline void umfpack_free_numeric(void **Numeric, double, UF_long)
|
||||
{ umfpack_dl_free_numeric(Numeric); *Numeric = 0; }
|
||||
|
||||
inline void umfpack_free_numeric(void **Numeric, std::complex<double>, UF_long)
|
||||
{ umfpack_zl_free_numeric(Numeric); *Numeric = 0; }
|
||||
|
||||
// Free symbolic
|
||||
inline void umfpack_free_symbolic(void **Symbolic, double, int)
|
||||
{ umfpack_di_free_symbolic(Symbolic); *Symbolic = 0; }
|
||||
|
||||
inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>)
|
||||
inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, int)
|
||||
{ umfpack_zi_free_symbolic(Symbolic); *Symbolic = 0; }
|
||||
|
||||
inline void umfpack_free_symbolic(void **Symbolic, double, UF_long)
|
||||
{ umfpack_dl_free_symbolic(Symbolic); *Symbolic = 0; }
|
||||
|
||||
inline void umfpack_free_symbolic(void **Symbolic, std::complex<double>, UF_long)
|
||||
{ umfpack_zl_free_symbolic(Symbolic); *Symbolic = 0; }
|
||||
|
||||
// Symbolic
|
||||
inline int 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])
|
||||
@ -66,7 +111,21 @@ inline int umfpack_symbolic(int n_row,int n_col,
|
||||
{
|
||||
return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
|
||||
}
|
||||
inline int umfpack_symbolic(UF_long n_row,UF_long n_col,
|
||||
const UF_long Ap[], const UF_long Ai[], const double Ax[], void **Symbolic,
|
||||
const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
|
||||
{
|
||||
return umfpack_dl_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
|
||||
}
|
||||
|
||||
inline int umfpack_symbolic(UF_long n_row,UF_long n_col,
|
||||
const UF_long Ap[], const UF_long Ai[], const std::complex<double> Ax[], void **Symbolic,
|
||||
const double Control [UMFPACK_CONTROL], double Info [UMFPACK_INFO])
|
||||
{
|
||||
return umfpack_zl_symbolic(n_row,n_col,Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Control,Info);
|
||||
}
|
||||
|
||||
// Numeric
|
||||
inline int umfpack_numeric( const int Ap[], const int Ai[], const double Ax[],
|
||||
void *Symbolic, void **Numeric,
|
||||
const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
|
||||
@ -80,7 +139,21 @@ inline int umfpack_numeric( const int Ap[], const int Ai[], const std::complex<d
|
||||
{
|
||||
return umfpack_zi_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
|
||||
}
|
||||
inline int umfpack_numeric( const UF_long Ap[], const UF_long Ai[], const double Ax[],
|
||||
void *Symbolic, void **Numeric,
|
||||
const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
|
||||
{
|
||||
return umfpack_dl_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
|
||||
}
|
||||
|
||||
inline int umfpack_numeric( const UF_long Ap[], const UF_long Ai[], const std::complex<double> Ax[],
|
||||
void *Symbolic, void **Numeric,
|
||||
const double Control[UMFPACK_CONTROL],double Info [UMFPACK_INFO])
|
||||
{
|
||||
return umfpack_zl_numeric(Ap,Ai,&numext::real_ref(Ax[0]),0,Symbolic,Numeric,Control,Info);
|
||||
}
|
||||
|
||||
// solve
|
||||
inline int 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])
|
||||
@ -95,6 +168,21 @@ inline int umfpack_solve( int sys, const int Ap[], const int Ai[], const std::co
|
||||
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);
|
||||
}
|
||||
|
||||
inline int umfpack_solve( int sys, const UF_long Ap[], const UF_long Ai[], const double Ax[],
|
||||
double X[], const double B[], void *Numeric,
|
||||
const double Control[UMFPACK_CONTROL], double Info[UMFPACK_INFO])
|
||||
{
|
||||
return umfpack_dl_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
|
||||
}
|
||||
|
||||
inline int umfpack_solve( int sys, const UF_long Ap[], const UF_long 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])
|
||||
{
|
||||
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);
|
||||
}
|
||||
|
||||
// Get Lunz
|
||||
inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_udiag, void *Numeric, double)
|
||||
{
|
||||
return umfpack_di_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
|
||||
@ -105,6 +193,17 @@ inline int umfpack_get_lunz(int *lnz, int *unz, int *n_row, int *n_col, int *nz_
|
||||
return umfpack_zi_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
|
||||
}
|
||||
|
||||
inline int umfpack_get_lunz(UF_long *lnz, UF_long *unz, UF_long *n_row, UF_long *n_col, UF_long *nz_udiag, void *Numeric, double)
|
||||
{
|
||||
return umfpack_dl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
|
||||
}
|
||||
|
||||
inline int umfpack_get_lunz(UF_long *lnz, UF_long *unz, UF_long *n_row, UF_long *n_col, UF_long *nz_udiag, void *Numeric, std::complex<double>)
|
||||
{
|
||||
return umfpack_zl_get_lunz(lnz,unz,n_row,n_col,nz_udiag,Numeric);
|
||||
}
|
||||
|
||||
// Get Numeric
|
||||
inline int 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)
|
||||
{
|
||||
@ -120,18 +219,45 @@ inline int umfpack_get_numeric(int Lp[], int Lj[], std::complex<double> Lx[], in
|
||||
return umfpack_zi_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
|
||||
Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
|
||||
}
|
||||
inline int umfpack_get_numeric(UF_long Lp[], UF_long Lj[], double Lx[], UF_long Up[], UF_long Ui[], double Ux[],
|
||||
UF_long P[], UF_long Q[], double Dx[], UF_long *do_recip, double Rs[], void *Numeric)
|
||||
{
|
||||
return umfpack_dl_get_numeric(Lp,Lj,Lx,Up,Ui,Ux,P,Q,Dx,do_recip,Rs,Numeric);
|
||||
}
|
||||
|
||||
inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
|
||||
inline int umfpack_get_numeric(UF_long Lp[], UF_long Lj[], std::complex<double> Lx[], UF_long Up[], UF_long Ui[], std::complex<double> Ux[],
|
||||
UF_long P[], UF_long Q[], std::complex<double> Dx[], UF_long *do_recip, double Rs[], void *Numeric)
|
||||
{
|
||||
double& lx0_real = numext::real_ref(Lx[0]);
|
||||
double& ux0_real = numext::real_ref(Ux[0]);
|
||||
double& dx0_real = numext::real_ref(Dx[0]);
|
||||
return umfpack_zl_get_numeric(Lp,Lj,Lx?&lx0_real:0,0,Up,Ui,Ux?&ux0_real:0,0,P,Q,
|
||||
Dx?&dx0_real:0,0,do_recip,Rs,Numeric);
|
||||
}
|
||||
|
||||
// Get Determinant
|
||||
inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int)
|
||||
{
|
||||
return umfpack_di_get_determinant(Mx,Ex,NumericHandle,User_Info);
|
||||
}
|
||||
|
||||
inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO])
|
||||
inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], int)
|
||||
{
|
||||
double& mx_real = numext::real_ref(*Mx);
|
||||
return umfpack_zi_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
|
||||
}
|
||||
|
||||
inline int umfpack_get_determinant(double *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], UF_long)
|
||||
{
|
||||
return umfpack_dl_get_determinant(Mx,Ex,NumericHandle,User_Info);
|
||||
}
|
||||
|
||||
inline int umfpack_get_determinant(std::complex<double> *Mx, double *Ex, void *NumericHandle, double User_Info [UMFPACK_INFO], UF_long)
|
||||
{
|
||||
double& mx_real = numext::real_ref(*Mx);
|
||||
return umfpack_zl_get_determinant(&mx_real,0,Ex,NumericHandle,User_Info);
|
||||
}
|
||||
|
||||
|
||||
/** \ingroup UmfPackSupport_Module
|
||||
* \brief A sparse LU factorization and solver based on UmfPack
|
||||
@ -164,7 +290,7 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
typedef Matrix<int, 1, MatrixType::ColsAtCompileTime> IntRowVectorType;
|
||||
typedef Matrix<int, MatrixType::RowsAtCompileTime, 1> IntColVectorType;
|
||||
typedef SparseMatrix<Scalar> LUMatrixType;
|
||||
typedef SparseMatrix<Scalar,ColMajor,int> UmfpackMatrixType;
|
||||
typedef SparseMatrix<Scalar,ColMajor,StorageIndex> UmfpackMatrixType;
|
||||
typedef Ref<const UmfpackMatrixType, StandardCompressedFormat> UmfpackMatrixRef;
|
||||
enum {
|
||||
ColsAtCompileTime = MatrixType::ColsAtCompileTime,
|
||||
@ -192,8 +318,8 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
|
||||
~UmfPackLU()
|
||||
{
|
||||
if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
|
||||
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
|
||||
if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(), StorageIndex());
|
||||
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(), StorageIndex());
|
||||
}
|
||||
|
||||
inline Index rows() const { return mp_matrix.rows(); }
|
||||
@ -241,8 +367,8 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
template<typename InputMatrixType>
|
||||
void compute(const InputMatrixType& matrix)
|
||||
{
|
||||
if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
|
||||
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
|
||||
if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
|
||||
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
|
||||
grab(matrix.derived());
|
||||
analyzePattern_impl();
|
||||
factorize_impl();
|
||||
@ -257,8 +383,8 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
template<typename InputMatrixType>
|
||||
void analyzePattern(const InputMatrixType& matrix)
|
||||
{
|
||||
if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar());
|
||||
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar());
|
||||
if(m_symbolic) umfpack_free_symbolic(&m_symbolic,Scalar(),StorageIndex());
|
||||
if(m_numeric) umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
|
||||
|
||||
grab(matrix.derived());
|
||||
|
||||
@ -309,7 +435,7 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
{
|
||||
eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
|
||||
if(m_numeric)
|
||||
umfpack_free_numeric(&m_numeric,Scalar());
|
||||
umfpack_free_numeric(&m_numeric,Scalar(),StorageIndex());
|
||||
|
||||
grab(matrix.derived());
|
||||
|
||||
@ -322,7 +448,7 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
*/
|
||||
void printUmfpackControl()
|
||||
{
|
||||
umfpack_report_control(m_control.data(), Scalar());
|
||||
umfpack_report_control(m_control.data(), Scalar(),StorageIndex());
|
||||
}
|
||||
|
||||
/** Prints statistics collected by UmfPack.
|
||||
@ -332,7 +458,7 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
void printUmfpackInfo()
|
||||
{
|
||||
eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
|
||||
umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar());
|
||||
umfpack_report_info(m_control.data(), m_umfpackInfo.data(), Scalar(),StorageIndex());
|
||||
}
|
||||
|
||||
/** Prints the status of the previous factorization operation performed by UmfPack (symbolic or numerical factorization).
|
||||
@ -341,7 +467,7 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
*/
|
||||
void printUmfpackStatus() {
|
||||
eigen_assert(m_analysisIsOk && "UmfPackLU: you must first call analyzePattern()");
|
||||
umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar());
|
||||
umfpack_report_status(m_control.data(), m_fact_errorCode, Scalar(),StorageIndex());
|
||||
}
|
||||
|
||||
/** \internal */
|
||||
@ -362,13 +488,13 @@ class UmfPackLU : public SparseSolverBase<UmfPackLU<_MatrixType> >
|
||||
m_symbolic = 0;
|
||||
m_extractedDataAreDirty = true;
|
||||
|
||||
umfpack_defaults(m_control.data(), Scalar());
|
||||
umfpack_defaults(m_control.data(), Scalar(),StorageIndex());
|
||||
}
|
||||
|
||||
void analyzePattern_impl()
|
||||
{
|
||||
m_fact_errorCode = umfpack_symbolic(internal::convert_index<int>(mp_matrix.rows()),
|
||||
internal::convert_index<int>(mp_matrix.cols()),
|
||||
m_fact_errorCode = umfpack_symbolic(internal::convert_index<StorageIndex>(mp_matrix.rows()),
|
||||
internal::convert_index<StorageIndex>(mp_matrix.cols()),
|
||||
mp_matrix.outerIndexPtr(), mp_matrix.innerIndexPtr(), mp_matrix.valuePtr(),
|
||||
&m_symbolic, m_control.data(), m_umfpackInfo.data());
|
||||
|
||||
@ -438,7 +564,7 @@ void UmfPackLU<MatrixType>::extractData() const
|
||||
if (m_extractedDataAreDirty)
|
||||
{
|
||||
// get size of the data
|
||||
int lnz, unz, rows, cols, nz_udiag;
|
||||
StorageIndex lnz, unz, rows, cols, nz_udiag;
|
||||
umfpack_get_lunz(&lnz, &unz, &rows, &cols, &nz_udiag, m_numeric, Scalar());
|
||||
|
||||
// allocate data
|
||||
@ -464,7 +590,7 @@ template<typename MatrixType>
|
||||
typename UmfPackLU<MatrixType>::Scalar UmfPackLU<MatrixType>::determinant() const
|
||||
{
|
||||
Scalar det;
|
||||
umfpack_get_determinant(&det, 0, m_numeric, 0);
|
||||
umfpack_get_determinant(&det, 0, m_numeric, 0, StorageIndex());
|
||||
return det;
|
||||
}
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user