DGMRES: fix null rhs, fix restart, fix m_isDeflInitialized for multiple solve

This commit is contained in:
Gael Guennebaud 2018-10-15 23:46:00 +02:00
parent d835a0bf53
commit 2747b98cfc

View File

@ -241,7 +241,18 @@ template<typename Rhs, typename Dest>
void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rhs& rhs, Dest& x,
const Preconditioner& precond) const
{
const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
RealScalar normRhs = rhs.norm();
if(normRhs <= considerAsZero)
{
x.setZero();
m_error = 0;
return;
}
//Initialization
m_isDeflInitialized = false;
Index n = mat.rows();
DenseVector r0(n);
Index nbIts = 0;
@ -249,10 +260,11 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh
m_Hes.resize(m_restart, m_restart);
m_V.resize(n,m_restart+1);
//Initial residual vector and initial norm
x = precond.solve(x);
if(x.squaredNorm()==0)
x = precond.solve(rhs);
r0 = rhs - mat * x;
RealScalar beta = r0.norm();
RealScalar normRhs = rhs.norm();
m_error = beta/normRhs;
if(m_error < m_tolerance)
m_info = Success;
@ -265,8 +277,10 @@ void DGMRES<_MatrixType, _Preconditioner>::dgmres(const MatrixType& mat,const Rh
dgmresCycle(mat, precond, x, r0, beta, normRhs, nbIts);
// Compute the new residual vector for the restart
if (nbIts < m_iterations && m_info == NoConvergence)
r0 = rhs - mat * x;
if (nbIts < m_iterations && m_info == NoConvergence) {
r0 = rhs - mat * x;
beta = r0.norm();
}
}
}