diff options
| -rw-r--r-- | Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 11 |
1 files changed, 6 insertions, 5 deletions
diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index e67f09184..be98993f0 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -59,20 +59,21 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, VectorType s(n), t(n); - RealScalar tol2 = tol*tol; + RealScalar tol2 = tol*tol*rhs_sqnorm; RealScalar eps2 = NumTraits<Scalar>::epsilon()*NumTraits<Scalar>::epsilon(); Index i = 0; Index restarts = 0; - while ( r.squaredNorm()/rhs_sqnorm > tol2 && i<maxIters ) + while ( r.squaredNorm() > tol2 && i<maxIters ) { Scalar rho_old = rho; rho = r0.dot(r); 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: + r = rhs - mat * x; r0 = r; rho = r0_sqnorm = r.squaredNorm(); if(restarts++ == 0) @@ -202,8 +203,8 @@ public: template<typename Rhs,typename Dest> void _solve_impl(const MatrixBase<Rhs>& b, Dest& x) const { - // x.setZero(); - x = b; + x.resize(this->rows(),b.cols()); + x.setZero(); _solve_with_guess_impl(b,x); } |