BiCGSTAB: set default guess to 0, and improve restart mechanism by recomputing the accurate residual.

This commit is contained in:
Gael Guennebaud 2015-06-05 14:37:57 +02:00
parent 98a8d43457
commit 56d4ef7ad6

View File

@ -59,20 +59,21 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
VectorType s(n), t(n); VectorType s(n), t(n);
RealScalar tol2 = tol*tol; RealScalar tol2 = tol*tol*rhs_sqnorm;
RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon(); RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon();
Index i = 0; Index i = 0;
Index restarts = 0; Index restarts = 0;
while ( r.squaredNorm()/rhs_sqnorm > tol2 && i<maxIters ) while ( r.squaredNorm() > tol2 && i<maxIters )
{ {
Scalar rho_old = rho; Scalar rho_old = rho;
rho = r0.dot(r); rho = r0.dot(r);
if (abs(rho) < eps2*r0_sqnorm) if (abs(rho) < eps2*r0_sqnorm)
{ {
// The new residual vector became too orthogonal to the arbitrarily choosen direction r0 // The new residual vector became too orthogonal to the arbitrarily chosen direction r0
// Let's restart with a new r0: // Let's restart with a new r0:
r = rhs - mat * x;
r0 = r; r0 = r;
rho = r0_sqnorm = r.squaredNorm(); rho = r0_sqnorm = r.squaredNorm();
if(restarts++ == 0) if(restarts++ == 0)
@ -202,8 +203,8 @@ public:
template<typename Rhs,typename Dest> template<typename Rhs,typename Dest>
void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const
{ {
// x.setZero(); x.resize(this->rows(),b.cols());
x = b; x.setZero();
_solve_with_guess_impl(b,x); _solve_with_guess_impl(b,x);
} }