BiCGSTABL.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2020 Chris Schoutrop <c.e.m.schoutrop@tue.nl>
// Copyright (C) 2020 Jens Wehner <j.wehner@esciencecenter.nl>
// Copyright (C) 2020 Jan van Dijk <j.v.dijk@tue.nl>
// Copyright (C) 2020 Adithya Vijaykumar
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
/*
 
  This implementation of BiCGStab(L) is based on the papers
      General algorithm:
      1. G.L.G. Sleijpen, D.R. Fokkema. (1993). BiCGstab(l) for linear equations
  involving unsymmetric matrices with complex spectrum. Electronic Transactions
  on Numerical Analysis. Polynomial step update:
      2. G.L.G. Sleijpen, M.B. Van Gijzen. (2010) Exploiting BiCGstab(l)
  strategies to induce dimension reduction SIAM Journal on Scientific Computing.
      3. Fokkema, Diederik R. Enhanced implementation of BiCGstab (l) for
  solving linear systems of equations. Universiteit Utrecht. Mathematisch
  Instituut, 1996
      4. Sleijpen, G. L., & van der Vorst, H. A. (1996). Reliable updated
  residuals in hybrid Bi-CG methods. Computing, 56(2), 141-163.
*/
 
#ifndef EIGEN_BICGSTABL_H
#define EIGEN_BICGSTABL_H
 
namespace Eigen {
 
namespace internal {
template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
bool bicgstabl(const MatrixType &mat, const Rhs &rhs, Dest &x, const Preconditioner &precond, Index &iters,
               typename Dest::RealScalar &tol_error, Index L) {
  using numext::abs;
  using numext::sqrt;
  typedef typename Dest::RealScalar RealScalar;
  typedef typename Dest::Scalar Scalar;
  const Index N = rhs.size();
  L = L < x.rows() ? L : x.rows();
 
  Index k = 0;
 
  const RealScalar tol = tol_error;
  const Index maxIters = iters;
 
  typedef Matrix<Scalar, Dynamic, 1> VectorType;
  typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> DenseMatrixType;
 
  DenseMatrixType rHat(N, L + 1);
  DenseMatrixType uHat(N, L + 1);
 
  // We start with an initial guess x_0 and let us set r_0 as (residual
  // calculated from x_0)
  VectorType x0 = x;
  rHat.col(0) = rhs - mat * x0;  // r_0
 
  x.setZero();  // This will contain the updates to the solution.
  // rShadow is arbitrary, but must never be orthogonal to any residual.
  VectorType rShadow = VectorType::Random(N);
 
  VectorType x_prime = x;
 
  // Redundant: x is already set to 0
  // x.setZero();
  VectorType b_prime = rHat.col(0);
 
  // Other vectors and scalars initialization
  Scalar rho0 = 1.0;
  Scalar alpha = 0.0;
  Scalar omega = 1.0;
 
  uHat.col(0).setZero();
 
  bool bicg_convergence = false;
 
  const RealScalar normb = rhs.stableNorm();
  if (internal::isApprox(normb, RealScalar(0))) {
    x.setZero();
    iters = 0;
    return true;
  }
  RealScalar normr = rHat.col(0).stableNorm();
  RealScalar Mx = normr;
  RealScalar Mr = normr;
 
  // Keep track of the solution with the lowest residual
  RealScalar normr_min = normr;
  VectorType x_min = x_prime + x;
 
  // Criterion for when to apply the group-wise update, conform ref 3.
  const RealScalar delta = 0.01;
 
  bool compute_res = false;
  bool update_app = false;
 
  while (normr > tol * normb && k < maxIters) {
    rho0 *= -omega;
 
    for (Index j = 0; j < L; ++j) {
      const Scalar rho1 = rShadow.dot(rHat.col(j));
 
      if (!(numext::isfinite)(rho1) || rho0 == RealScalar(0.0)) {
        // We cannot continue computing, return the best solution found.
        x += x_prime;
 
        // Check if x is better than the best stored solution thus far.
        normr = (rhs - mat * (precond.solve(x) + x0)).stableNorm();
 
        if (normr > normr_min || !(numext::isfinite)(normr)) {
          // x_min is a better solution than x, return x_min
          x = x_min;
          normr = normr_min;
        }
        tol_error = normr / normb;
        iters = k;
        // x contains the updates to x0, add those back to obtain the solution
        x = precond.solve(x);
        x += x0;
        return (normr < tol * normb);
      }
 
      const Scalar beta = alpha * (rho1 / rho0);
      rho0 = rho1;
      // Update search directions
      uHat.leftCols(j + 1) = rHat.leftCols(j + 1) - beta * uHat.leftCols(j + 1);
      uHat.col(j + 1) = mat * precond.solve(uHat.col(j));
      const Scalar sigma = rShadow.dot(uHat.col(j + 1));
      alpha = rho1 / sigma;
      // Update residuals
      rHat.leftCols(j + 1) -= alpha * uHat.middleCols(1, j + 1);
      rHat.col(j + 1) = mat * precond.solve(rHat.col(j));
      // Complete BiCG iteration by updating x
      x += alpha * uHat.col(0);
      normr = rHat.col(0).stableNorm();
      // Check for early exit
      if (normr < tol * normb) {
        /*
          Convergence was achieved during BiCG step.
          Without this check BiCGStab(L) fails for trivial matrices, such as
          when the preconditioner already is the inverse, or the input matrix is
          identity.
        */
        bicg_convergence = true;
        break;
      } else if (normr < normr_min) {
        // We found an x with lower residual, keep this one.
        x_min = x + x_prime;
        normr_min = normr;
      }
    }
    if (!bicg_convergence) {
      /*
        The polynomial/minimize residual step.
 
        QR Householder method for argmin is more stable than (modified)
        Gram-Schmidt, in the sense that there is less loss of orthogonality. It
        is more accurate than solving the normal equations, since the normal
        equations scale with condition number squared.
      */
      const VectorType gamma = rHat.rightCols(L).householderQr().solve(rHat.col(0));
      x += rHat.leftCols(L) * gamma;
      rHat.col(0) -= rHat.rightCols(L) * gamma;
      uHat.col(0) -= uHat.rightCols(L) * gamma;
      normr = rHat.col(0).stableNorm();
      omega = gamma(L - 1);
    }
    if (normr < normr_min) {
      // We found an x with lower residual, keep this one.
      x_min = x + x_prime;
      normr_min = normr;
    }
 
    k++;
 
    /*
      Reliable update part
 
      The recursively computed residual can deviate from the actual residual
      after several iterations. However, computing the residual from the
      definition costs extra MVs and should not be done at each iteration. The
      reliable update strategy computes the true residual from the definition:
      r=b-A*x at strategic intervals. Furthermore a "group wise update" strategy
      is used to combine updates, which improves accuracy.
    */
 
    // Maximum norm of residuals since last update of x.
    Mx = numext::maxi(Mx, normr);
    // Maximum norm of residuals since last computation of the true residual.
    Mr = numext::maxi(Mr, normr);
 
    if (normr < delta * normb && normb <= Mx) {
      update_app = true;
    }
 
    if (update_app || (normr < delta * Mr && normb <= Mr)) {
      compute_res = true;
    }
 
    if (bicg_convergence) {
      update_app = true;
      compute_res = true;
      bicg_convergence = false;
    }
 
    if (compute_res) {
      // Explicitly compute residual from the definition
 
      // This is equivalent to the shifted version of rhs - mat *
      // (precond.solve(x)+x0)
      rHat.col(0) = b_prime - mat * precond.solve(x);
      normr = rHat.col(0).stableNorm();
      Mr = normr;
 
      if (update_app) {
        // After the group wise update, the original problem is translated to a
        // shifted one.
        x_prime += x;
        x.setZero();
        b_prime = rHat.col(0);
        Mx = normr;
      }
    }
    if (normr < normr_min) {
      // We found an x with lower residual, keep this one.
      x_min = x + x_prime;
      normr_min = normr;
    }
 
    compute_res = false;
    update_app = false;
  }
 
  // Convert internal variable to the true solution vector x
  x += x_prime;
 
  normr = (rhs - mat * (precond.solve(x) + x0)).stableNorm();
  if (normr > normr_min || !(numext::isfinite)(normr)) {
    // x_min is a better solution than x, return x_min
    x = x_min;
    normr = normr_min;
  }
  tol_error = normr / normb;
  iters = k;
 
  // x contains the updates to x0, add those back to obtain the solution
  x = precond.solve(x);
  x += x0;
  return true;
}
 
}  // namespace internal
 
template <typename MatrixType_, typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar>>
class BiCGSTABL;
 
namespace internal {
 
template <typename MatrixType_, typename Preconditioner_>
struct traits<Eigen::BiCGSTABL<MatrixType_, Preconditioner_>> {
  typedef MatrixType_ MatrixType;
  typedef Preconditioner_ Preconditioner;
};
 
}  // namespace internal
 
template <typename MatrixType_, typename Preconditioner_>
class BiCGSTABL : public IterativeSolverBase<BiCGSTABL<MatrixType_, Preconditioner_>> {
  typedef IterativeSolverBase<BiCGSTABL> Base;
  using Base::m_error;
  using Base::m_info;
  using Base::m_isInitialized;
  using Base::m_iterations;
  using Base::matrix;
  Index m_L;
 
 public:
  typedef MatrixType_ MatrixType;
  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::RealScalar RealScalar;
  typedef Preconditioner_ Preconditioner;

  BiCGSTABL() : m_L(2) {}

  template <typename MatrixDerived>
  explicit BiCGSTABL(const EigenBase<MatrixDerived> &A) : Base(A.derived()), m_L(2) {}

  template <typename Rhs, typename Dest>
  void _solve_vector_with_guess_impl(const Rhs &b, Dest &x) const {
    m_iterations = Base::maxIterations();
 
    m_error = Base::m_tolerance;
 
    bool ret = internal::bicgstabl(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, m_L);
    m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence;
  }

  void setL(Index L) {
    eigen_assert(L >= 1 && "L needs to be positive");
    m_L = L;
  }
};
 
}  // namespace Eigen
 
#endif /* EIGEN_BICGSTABL_H */