From 809d266b493cb222bab1f434473b0c67175bc380 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Antonio=20S=C3=A1nchez?= Date: Tue, 11 Feb 2025 19:41:59 +0000 Subject: [PATCH] Fix numerical issues with BiCGSTAB. --- Eigen/src/IterativeLinearSolvers/BiCGSTAB.h | 49 +++++++++++------- test/bicgstab.cpp | 55 +++++++++++++++++++++ 2 files changed, 87 insertions(+), 17 deletions(-) diff --git a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h index e3154b497..8fdeb849b 100644 --- a/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h +++ b/Eigen/src/IterativeLinearSolvers/BiCGSTAB.h @@ -31,8 +31,6 @@ namespace internal { template bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, const Preconditioner& precond, Index& iters, typename Dest::RealScalar& tol_error) { - using std::abs; - using std::sqrt; typedef typename Dest::RealScalar RealScalar; typedef typename Dest::Scalar Scalar; typedef Matrix VectorType; @@ -43,14 +41,15 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, const Precondition VectorType r = rhs - mat * x; VectorType r0 = r; - RealScalar r0_sqnorm = r0.squaredNorm(); - RealScalar rhs_sqnorm = rhs.squaredNorm(); - if (rhs_sqnorm == 0) { + RealScalar r0_norm = r0.stableNorm(); + RealScalar r_norm = r0_norm; + RealScalar rhs_norm = rhs.stableNorm(); + if (rhs_norm == 0) { x.setZero(); return true; } Scalar rho(1); - Scalar alpha(1); + Scalar alpha(0); Scalar w(1); VectorType v = VectorType::Zero(n), p = VectorType::Zero(n); @@ -59,21 +58,22 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, const Precondition VectorType s(n), t(n); - RealScalar tol2 = tol * tol * rhs_sqnorm; - RealScalar eps2 = NumTraits::epsilon() * NumTraits::epsilon(); + RealScalar eps = NumTraits::epsilon(); Index i = 0; Index restarts = 0; - while (r.squaredNorm() > tol2 && i < maxIters) { + while (r_norm > tol && i < maxIters) { Scalar rho_old = rho; - rho = r0.dot(r); - if (abs(rho) < eps2 * r0_sqnorm) { + if (Eigen::numext::abs(rho) / Eigen::numext::maxi(r0_norm, r_norm) < eps * Eigen::numext::mini(r0_norm, r_norm)) { // 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(); + rho = r.squaredNorm(); + r0_norm = r.stableNorm(); + alpha = Scalar(0); + w = Scalar(1); if (restarts++ == 0) i = 0; } Scalar beta = (rho / rho_old) * (alpha / w); @@ -82,23 +82,38 @@ bool bicgstab(const MatrixType& mat, const Rhs& rhs, Dest& x, const Precondition y = precond.solve(p); v.noalias() = mat * y; - - alpha = rho / r0.dot(v); + Scalar theta = r0.dot(v); + // For small angles ∠(r0, v) < eps, random restart. + RealScalar v_norm = v.stableNorm(); + if (Eigen::numext::abs(theta) / Eigen::numext::maxi(r0_norm, v_norm) < eps * Eigen::numext::mini(r0_norm, v_norm)) { + r = rhs - mat * x; + r0.setRandom(); + r0_norm = r0.stableNorm(); + rho = Scalar(1); + alpha = Scalar(0); + w = Scalar(1); + if (restarts++ == 0) i = 0; + continue; + } + alpha = rho / theta; s = r - alpha * v; z = precond.solve(s); t.noalias() = mat * z; RealScalar tmp = t.squaredNorm(); - if (tmp > RealScalar(0)) + if (tmp > RealScalar(0)) { w = t.dot(s) / tmp; - else + } else { w = Scalar(0); + } x += alpha * y + w * z; r = s - w * t; + r_norm = r.stableNorm(); ++i; } - tol_error = sqrt(r.squaredNorm() / rhs_sqnorm); + + tol_error = r_norm / rhs_norm; iters = i; return true; } diff --git a/test/bicgstab.cpp b/test/bicgstab.cpp index 3f53e3ee9..7ff2f3dee 100644 --- a/test/bicgstab.cpp +++ b/test/bicgstab.cpp @@ -26,8 +26,63 @@ void test_bicgstab_T() { // CALL_SUBTEST( check_sparse_square_solving(bicgstab_colmajor_ssor) ); } +// https://gitlab.com/libeigen/eigen/-/issues/2856 +void test_2856() { + Eigen::MatrixXd D = Eigen::MatrixXd::Identity(14, 14); + D(6, 13) = 1; + D(13, 12) = 1; + using CSRMatrix = Eigen::SparseMatrix; + CSRMatrix A = D.sparseView(); + + Eigen::VectorXd b = Eigen::VectorXd::Zero(14); + b(12) = -1001; + + Eigen::BiCGSTAB solver; + solver.compute(A); + Eigen::VectorXd x = solver.solve(b); + Eigen::VectorXd expected = Eigen::VectorXd::Zero(14); + expected(6) = -1001; + expected(12) = -1001; + expected(13) = 1001; + VERIFY_IS_EQUAL(x, expected); + + Eigen::VectorXd residual = b - A * x; + VERIFY(residual.isZero()); +} + +// https://gitlab.com/libeigen/eigen/-/issues/2899 +void test_2899() { + Eigen::MatrixXd A(4, 4); + A(0, 0) = 1; + A(1, 0) = -1.0 / 6; + A(1, 1) = 2.0 / 3; + A(1, 2) = -1.0 / 6; + A(1, 3) = -1.0 / 3; + A(2, 1) = -1.0 / 3; + A(2, 2) = 1; + A(2, 3) = -2.0 / 3; + A(3, 1) = -1.0 / 3; + A(3, 2) = -1.0 / 3; + A(3, 3) = 2.0 / 3; + Eigen::VectorXd b(4); + b(0) = 0; + b(1) = 1; + b(2) = 1; + b(3) = 1; + Eigen::BiCGSTAB solver; + solver.compute(A); + Eigen::VectorXd x = solver.solve(b); + Eigen::VectorXd expected(4); + expected << 0, 15, 18, 18; + VERIFY_IS_APPROX(x, expected); + Eigen::VectorXd residual = b - A * x; + VERIFY(residual.isZero()); +} + EIGEN_DECLARE_TEST(bicgstab) { CALL_SUBTEST_1((test_bicgstab_T())); CALL_SUBTEST_2((test_bicgstab_T, int>())); CALL_SUBTEST_3((test_bicgstab_T())); + CALL_SUBTEST_4(test_2856()); + CALL_SUBTEST_5(test_2899()); }