Fix numerical issues with BiCGSTAB.

Eigen/src/IterativeLinearSolvers/BiCGSTAB.h

@@ -31,8 +31,6 @@ namespace internal {
 template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
 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<Scalar, Dynamic, 1> 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 r0_norm = r0.stableNorm();
-  RealScalar rhs_sqnorm = rhs.squaredNorm();
+  RealScalar r_norm = r0_norm;
-  if (rhs_sqnorm == 0) {
+  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 eps = NumTraits<Scalar>::epsilon();
-  RealScalar eps2 = NumTraits<Scalar>::epsilon() * NumTraits<Scalar>::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;
-
+    Scalar theta = r0.dot(v);
-    alpha = rho / 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;
 }

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<double, Eigen::RowMajor>;
+  CSRMatrix A = D.sparseView();
+
+  Eigen::VectorXd b = Eigen::VectorXd::Zero(14);
+  b(12) = -1001;
+
+  Eigen::BiCGSTAB<CSRMatrix> 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<Eigen::MatrixXd> 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<double, int>()));
   CALL_SUBTEST_2((test_bicgstab_T<std::complex<double>, int>()));
   CALL_SUBTEST_3((test_bicgstab_T<double, long int>()));
+  CALL_SUBTEST_4(test_2856());
+  CALL_SUBTEST_5(test_2899());
 }