GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 11:49:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

Added support for SuperLU's ILU factorization

Browse Source
This commit is contained in:
Peter Román 2009-08-21 11:14:45 +02:00
parent b0aa2520f1
commit 80179e9549
1 changed files with 81 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
Eigen/src/Sparse/SuperLUSupport.h
View File

@ -48,6 +48,29 @@ DECL_GSSVX(SuperLU_C,cgssvx,float,std::complex<float>)
DECL_GSSVX(SuperLU_D,dgssvx,double,double) DECL_GSSVX(SuperLU_D,dgssvx,double,double)
DECL_GSSVX(SuperLU_Z,zgssvx,double,std::complex<double>) DECL_GSSVX(SuperLU_Z,zgssvx,double,std::complex<double>)
// similarly for the incomplete factorization using gsisx
#define DECL_GSISX(NAMESPACE,FNAME,FLOATTYPE,KEYTYPE) \
inline float SuperLU_gsisx(superlu_options_t *options, SuperMatrix *A, \
int *perm_c, int *perm_r, int *etree, char *equed, \
FLOATTYPE *R, FLOATTYPE *C, SuperMatrix *L, \
SuperMatrix *U, void *work, int lwork, \
SuperMatrix *B, SuperMatrix *X, \
FLOATTYPE *recip_pivot_growth, \
FLOATTYPE *rcond, \
SuperLUStat_t *stats, int *info, KEYTYPE) { \
using namespace NAMESPACE; \
mem_usage_t mem_usage; \
NAMESPACE::FNAME(options, A, perm_c, perm_r, etree, equed, R, C, L, \
U, work, lwork, B, X, recip_pivot_growth, rcond, \
&mem_usage, stats, info); \
return mem_usage.for_lu; /* bytes used by the factor storage */ \
}
DECL_GSISX(SuperLU_S,sgsisx,float,float)
DECL_GSISX(SuperLU_C,cgsisx,float,std::complex<float>)
DECL_GSISX(SuperLU_D,dgsisx,double,double)
DECL_GSISX(SuperLU_Z,zgsisx,double,std::complex<double>)
template<typename MatrixType> template<typename MatrixType>
struct SluMatrixMapHelper; struct SluMatrixMapHelper;
@ -373,7 +396,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
m_sluA = m_matrix.asSluMatrix(); m_sluA = m_matrix.asSluMatrix();
memset(&m_sluL,0,sizeof m_sluL); memset(&m_sluL,0,sizeof m_sluL);
memset(&m_sluU,0,sizeof m_sluU); memset(&m_sluU,0,sizeof m_sluU);
m_sluEqued = 'B'; //m_sluEqued = 'B';
int info = 0; int info = 0;
m_p.resize(size); m_p.resize(size);
@ -395,6 +418,29 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
m_sluX = m_sluB; m_sluX = m_sluB;
StatInit(&m_sluStat); StatInit(&m_sluStat);
if (m_flags&IncompleteFactorization)
{
ilu_set_default_options(&m_sluOptions);
// no attempt to preserve column sum
m_sluOptions.ILU_MILU = SILU;
// only basic ILU(k) support -- no direct control over memory consumption
// better to use ILU_DropRule = DROP_BASIC | DROP_AREA
// and set ILU_FillFactor to max memory growth
m_sluOptions.ILU_DropRule = DROP_BASIC;
m_sluOptions.ILU_DropTol = Base::m_precision;
SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&m_sluStat, &info, Scalar());
}
else
{
SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0], SuperLU_gssvx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0], &m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU, &m_sluL, &m_sluU,
@ -403,6 +449,7 @@ void SparseLU<MatrixType,SuperLU>::compute(const MatrixType& a)
&recip_pivot_gross, &rcond, &recip_pivot_gross, &rcond,
&ferr, &berr, &ferr, &berr,
&m_sluStat, &info, Scalar()); &m_sluStat, &info, Scalar());
}
StatFree(&m_sluStat); StatFree(&m_sluStat);
m_extractedDataAreDirty = true; m_extractedDataAreDirty = true;
@ -440,6 +487,19 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b,
StatInit(&m_sluStat); StatInit(&m_sluStat);
int info = 0; int info = 0;
RealScalar recip_pivot_gross, rcond; RealScalar recip_pivot_gross, rcond;
if (m_flags&IncompleteFactorization)
{
SuperLU_gsisx(&m_sluOptions, &m_sluA, m_q.data(), m_p.data(), &m_sluEtree[0],
&m_sluEqued, &m_sluRscale[0], &m_sluCscale[0],
&m_sluL, &m_sluU,
NULL, 0,
&m_sluB, &m_sluX,
&recip_pivot_gross, &rcond,
&m_sluStat, &info, Scalar());
}
else
{
SuperLU_gssvx( SuperLU_gssvx(
&m_sluOptions, &m_sluA, &m_sluOptions, &m_sluA,
m_q.data(), m_p.data(), m_q.data(), m_p.data(),
@ -451,6 +511,7 @@ bool SparseLU<MatrixType,SuperLU>::solve(const MatrixBase<BDerived> &b,
&recip_pivot_gross, &rcond, &recip_pivot_gross, &rcond,
&m_sluFerr[0], &m_sluBerr[0], &m_sluFerr[0], &m_sluBerr[0],
&m_sluStat, &info, Scalar()); &m_sluStat, &info, Scalar());
}
StatFree(&m_sluStat); StatFree(&m_sluStat);
// reset to previous state // reset to previous state