AMG.h

// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
//-----------------------------------------------------------------------------
// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: Apache-2.0
//-----------------------------------------------------------------------------
/*---------------------------------------------------------------------------*/
/* AMG.h                                                       (C) 2000-2026 */
/*                                                                           */
/* An AMG preconditioner.                                     .              */
/*---------------------------------------------------------------------------*/
#ifndef ARCCORE_ALINA_AMG_H
#define ARCCORE_ALINA_AMG_H
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
/*
 * This file is based on the work on AMGCL library (version march 2026)
 * which can be found at https://github.com/ddemidov/amgcl.
 *
 * Copyright (c) 2012-2022 Denis Demidov <dennis.demidov@gmail.com>
 * SPDX-License-Identifier: MIT
 */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
// Temporarily add this pragma because there is a lot of conversion warnings.
#pragma GCC diagnostic ignored "-Wconversion"
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
#include "arccore/alina/BuiltinBackend.h"
#include "arccore/alina/SolverUtils.h"
#include "arccore/alina/AlinaUtils.h"
 
#include <iostream>
#include <iomanip>
#include <list>
#include <memory>
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
namespace Arcane::Alina
{
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
template <class Backend,
          template <class> class Coarsening,
          template <class> class Relax>
class AMG
{
 public:
 
  typedef Backend backend_type;
 
  typedef typename Backend::value_type value_type;
  typedef typename Backend::col_type col_type;
  typedef typename Backend::ptr_type ptr_type;
  typedef typename Backend::matrix matrix;
  typedef typename Backend::vector vector;
 
  typedef Coarsening<Backend> coarsening_type;
  typedef Relax<Backend> relax_type;
 
  typedef typename BuiltinBackend<value_type, col_type, ptr_type>::matrix build_matrix;
 
  typedef typename math::scalar_of<value_type>::type scalar_type;

  typedef typename Backend::params backend_params;

  struct params
  {
    typedef coarsening_type::params coarsening_params;
    typedef relax_type::params relax_params;
 
    coarsening_params coarsening; 
    relax_params relax; 

    Int32 coarse_enough = Backend::direct_solver::coarse_enough();

    bool direct_coarse = true;

    Int32 max_levels = std::numeric_limits<Int32>::max();

    Int32 npre = 1;

    Int32 npost = 1;

    Int32 ncycle = 1;

    Int32 pre_cycles = 1;

    bool allow_rebuild = std::is_same<matrix, build_matrix>::value;
 
    params() = default;
 
    params(const PropertyTree& p)
    : ARCCORE_ALINA_PARAMS_IMPORT_CHILD(p, coarsening)
    , ARCCORE_ALINA_PARAMS_IMPORT_CHILD(p, relax)
    , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, coarse_enough)
    , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, direct_coarse)
    , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, max_levels)
    , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, npre)
    , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, npost)
    , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, ncycle)
    , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, pre_cycles)
    , ARCCORE_ALINA_PARAMS_IMPORT_VALUE(p, allow_rebuild)
    {
      p.check_params( { "coarsening", "relax", "coarse_enough", "direct_coarse", "max_levels", "npre", "npost", "ncycle", "pre_cycles", "allow_rebuild" });
 
      precondition(max_levels > 0, "max_levels should be positive");
    }
 
    void get(PropertyTree& p, const std::string& path = "") const
    {
      ARCCORE_ALINA_PARAMS_EXPORT_CHILD(p, path, coarsening);
      ARCCORE_ALINA_PARAMS_EXPORT_CHILD(p, path, relax);
      ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, coarse_enough);
      ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, direct_coarse);
      ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, max_levels);
      ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, npre);
      ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, npost);
      ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, ncycle);
      ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, pre_cycles);
      ARCCORE_ALINA_PARAMS_EXPORT_VALUE(p, path, allow_rebuild);
    }
  };

  template <class Matrix>
  explicit AMG([[maybe_unused]] const AlinaContext& context,
               const Matrix& M, const params& p = params(),
               const backend_params& bprm = backend_params())
  : prm(p)
  {
    auto A = std::make_shared<build_matrix>(M);
    sort_rows(*A);
 
    _initialize(A, bprm);
  }

  explicit AMG([[maybe_unused]] const AlinaContext& context,
               std::shared_ptr<build_matrix> A,
               const params& p = params(),
               const backend_params& bprm = backend_params())
  : prm(p)
  {
    _initialize(A, bprm);
  }

  template <class Matrix>
  void rebuild(const Matrix& M,
               const backend_params& bprm = backend_params())
  {
    auto A = std::make_shared<build_matrix>(M);
    sort_rows(*A);
    rebuild(A, bprm);
  }

  void rebuild(std::shared_ptr<build_matrix> A,
               const backend_params& bprm = backend_params())
  {
    precondition(prm.allow_rebuild, "allow_rebuild is not set!");
    precondition(backend::nbRow(*A) == backend::nbRow(system_matrix()) &&
                 backend::nbColumn(*A) == backend::nbRow(*A),
                 "Matrix dimensions differ from the original ones!");
 
    ARCCORE_ALINA_TIC("rebuild");
    coarsening_type C(prm.coarsening);
    for (auto& level : levels) {
      A = level.rebuild(A, C, prm, bprm);
    }
    ARCCORE_ALINA_TOC("rebuild");
  }

  template <class Vec1, class Vec2>
  void cycle(const Vec1& rhs, Vec2&& x) const
  {
    cycle(levels.begin(), rhs, x);
  }

  template <class Vec1, class Vec2>
  void apply(const Vec1& rhs, Vec2&& x) const
  {
    if (prm.pre_cycles) {
      backend::clear(x);
      for (unsigned i = 0; i < prm.pre_cycles; ++i)
        cycle(rhs, x);
    }
    else {
      backend::copy(rhs, x);
    }
  }

  std::shared_ptr<matrix> system_matrix_ptr() const
  {
    return levels.front().A;
  }
 
  const matrix& system_matrix() const
  {
    return *system_matrix_ptr();
  }
 
  size_t bytes() const
  {
    size_t b = 0;
    for (const auto& lvl : levels)
      b += lvl.bytes();
    return b;
  }
 
 public:
 
  params prm;
 
 private:
 
  struct AMGLevel
  {
    size_t m_rows = 0;
    size_t m_nonzeros = 0;
 
    std::shared_ptr<vector> f;
    std::shared_ptr<vector> u;
    std::shared_ptr<vector> t;
 
    std::shared_ptr<matrix> A;
    std::shared_ptr<matrix> P;
    std::shared_ptr<matrix> R;
 
    std::shared_ptr<build_matrix> bP;
    std::shared_ptr<build_matrix> bR;
 
    std::shared_ptr<typename Backend::direct_solver> solve;
 
    std::shared_ptr<relax_type> relax;
 
    size_t bytes() const
    {
      size_t b = 0;
 
      if (f)
        b += backend::bytes(*f);
      if (u)
        b += backend::bytes(*u);
      if (t)
        b += backend::bytes(*t);
 
      if (A)
        b += backend::bytes(*A);
      if (P)
        b += backend::bytes(*P);
      if (R)
        b += backend::bytes(*R);
 
      if (solve)
        b += backend::bytes(*solve);
      if (relax)
        b += backend::bytes(*relax);
 
      return b;
    }
 
    AMGLevel() = default;
 
    AMGLevel(std::shared_ptr<build_matrix> A, params& prm, const backend_params& bprm)
    : m_rows(backend::nbRow(*A))
    , m_nonzeros(backend::nonzeros(*A))
    {
      ARCCORE_ALINA_TIC("move to backend");
      f = Backend::create_vector(m_rows, bprm);
      u = Backend::create_vector(m_rows, bprm);
      t = Backend::create_vector(m_rows, bprm);
      this->A = Backend::copy_matrix(A, bprm);
      ARCCORE_ALINA_TOC("move to backend");
 
      ARCCORE_ALINA_TIC("relaxation");
      relax = std::make_shared<relax_type>(*A, prm.relax, bprm);
      ARCCORE_ALINA_TOC("relaxation");
    }
 
    std::shared_ptr<build_matrix>
    step_down(std::shared_ptr<build_matrix> A, coarsening_type& C,
              const backend_params& bprm, bool allow_rebuild)
    {
      ARCCORE_ALINA_TIC("transfer operators");
      std::shared_ptr<build_matrix> P, R;
 
      try {
        std::tie(P, R) = C.transfer_operators(*A);
      }
      catch (error::empty_level) {
        ARCCORE_ALINA_TOC("transfer operators");
        return std::shared_ptr<build_matrix>();
      }
 
      sort_rows(*P);
      sort_rows(*R);
 
      if (allow_rebuild) {
        bP = P;
        bR = R;
      }
      ARCCORE_ALINA_TOC("transfer operators");
 
      ARCCORE_ALINA_TIC("move to backend");
      this->P = Backend::copy_matrix(P, bprm);
      this->R = Backend::copy_matrix(R, bprm);
      ARCCORE_ALINA_TOC("move to backend");
 
      ARCCORE_ALINA_TIC("coarse operator");
      A = C.coarse_operator(*A, *P, *R);
      sort_rows(*A);
      ARCCORE_ALINA_TOC("coarse operator");
 
      return A;
    }
 
    void create_coarse(std::shared_ptr<build_matrix> A,
                       const backend_params& bprm, bool single_level)
    {
      m_rows = backend::nbRow(*A);
      m_nonzeros = backend::nonzeros(*A);
 
      u = Backend::create_vector(m_rows, bprm);
      f = Backend::create_vector(m_rows, bprm);
 
      solve = Backend::create_solver(A, bprm);
      if (single_level)
        this->A = Backend::copy_matrix(A, bprm);
    }
 
    std::shared_ptr<build_matrix>
    rebuild(std::shared_ptr<build_matrix> A,
            const coarsening_type& C,
            const params& prm,
            const backend_params& bprm)
    {
      if (this->A) {
        ARCCORE_ALINA_TIC("move to backend");
        this->A = Backend::copy_matrix(A, bprm);
        ARCCORE_ALINA_TOC("move to backend");
      }
 
      if (relax) {
        ARCCORE_ALINA_TIC("relaxation");
        relax = std::make_shared<relax_type>(*A, prm.relax, bprm);
        ARCCORE_ALINA_TOC("relaxation");
      }
 
      if (solve) {
        ARCCORE_ALINA_TIC("coarsest level");
        solve = Backend::create_solver(A, bprm);
        ARCCORE_ALINA_TOC("coarsest level");
      }
 
      if (bP && bR) {
        ARCCORE_ALINA_TIC("coarse operator");
        A = C.coarse_operator(*A, *bP, *bR);
        sort_rows(*A);
        ARCCORE_ALINA_TOC("coarse operator");
      }
 
      return A;
    }
 
    size_t rows() const
    {
      return m_rows;
    }
 
    size_t nonzeros() const
    {
      return m_nonzeros;
    }
  };
 
  typedef std::list<AMGLevel>::const_iterator level_iterator;
 
  std::list<AMGLevel> levels;
 
  void _initialize(std::shared_ptr<build_matrix> A,
                   const backend_params& bprm = backend_params())
  {
    precondition(backend::nbRow(*A) == backend::nbColumn(*A), "Matrix should be square!");
 
    bool direct_coarse_solve = true;
 
    coarsening_type C(prm.coarsening);
    std::cout << "DoInit AMG coarse_enough=" << prm.coarse_enough << "\n";
    while (backend::nbRow(*A) > prm.coarse_enough) {
      std::cout << "DoIteration nb_row=" << backend::nbRow(*A) << "\n";
      levels.push_back(AMGLevel(A, prm, bprm));
 
      if (levels.size() >= prm.max_levels)
        break;
 
      A = levels.back().step_down(A, C, bprm, prm.allow_rebuild);
      if (!A) {
        // Zero-sized coarse level. Probably the system matrix on
        // this level is diagonal, should be easily solvable with a
        // couple of smoother iterations.
        direct_coarse_solve = false;
        break;
      }
    }
 
    if (!A || backend::nbRow(*A) > prm.coarse_enough) {
      // The coarse matrix is still too big to be solved directly.
      direct_coarse_solve = false;
    }
 
    if (direct_coarse_solve) {
      ARCCORE_ALINA_TIC("coarsest level");
      if (prm.direct_coarse) {
        AMGLevel l;
        l.create_coarse(A, bprm, levels.empty());
        levels.push_back(l);
      }
      else {
        levels.push_back(AMGLevel(A, prm, bprm));
      }
      ARCCORE_ALINA_TOC("coarsest level");
    }
  }
 
  template <class Vec1, class Vec2>
  void cycle(level_iterator lvl, const Vec1& rhs, Vec2& x) const
  {
    level_iterator nxt = lvl, end = levels.end();
    ++nxt;
 
    if (nxt == end) {
      if (lvl->solve) {
        ARCCORE_ALINA_TIC("coarse");
        (*lvl->solve)(rhs, x);
        ARCCORE_ALINA_TOC("coarse");
      }
      else {
        ARCCORE_ALINA_TIC("relax");
        for (size_t i = 0; i < prm.npre; ++i)
          lvl->relax->apply_pre(*lvl->A, rhs, x, *lvl->t);
        for (size_t i = 0; i < prm.npost; ++i)
          lvl->relax->apply_post(*lvl->A, rhs, x, *lvl->t);
        ARCCORE_ALINA_TOC("relax");
      }
    }
    else {
      for (size_t j = 0; j < prm.ncycle; ++j) {
        ARCCORE_ALINA_TIC("relax");
        for (size_t i = 0; i < prm.npre; ++i)
          lvl->relax->apply_pre(*lvl->A, rhs, x, *lvl->t);
        ARCCORE_ALINA_TOC("relax");
 
        backend::residual(rhs, *lvl->A, x, *lvl->t);
 
        backend::spmv(math::identity<scalar_type>(), *lvl->R, *lvl->t, math::zero<scalar_type>(), *nxt->f);
 
        backend::clear(*nxt->u);
        cycle(nxt, *nxt->f, *nxt->u);
 
        backend::spmv(math::identity<scalar_type>(), *lvl->P, *nxt->u, math::identity<scalar_type>(), x);
 
        ARCCORE_ALINA_TIC("relax");
        for (size_t i = 0; i < prm.npost; ++i)
          lvl->relax->apply_post(*lvl->A, rhs, x, *lvl->t);
        ARCCORE_ALINA_TOC("relax");
      }
    }
  }
 
  template <class B, template <class> class C, template <class> class R>
  friend std::ostream& operator<<(std::ostream& os, const AMG<B, C, R>& a);
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

template <class B, template <class> class C, template <class> class R>
std::ostream& operator<<(std::ostream& os, const AMG<B, C, R>& a)
{
  typedef typename AMG<B, C, R>::AMGLevel level;
  ScopedStreamModifier ss(os);
 
  size_t sum_dof = 0;
  size_t sum_nnz = 0;
  size_t sum_mem = 0;
 
  for (const level& lvl : a.levels) {
    sum_dof += lvl.rows();
    sum_nnz += lvl.nonzeros();
    sum_mem += lvl.bytes();
  }
 
  os << "Number of levels:    " << a.levels.size()
     << "\nOperator complexity: " << std::fixed << std::setprecision(2)
     << 1.0 * sum_nnz / a.levels.front().nonzeros()
     << "\nGrid complexity:     " << std::fixed << std::setprecision(2)
     << 1.0 * sum_dof / a.levels.front().rows()
     << "\nMemory footprint:    " << human_readable_memory(sum_mem)
     << "\n\n"
        "level     unknowns       nonzeros      memory\n"
        "---------------------------------------------\n";
 
  size_t depth = 0;
  for (const level& lvl : a.levels) {
    os << std::setw(5) << depth++
       << std::setw(13) << lvl.rows()
       << std::setw(15) << lvl.nonzeros()
       << std::setw(12) << human_readable_memory(lvl.bytes())
       << " (" << std::setw(5) << std::fixed << std::setprecision(2)
       << 100.0 * lvl.nonzeros() / sum_nnz
       << "%)" << std::endl;
  }
 
  return os;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
} // namespace Arcane::Alina
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
#endif