IDRS.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>
//
// 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/.

#ifndef EIGEN_IDRS_H
#define EIGEN_IDRS_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

namespace internal {
template <typename Vector, typename RealScalar>
typename Vector::Scalar omega(const Vector& t, const Vector& s, RealScalar angle) {
  using numext::abs;
  typedef typename Vector::Scalar Scalar;
  const RealScalar ns = s.stableNorm();
  const RealScalar nt = t.stableNorm();
  const Scalar ts = t.dot(s);
  const RealScalar rho = abs(ts / (nt * ns));

  if (rho < angle) {
    if (ts == Scalar(0)) {
      return Scalar(0);
    }
    // Original relation for om is given by
    // om = om * angle / rho;
    // To alleviate potential (near) division by zero this can be rewritten as
    // om = angle * (ns / nt) * (ts / abs(ts)) = angle * (ns / nt) * sgn(ts)
    return angle * (ns / nt) * (ts / abs(ts));
  }
  return ts / (nt * nt);
}

template <typename MatrixType, typename Rhs, typename Dest, typename Preconditioner>
bool idrs(const MatrixType& A, const Rhs& b, Dest& x, const Preconditioner& precond, Index& iter,
          typename Dest::RealScalar& relres, Index S, bool smoothing, typename Dest::RealScalar angle,
          bool replacement) {
  typedef typename Dest::RealScalar RealScalar;
  typedef typename Dest::Scalar Scalar;
  typedef Matrix<Scalar, Dynamic, 1> VectorType;
  typedef Matrix<Scalar, Dynamic, Dynamic, ColMajor> DenseMatrixType;
  const Index N = b.size();
  S = S < x.rows() ? S : x.rows();
  const RealScalar tol = relres;
  const Index maxit = iter;

  bool trueres = false;

  FullPivLU<DenseMatrixType> lu_solver;

  DenseMatrixType P;
  {
    HouseholderQR<DenseMatrixType> qr(DenseMatrixType::Random(N, S));
    P = (qr.householderQ() * DenseMatrixType::Identity(N, S));
  }

  const RealScalar normb = b.stableNorm();

  if (internal::isApprox(normb, RealScalar(0))) {
    // Solution is the zero vector
    x.setZero();
    iter = 0;
    relres = 0;
    return true;
  }
  // from http://homepage.tudelft.nl/1w5b5/IDRS/manual.pdf
  // A peak in the residual is considered dangerously high if‖ri‖/‖b‖> C(tol/epsilon).
  // With epsilon the relative machine precision. The factor tol/epsilon corresponds
  // to the size of a finite precision number that is so large that the absolute
  // round-off error in this number, when propagated through the process, makes it
  // impossible to achieve the required accuracy. The factor C accounts for the
  // accumulation of round-off errors. This parameter has been set to 10^{-3}.
  // mp is epsilon/C 10^3 * eps is very conservative, so normally no residual
  // replacements will take place. It only happens if things go very wrong. Too many
  // restarts may ruin the convergence.
  const RealScalar mp = RealScalar(1e3) * NumTraits<Scalar>::epsilon();

  // Compute initial residual
  const RealScalar tolb = tol * normb;  // Relative tolerance
  VectorType r = b - A * x;

  VectorType x_s, r_s;

  if (smoothing) {
    x_s = x;
    r_s = r;
  }

  RealScalar normr = r.stableNorm();

  if (normr <= tolb) {
    // Initial guess is a good enough solution
    iter = 0;
    relres = normr / normb;
    return true;
  }

  DenseMatrixType G = DenseMatrixType::Zero(N, S);
  DenseMatrixType U = DenseMatrixType::Zero(N, S);
  DenseMatrixType M = DenseMatrixType::Identity(S, S);
  VectorType t(N), v(N);
  Scalar om = 1.;

  // Main iteration loop, guild G-spaces:
  iter = 0;

  while (normr > tolb && iter < maxit) {
    // New right hand size for small system:
    VectorType f = (r.adjoint() * P).adjoint();

    for (Index k = 0; k < S; ++k) {
      // Solve small system and make v orthogonal to P:
      // c = M(k:s,k:s)\f(k:s);
      lu_solver.compute(M.block(k, k, S - k, S - k));
      VectorType c = lu_solver.solve(f.segment(k, S - k));
      // v = r - G(:,k:s)*c;
      v = r - G.rightCols(S - k) * c;
      // Preconditioning
      v = precond.solve(v);

      // Compute new U(:,k) and G(:,k), G(:,k) is in space G_j
      U.col(k) = U.rightCols(S - k) * c + om * v;
      G.col(k) = A * U.col(k);

      // Bi-Orthogonalise the new basis vectors:
      for (Index i = 0; i < k - 1; ++i) {
        // alpha =  ( P(:,i)'*G(:,k) )/M(i,i);
        Scalar alpha = P.col(i).dot(G.col(k)) / M(i, i);
        G.col(k) = G.col(k) - alpha * G.col(i);
        U.col(k) = U.col(k) - alpha * U.col(i);
      }

      // New column of M = P'*G  (first k-1 entries are zero)
      // M(k:s,k) = (G(:,k)'*P(:,k:s))';
      M.block(k, k, S - k, 1) = (G.col(k).adjoint() * P.rightCols(S - k)).adjoint();

      if (internal::isApprox(M(k, k), Scalar(0))) {
        return false;
      }

      // Make r orthogonal to q_i, i = 0..k-1
      Scalar beta = f(k) / M(k, k);
      r = r - beta * G.col(k);
      x = x + beta * U.col(k);
      normr = r.stableNorm();

      if (replacement && normr > tolb / mp) {
        trueres = true;
      }

      // Smoothing:
      if (smoothing) {
        t = r_s - r;
        // gamma is a Scalar, but the conversion is not allowed
        Scalar gamma = t.dot(r_s) / t.stableNorm();
        r_s = r_s - gamma * t;
        x_s = x_s - gamma * (x_s - x);
        normr = r_s.stableNorm();
      }

      if (normr < tolb || iter == maxit) {
        break;
      }

      // New f = P'*r (first k  components are zero)
      if (k < S - 1) {
        f.segment(k + 1, S - (k + 1)) = f.segment(k + 1, S - (k + 1)) - beta * M.block(k + 1, k, S - (k + 1), 1);
      }
    }  // end for

    if (normr < tolb || iter == maxit) {
      break;
    }

    // Now we have sufficient vectors in G_j to compute residual in G_j+1
    // Note: r is already perpendicular to P so v = r
    // Preconditioning
    v = r;
    v = precond.solve(v);

    // Matrix-vector multiplication:
    t = A * v;

    // Computation of a new omega
    om = internal::omega(t, r, angle);

    if (om == RealScalar(0.0)) {
      return false;
    }

    r = r - om * t;
    x = x + om * v;
    normr = r.stableNorm();

    if (replacement && normr > tolb / mp) {
      trueres = true;
    }

    // Residual replacement?
    if (trueres && normr < normb) {
      r = b - A * x;
      trueres = false;
    }

    // Smoothing:
    if (smoothing) {
      t = r_s - r;
      Scalar gamma = t.dot(r_s) / t.stableNorm();
      r_s = r_s - gamma * t;
      x_s = x_s - gamma * (x_s - x);
      normr = r_s.stableNorm();
    }

    iter++;

  }  // end while

  if (smoothing) {
    x = x_s;
  }
  relres = normr / normb;
  return true;
}

}  // namespace internal

template <typename MatrixType_, typename Preconditioner_ = DiagonalPreconditioner<typename MatrixType_::Scalar> >
class IDRS;

namespace internal {

template <typename MatrixType_, typename Preconditioner_>
struct traits<Eigen::IDRS<MatrixType_, Preconditioner_> > {
  typedef MatrixType_ MatrixType;
  typedef Preconditioner_ Preconditioner;
};

}  // namespace internal

template <typename MatrixType_, typename Preconditioner_>
class IDRS : public IterativeSolverBase<IDRS<MatrixType_, Preconditioner_> > {
 public:
  typedef MatrixType_ MatrixType;
  typedef typename MatrixType::Scalar Scalar;
  typedef typename MatrixType::RealScalar RealScalar;
  typedef Preconditioner_ Preconditioner;

 private:
  typedef IterativeSolverBase<IDRS> Base;
  using Base::m_error;
  using Base::m_info;
  using Base::m_isInitialized;
  using Base::m_iterations;
  using Base::matrix;
  Index m_S;
  bool m_smoothing;
  RealScalar m_angle;
  bool m_residual;

 public:
  IDRS() : m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {}

  template <typename MatrixDerived>
  explicit IDRS(const EigenBase<MatrixDerived>& A)
      : Base(A.derived()), m_S(4), m_smoothing(false), m_angle(RealScalar(0.7)), m_residual(false) {}

  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::idrs(matrix(), b, x, Base::m_preconditioner, m_iterations, m_error, m_S, m_smoothing, m_angle,
                              m_residual);

    m_info = (!ret) ? NumericalIssue : m_error <= Base::m_tolerance ? Success : NoConvergence;
  }

  void setS(Index S) {
    if (S < 1) {
      S = 4;
    }

    m_S = S;
  }

  void setSmoothing(bool smoothing) { m_smoothing = smoothing; }

  void setAngle(RealScalar angle) { m_angle = angle; }

  void setResidualUpdate(bool update) { m_residual = update; }
};

}  // namespace Eigen

#endif /* EIGEN_IDRS_H */