smoother.hpp

Go to the documentation of this file.

#pragma once
 
#include <string>
 
#include <Eigen/Core>
#include <Eigen/Sparse>
 
#include <amg/common.hpp>
 
namespace AMG {
 
template <class EleType>
class SmootherBase {
 public:
  EleType tolerance{1e-9};
 
  size_t compute_error_every_n_iters{100};
 
  size_t n_iters{1};
 
  SmootherBase() {}
 
  SmootherBase(size_t n_iters_) : n_iters(n_iters_) {}
 
  SmootherBase(double tolerance_, size_t compute_error_every_n_iters_,
               size_t n_iters_)
      : tolerance(tolerance_),
        compute_error_every_n_iters(compute_error_every_n_iters_),
        n_iters(n_iters_) {}
 
  virtual void smooth(const Eigen::SparseMatrix<EleType>& A,
                      Eigen::Matrix<EleType, -1, 1>& u,
                      const Eigen::Matrix<EleType, -1, 1>& b) = 0;
};
 
template <class EleType>
class SparseGaussSeidel : public SmootherBase<EleType> {
 private:
  void matvecprod(int col, EleType& rsum, EleType& diag, const EleType& z,
                  const Eigen::SparseMatrix<EleType>& A,
                  const Eigen::Matrix<EleType, -1, 1>& u) {
    int row;
    EleType val;
    for (typename Eigen::SparseMatrix<EleType>::InnerIterator it(A, col); it;
         ++it) {
      row = it.row();
      val = it.value();
      diag = (col == row) // if you found a value on the diagonal, update it for other iters
                 ? val
                 : diag;  // column == row therefore val == A_ii on diag
      rsum += (col == row)
                  ? z
                  : val * u[row];  // contrib to sum should be 0 for i==j by def
    }
  }
 
  void update(int col, const Eigen::SparseMatrix<EleType>& A,
              const Eigen::Matrix<EleType, -1, 1>& b,
              Eigen::Matrix<EleType, -1, 1>& u) {
    EleType z = 0;
    EleType rsum = z;
    EleType diag = z;
    matvecprod(col, rsum, diag, z, A, u);
    u[col] = diag == z ? u[col] : (b[col] - rsum) / diag;
    return;
  }
 
  void forwardsweep(const Eigen::SparseMatrix<EleType>& A,
                    const Eigen::Matrix<EleType, -1, 1>& b,
                    Eigen::Matrix<EleType, -1, 1>& u, const int& ncols) {
 
    // iterate through cols of A in forward direction
    for (int col = 0; col < ncols; ++col) {
      update(col, A, b, u);
    }
    return;
  }
 
  void backwardsweep(const Eigen::SparseMatrix<EleType>& A,
                     const Eigen::Matrix<EleType, -1, 1>& b,
                     Eigen::Matrix<EleType, -1, 1>& u, const int& ncols) {
    // iterate through cols A in the backward direction
    for (int col = ncols - 1; col >= 0; --col) {
      update(col, A, b, u);
    }
  }
 
 public:
  using SmootherBase<EleType>::SmootherBase;
 
  SparseGaussSeidel() {
    this->tolerance = 1e-9;
    this->compute_error_every_n_iters = 0;
    this->n_iters = 1;
  }
 
  void smooth(const Eigen::SparseMatrix<EleType>& A,
              Eigen::Matrix<EleType, -1, 1>& u,
              const Eigen::Matrix<EleType, -1, 1>& b) {
    int ncols = A.cols();
    size_t iter = 0;
    EleType error = 100;
    while (iter < this->n_iters && error > this->tolerance) {
      forwardsweep(A, b, u, ncols);
      backwardsweep(A, b, u, ncols);
      iter += 1;
      if (this->compute_error_every_n_iters != 0 &&
          iter % this->compute_error_every_n_iters == 0) {
        error = rss(A, u, b);
      }
    }
 
    if (this->compute_error_every_n_iters != 0) {
      if (error <= this->tolerance)
        std::cout << "SPGS converged after " << iter << " iterations."
                  << std::endl;
      else
        std::cout << "SPGS did not converge after " << iter << " iterations."
                  << std::endl;
    }
 
    return;
  }
};
 
template <class EleType>
class Jacobi : public SmootherBase<EleType> {
 public:
  // C++11
  // https://stackoverflow.com/questions/8093882/using-c-base-class-constructors
  using SmootherBase<EleType>::SmootherBase;
 
  Jacobi() {}
 
  void smooth(const Eigen::SparseMatrix<EleType>& A,
              Eigen::Matrix<EleType, -1, 1>& u,
              const Eigen::Matrix<EleType, -1, 1>& b) {
    size_t iter = 0;
    size_t ndofs = b.size();
    EleType error = 100;
    EleType sigma;
    while (iter < this->n_iters && error > this->tolerance) {
      for (size_t i = 0; i < ndofs; ++i) {
        sigma = 0;
        for (size_t j = 0; j < ndofs; ++j) {
          if (j != i) {
            sigma += A.coeff(i, j) * u[j];
          }
        }
        EleType aii = A.coeff(i, i);
        u[i] = (b[i] - sigma) / aii;
      }
      iter += 1;
      if (iter % this->compute_error_every_n_iters == 0) {
        error = rss(A, u, b);
      }
    }
    return;
  }
};
 
template <class EleType>
class SuccessiveOverRelaxation : public SmootherBase<EleType> {
 private:
  double omega{1.0};
 
  void validate_omega() {
    if (omega > 2 || omega < 0) {
      std::string msg =
          "`omega` must be in [0, 2] but got omega=" + std::to_string(omega) +
          "\n";
      throw std::invalid_argument(msg);
    }
  }
 
 public:
  // C++11
  // https://stackoverflow.com/questions/8093882/using-c-base-class-constructors
  using SmootherBase<EleType>::SmootherBase;
 
  SuccessiveOverRelaxation() {}
  SuccessiveOverRelaxation(double omega_) : omega(omega_) { validate_omega(); }
 
  SuccessiveOverRelaxation(double omega_, double tolerance_,
                           size_t compute_error_every_n_iters_, size_t n_iters_)
      : SmootherBase<EleType>(tolerance_, compute_error_every_n_iters_,
                              n_iters_),
        omega(omega_) {
    validate_omega();
  }
 
  void smooth(const Eigen::SparseMatrix<EleType>& A,
              Eigen::Matrix<EleType, -1, 1>& u,
              const Eigen::Matrix<EleType, -1, 1>& b) {
    size_t iter = 0;
    size_t ndofs = b.size();
    EleType error = 100;
    EleType sigma_j_less_i;
    EleType sigma_j_greater_i;
    while (iter < this->n_iters && error > this->tolerance) {
      for (size_t i = 0; i < ndofs; ++i) {
        sigma_j_less_i = 0;
        for (size_t j = 0; j < i; ++j) {
          sigma_j_less_i += A.coeff(i, j) * u[j];
        }
 
        sigma_j_greater_i = 0;
        for (size_t j = i + 1; j < ndofs; ++j) {
          sigma_j_greater_i += A.coeff(i, j) * u[j];
        }
 
        EleType aii = A.coeff(i, i);
        EleType uk_plus_one_gauss_seidel =
            (b[i] - sigma_j_less_i - sigma_j_greater_i) / aii;
        EleType uk = u[i];
 
        u[i] = uk + omega * (uk_plus_one_gauss_seidel - uk);
      }
      iter += 1;
      if (iter % this->compute_error_every_n_iters == 0) {
        error = rss(A, u, b);
      }
    }
    return;
  }
};
 
}  // end namespace AMG