GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-05-10 23:09:06 +08:00
Code Issues Packages Projects Releases Wiki Activity

Improve BiCGSTAB robustness: fix a divide by zero and allow to restart with a new initial residual reference.

Browse Source
This commit is contained in:
Gael Guennebaud 2013-07-01 11:49:23 +02:00
parent 99bef0957b
commit 22820e950e
1 changed files with 18 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
Eigen/src/IterativeLinearSolvers/BiCGSTAB.h
View File

@ -43,8 +43,9 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
VectorType r = rhs - mat * x; VectorType r = rhs - mat * x;
VectorType r0 = r; VectorType r0 = r;
RealScalar r0_sqnorm = rhs.squaredNorm(); RealScalar r0_sqnorm = r0.squaredNorm();
if(r0_sqnorm == 0) RealScalar rhs_sqnorm = rhs.squaredNorm();
if(rhs_sqnorm == 0)
{ {
x.setZero(); x.setZero();
return true; return true;
@ -61,13 +62,22 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
RealScalar tol2 = tol*tol; RealScalar tol2 = tol*tol;
int i = 0; int i = 0;
int restarts = 0;
while ( r.squaredNorm()/r0_sqnorm > tol2 && i<maxIters ) while ( r.squaredNorm()/rhs_sqnorm > tol2 && i<maxIters )
{ {
Scalar rho_old = rho; Scalar rho_old = rho;
rho = r0.dot(r); rho = r0.dot(r);
if (rho == Scalar(0)) return false; /* New search directions cannot be found */ if (internal::isMuchSmallerThan(rho,r0_sqnorm))
{
// The new residual vector became too orthogonal to the arbitrarily choosen direction r0
// Let's restart with a new r0:
r0 = r;
r0_sqnorm = rho = r0.dot(r);
if(restarts++ == 0)
i = 0;
}
Scalar beta = (rho/rho_old) * (alpha / w); Scalar beta = (rho/rho_old) * (alpha / w);
p = r + beta * (p - w * v); p = r + beta * (p - w * v);
@ -81,12 +91,14 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x,
z = precond.solve(s); z = precond.solve(s);
t.noalias() = mat * z; t.noalias() = mat * z;
w = t.dot(s) / t.squaredNorm(); w = t.squaredNorm();
if(w>RealScalar(0))
w = t.dot(s) / t.squaredNorm();
x += alpha * y + w * z; x += alpha * y + w * z;
r = s - w * t; r = s - w * t;
++i; ++i;
} }
tol_error = sqrt(r.squaredNorm()/r0_sqnorm); tol_error = sqrt(r.squaredNorm()/rhs_sqnorm);
iters = i; iters = i;
return true; return true;
} }