mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-31 19:00:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

UmfPack support: add support for complex<double>

Browse Source
This commit is contained in:
Gael Guennebaud 2008-10-20 00:39:11 +00:00
parent 3a231c2349
commit f44316e5f8
1 changed files with 51 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
Eigen/src/Sparse/UmfPackSupport.h
View File

@ -25,6 +25,51 @@
#ifndef EIGEN_UMFPACKSUPPORT_H
#define EIGEN_UMFPACKSUPPORT_H
/* TODO extract L, extrac U, compute det, etc... */
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])
{
return umfpack_di_symbolic(n_row,n_col,Ap,Ai,Ax,Symbolic,Control,Info);
}
inline int 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])
{
return umfpack_zi_symbolic(n_row,n_col,Ap,Ai,&Ax[0].real(),0,Symbolic,Control,Info);
}
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])
{
return umfpack_di_numeric(Ap,Ai,Ax,Symbolic,Numeric,Control,Info);
}
inline int 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])
{
return umfpack_zi_numeric(Ap,Ai,&Ax[0].real(),0,Symbolic,Numeric,Control,Info);
}
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])
{
return umfpack_di_solve(sys,Ap,Ai,Ax,X,B,Numeric,Control,Info);
}
inline int 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])
{
return umfpack_zi_solve(sys,Ap,Ai,&Ax[0].real(),0,&X[0].real(),0,&B[0].real(),0,Numeric,Control,Info);
}
template<typename MatrixType>
class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
{
@ -69,8 +114,9 @@ class SparseLU<MatrixType,UmfPack> : public SparseLU<MatrixType>
template<typename MatrixType>
void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
{
const int size = a.rows();
ei_assert((MatrixType::Flags&RowMajorBit)==0);
const int rows = a.rows();
const int cols = a.cols();
ei_assert((MatrixType::Flags&RowMajorBit)==0 && "Row major matrices are not supported yet");
m_matrixRef = &a;
@ -79,10 +125,10 @@ void SparseLU<MatrixType,UmfPack>::compute(const MatrixType& a)
void* symbolic;
int errorCode = 0;
errorCode = umfpack_di_symbolic(size, size, a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
errorCode = umfpack_symbolic(rows, cols, a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
&symbolic, 0, 0);
if (errorCode==0)
errorCode = umfpack_di_numeric(a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
errorCode = umfpack_numeric(a._outerIndexPtr(), a._innerIndexPtr(), a._valuePtr(),
symbolic, &m_numeric, 0, 0);
umfpack_di_free_symbolic(&symbolic);
@ -119,7 +165,7 @@ bool SparseLU<MatrixType,UmfPack>::solve(const MatrixBase<BDerived> &b, MatrixBa
int errorCode;
for (int j=0; j<rhsCols; ++j)
{
errorCode = umfpack_di_solve(UMFPACK_A,
errorCode = umfpack_solve(UMFPACK_A,
m_matrixRef->_outerIndexPtr(), m_matrixRef->_innerIndexPtr(), m_matrixRef->_valuePtr(),
&x->col(j).coeffRef(0), &b.const_cast_derived().col(j).coeffRef(0), m_numeric, 0, 0);
if (errorCode!=0)