AMG.cc

// -*- 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.cc                                                      (C) 2000-2026 */
/*                                                                           */
/* Multi-grille algébrique.                                                  */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
#include "arcane/utils/Array.h"
#include "arcane/utils/ArgumentException.h"
#include "arcane/utils/FatalErrorException.h"
#include "arcane/utils/TraceAccessor.h"
#include "arcane/utils/OStringStream.h"
#include "arcane/utils/StringBuilder.h"
 
#include "arcane/core/matvec/Matrix.h"
 
#include <set>
#include <fstream>
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
namespace Arcane::math
{
Real divide(Real a, Real b)
{
  if (b == 0.0)
    ARCANE_FATAL("Division by zero");
  return a / b;
}
} // namespace Arcane::math
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
namespace Arcane::MatVec
{
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void DirectSolver::
solve(const Matrix& matrix, const Vector& vector_b, Vector& vector_x)
{
  IntegerConstArrayView rows = matrix.rowsIndex();
  IntegerConstArrayView columns = matrix.columns();
  RealConstArrayView mat_values = matrix.values();
 
  Integer nb_row = matrix.nbRow();
  RealUniqueArray solution_values(nb_row);
  RealUniqueArray full_matrix_values(nb_row * nb_row);
  full_matrix_values.fill(0.0);
  solution_values.copy(vector_b.values());
  for (Integer row = 0; row < nb_row; ++row) {
    for (Integer j = rows[row]; j < rows[row + 1]; ++j) {
      full_matrix_values[row * nb_row + columns[j]] = mat_values[j];
    }
  }
  _solve(full_matrix_values, solution_values, nb_row);
  vector_x.values().copy(solution_values);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void DirectSolver::
_solve(RealArrayView mat_values, RealArrayView vec_values, Integer size)
{
  if (size == 1) {
    if (math::isZero(mat_values[0]))
      throw FatalErrorException("DirectSolver", "Null matrix");
    vec_values[0] /= mat_values[0];
    return;
  }
 
  for (Integer k = 0; k < size - 1; ++k) {
    if (!math::isZero(mat_values[k * size + k])) {
      for (Integer j = k + 1; j < size; ++j) {
        if (!math::isZero(mat_values[j * size + k])) {
          Real factor = mat_values[j * size + k] / mat_values[k * size + k];
          for (Integer m = k + 1; m < size; ++m)
            mat_values[j * size + m] -= factor * mat_values[k * size + m];
          vec_values[j] -= factor * vec_values[k];
        }
      }
    }
  }
 
  for (Integer k = (size - 1); k > 0; --k) {
    vec_values[k] /= mat_values[k * size + k];
    for (Integer j = 0; j < k; ++j) {
      if (!math::isZero(mat_values[j * size + k]))
        vec_values[j] -= vec_values[k] * mat_values[j * size + k];
    }
  }
 
  vec_values[0] /= mat_values[0];
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix MatrixOperation2::
matrixMatrixProduct(const Matrix& left_matrix, const Matrix& right_matrix)
{
  Integer nb_left_col = left_matrix.nbColumn();
  Integer nb_right_col = right_matrix.nbColumn();
  Integer nb_right_row = right_matrix.nbRow();
  Integer nb_left_row = left_matrix.nbRow();
  if (nb_left_col != nb_right_row)
    ARCANE_THROW(ArgumentException, "Bad size nb_left_column={0} nb_right_row={1}",
                 nb_left_col, nb_right_row);
  Integer nb_row_col = nb_left_col;
 
  Matrix new_matrix(nb_left_row, nb_right_col);
  IntegerUniqueArray new_matrix_rows_size(nb_left_row);
  RealUniqueArray new_matrix_values;
  IntegerUniqueArray new_matrix_columns;
 
  for (Integer i = 0; i < nb_left_row; ++i) {
    Integer local_nb_col = 0;
    for (Integer j = 0; j < nb_right_col; ++j) {
      Real v = 0.0;
      for (Integer k = 0; k < nb_row_col; ++k) {
        //if (i==1 && j==0){
        // Real v0 = left_matrix.value(i,k) * right_matrix.value(k,j);
        //  cout << "** CHECK CONTRIBUTION k=" << k
        //       << " l=" << left_matrix.value(i,k) << " r=" << right_matrix.value(k,j) << '\n';
        //  if (!math::isZero(v0))
        //    cout << "** ADD CONTRIBUTION k=" << k << " v0=" << v0
        //         << " l=" << left_matrix.value(i,k) << " r=" << right_matrix.value(k,j) << '\n';
        // }
        v += left_matrix.value(i, k) * right_matrix.value(k, j);
      }
      if (!math::isZero(v)) {
        ++local_nb_col;
        new_matrix_columns.add(j);
        new_matrix_values.add(v);
      }
    }
    new_matrix_rows_size[i] = local_nb_col;
  }
  new_matrix.setRowsSize(new_matrix_rows_size);
  new_matrix.setValues(new_matrix_columns, new_matrix_values);
  return new_matrix;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix MatrixOperation2::
matrixMatrixProductFast(const Matrix& left_matrix, const Matrix& right_matrix)
{
  Integer nb_left_col = left_matrix.nbColumn();
  Integer nb_right_col = right_matrix.nbColumn();
  Integer nb_right_row = right_matrix.nbRow();
  Integer nb_left_row = left_matrix.nbRow();
  if (nb_left_col != nb_right_row)
    ARCANE_THROW(ArgumentException, "Bad size nb_left_column={0} nb_right_row={1}",
                 nb_left_col, nb_right_row);
  //Integer nb_row_col = nb_left_col;
 
  IntegerConstArrayView left_rows_index = left_matrix.rowsIndex();
  IntegerConstArrayView left_columns = left_matrix.columns();
  RealConstArrayView left_values = left_matrix.values();
 
  IntegerConstArrayView right_rows_index = right_matrix.rowsIndex();
  IntegerConstArrayView right_columns = right_matrix.columns();
  RealConstArrayView right_values = right_matrix.values();
 
  Matrix new_matrix(nb_left_row, nb_right_col);
  IntegerUniqueArray new_matrix_rows_size(nb_left_row);
  RealUniqueArray new_matrix_values;
  IntegerUniqueArray new_matrix_columns;
 
  IntegerUniqueArray col_right_columns_index(nb_right_col + 1);
  IntegerUniqueArray col_right_rows;
  RealUniqueArray col_right_values;
  IntegerUniqueArray col_right_columns_size(nb_right_col);
  {
    // Calcule le nombre d'éléments de chaque colonne
    col_right_columns_size.fill(0);
    for (Integer i = 0; i < nb_right_row; ++i) {
      for (Integer j = right_rows_index[i]; j < right_rows_index[i + 1]; ++j) {
        ++col_right_columns_size[right_columns[j]];
      }
    }
    // Calcule l'index du premier élément de chaque colonne.
    Integer index = 0;
    for (Integer j = 0; j < nb_right_col; ++j) {
      col_right_columns_index[j] = index;
      index += col_right_columns_size[j];
    }
    col_right_columns_index[nb_right_col] = index;
 
    col_right_rows.resize(index);
    col_right_values.resize(index);
    index = 0;
    // Remplit les valeurs par colonne
    col_right_columns_size.fill(0);
    for (Integer i = 0; i < nb_right_row; ++i) {
      for (Integer j = right_rows_index[i]; j < right_rows_index[i + 1]; ++j) {
        Integer col = right_columns[j];
        Real value = right_values[j];
        Integer col_index = col_right_columns_size[col] + col_right_columns_index[col];
        ++col_right_columns_size[col];
        col_right_rows[col_index] = i;
        col_right_values[col_index] = value;
      }
    }
  }
  //_dumpColumnMatrix(cout,col_right_columns_index,col_right_rows,col_right_values);
  //cout << '\n';
  RealUniqueArray current_row_values(nb_left_col);
  current_row_values.fill(0.0);
  for (Integer i = 0; i < nb_left_row; ++i) {
    Integer local_nb_col = 0;
    // Remplit la ligne avec les valeurs courantes
    for (Integer z = left_rows_index[i], zs = left_rows_index[i + 1]; z < zs; ++z) {
      current_row_values[left_columns[z]] = left_values[z];
      //if (i==1)
      //cout << " ** FILL VALUE col=" << left_columns[z] << " v=" << left_values[z] << '\n';
    }
    //if (i==1){
    //for( Integer z=0; z<nb_left_col; ++z ){
    //current_row_values[ left_columns[z] ] = left_values[z];
    //  cout << " ** VALUE col=" << z << " v=" << current_row_values[z] << '\n';
    //}
    //}
 
    for (Integer j = 0; j < nb_right_col; ++j) {
      Real v = 0.0;
      for (Integer zj = col_right_columns_index[j]; zj < col_right_columns_index[j + 1]; ++zj) {
        //if (i==1 && j==0){
        //  Real v0 = col_right_values[zj] * current_row_values[ col_right_rows[zj] ];
        //  cout << "** CHECK CONTRIBUTION2 k=" << col_right_rows[zj]
        //       << " l=" << current_row_values[ col_right_rows[zj] ] << " r=" << col_right_values[zj] << '\n';
        //  if (!math::isZero(v0))
        //    cout << "** ADD CONTRIBUTION2 k=" << col_right_rows[zj] << " v0=" << v0
        //         << " l=" << current_row_values[ col_right_rows[zj] ] << " r=" << col_right_values[zj] << '\n';
        // }
        v += col_right_values[zj] * current_row_values[col_right_rows[zj]];
      }
      if (!math::isZero(v)) {
        ++local_nb_col;
        new_matrix_columns.add(j);
        new_matrix_values.add(v);
      }
    }
 
    new_matrix_rows_size[i] = local_nb_col;
 
    // Remet des zeros dans la ligne courante.
    for (Integer z = left_rows_index[i], zs = left_rows_index[i + 1]; z < zs; ++z)
      current_row_values[left_columns[z]] = 0.0;
  }
  new_matrix.setRowsSize(new_matrix_rows_size);
  new_matrix.setValues(new_matrix_columns, new_matrix_values);
  return new_matrix;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixOperation2::
_dumpColumnMatrix(std::ostream& o, IntegerConstArrayView columns_index, IntegerConstArrayView rows,
                  RealConstArrayView values)
{
  Integer nb_col = columns_index.size() - 1;
  o << "(ColumnMatrix nb_col=" << nb_col;
  for (Integer j = 0; j < nb_col; ++j) {
    for (Integer z = columns_index[j], zs = columns_index[j + 1]; z < zs; ++z) {
      Integer i = rows[z];
      Real v = values[z];
      o << " [" << i << "," << j << "]=" << v;
    }
  }
  o << ")";
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix MatrixOperation2::
transpose(const Matrix& matrix)
{
  Integer nb_column = matrix.nbColumn();
  Integer nb_row = matrix.nbRow();
 
  Integer new_matrix_nb_row = nb_column;
  Integer new_matrix_nb_column = nb_row;
  Matrix new_matrix(new_matrix_nb_row, new_matrix_nb_column);
  IntegerUniqueArray new_matrix_rows_size(new_matrix_nb_row);
  RealUniqueArray new_matrix_values;
  IntegerUniqueArray new_matrix_columns;
 
  for (Integer i = 0; i < new_matrix_nb_row; ++i) {
    Integer local_nb_col = 0;
    for (Integer j = 0; j < new_matrix_nb_column; ++j) {
      Real v = matrix.value(j, i);
      if (!math::isZero(v)) {
        ++local_nb_col;
        new_matrix_columns.add(j);
        new_matrix_values.add(v);
      }
    }
    new_matrix_rows_size[i] = local_nb_col;
  }
  new_matrix.setRowsSize(new_matrix_rows_size);
  new_matrix.setValues(new_matrix_columns, new_matrix_values);
  return new_matrix;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix MatrixOperation2::
transposeFast(const Matrix& matrix)
{
  Integer nb_column = matrix.nbColumn();
  Integer nb_row = matrix.nbRow();
 
  IntegerConstArrayView rows_index = matrix.rowsIndex();
  IntegerConstArrayView columns = matrix.columns();
  RealConstArrayView values = matrix.values();
 
  Integer new_matrix_nb_row = nb_column;
  Integer new_matrix_nb_column = nb_row;
  Matrix new_matrix(new_matrix_nb_row, new_matrix_nb_column);
 
  IntegerUniqueArray new_matrix_rows_size(new_matrix_nb_row);
 
  // Calcul le nombre de colonnes de chaque ligne de la transposee.
  new_matrix_rows_size.fill(0);
  Integer nb_element = values.size();
  for (Integer i = 0, is = columns.size(); i < is; ++i) {
    ++new_matrix_rows_size[columns[i]];
  }
  new_matrix.setRowsSize(new_matrix_rows_size);
 
  IntegerConstArrayView new_matrix_rows_index = new_matrix.rowsIndex();
  new_matrix_rows_size.fill(0);
 
  RealUniqueArray new_matrix_values(nb_element);
  IntegerUniqueArray new_matrix_columns(nb_element);
 
  for (Integer row = 0, is = nb_row; row < is; ++row) {
    for (Integer j = rows_index[row]; j < rows_index[row + 1]; ++j) {
      Integer col_index = columns[j];
      Integer pos = new_matrix_rows_index[col_index] + new_matrix_rows_size[col_index];
      //cout << "** CURRENT row=" << row << " col=" << col_index << " v=" << values[col_index]
      new_matrix_columns[pos] = row;
      new_matrix_values[pos] = values[j];
      ++new_matrix_rows_size[col_index];
    }
  }
 
  new_matrix.setValues(new_matrix_columns, new_matrix_values);
  return new_matrix;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix MatrixOperation2::
applyGalerkinOperator(const Matrix& left_matrix, const Matrix& matrix,
                      const Matrix& right_matrix)
{
  Integer nb_original_row = matrix.nbRow();
  Integer nb_final_row = left_matrix.nbRow();
  IntegerUniqueArray p_marker(nb_final_row);
  IntegerUniqueArray a_marker(nb_original_row);
  p_marker.fill(-1);
  a_marker.fill(-1);
 
  IntegerConstArrayView left_matrix_rows = left_matrix.rowsIndex();
  IntegerConstArrayView left_matrix_columns = left_matrix.columns();
  RealConstArrayView left_matrix_values = left_matrix.values();
 
  IntegerConstArrayView right_matrix_rows = right_matrix.rowsIndex();
  IntegerConstArrayView right_matrix_columns = right_matrix.columns();
  RealConstArrayView right_matrix_values = right_matrix.values();
 
  IntegerConstArrayView matrix_rows = matrix.rowsIndex();
  IntegerConstArrayView matrix_columns = matrix.columns();
  RealConstArrayView matrix_values = matrix.values();
 
  Integer jj_counter = 0;
  Integer jj_row_begining = 0;
 
  IntegerUniqueArray new_matrix_rows_size(nb_final_row);
 
  // D'abord, détermine le nombre de colonnes de chaque ligne de la
  // matrice finale
  for (Integer ic = 0; ic < nb_final_row; ++ic) {
    // Ajoute la diagonale
    p_marker[ic] = jj_counter;
    jj_row_begining = jj_counter;
    ++jj_counter;
 
    // Boucle sur les colonnes de la ligne \a ic de \a matrix
    for (Integer jj1 = left_matrix_rows[ic]; jj1 < left_matrix_rows[ic + 1]; ++jj1) {
      Integer i1 = left_matrix_columns[jj1];
 
      // Boucle sur les colonnes de la ligne \a i1 de \a matrix
      for (Integer jj2 = matrix_rows[i1]; jj2 < matrix_rows[i1 + 1]; ++jj2) {
        Integer i2 = matrix_columns[jj2];
        /*--------------------------------------------------------------
         *  Vérifie A_marker pour voir si le point i2 a déjà été
         *  visité. Les nouvelles entrées dans RAP ne proviennent que
         *  des points non marqués.
         *--------------------------------------------------------------*/
        if (a_marker[i2] != ic) {
          a_marker[i2] = ic;
          /*-----------------------------------------------------------
           *  Boucle sur les entrées de la ligne i2 de P.
           *-----------------------------------------------------------*/
          for (Integer jj3 = right_matrix_rows[i2]; jj3 < right_matrix_rows[i2 + 1]; ++jj3) {
            Integer i3 = right_matrix_columns[jj3];
            /*--------------------------------------------------------
             *  Vérifie P_marker pour s'assurer que RAP_{ic,i3} n'a pas déjà
             *  été pris en compte. Si ce n'est pas le cas, le marque et incrémente
             *  le compteur.
             *--------------------------------------------------------*/
            if (p_marker[i3] < jj_row_begining) {
              p_marker[i3] = jj_counter;
              ++jj_counter;
            }
          }
        }
      }
    }
    new_matrix_rows_size[ic] = jj_counter - jj_row_begining;
  }
  static Integer total_rap_size = 0;
  total_rap_size += jj_counter;
 
  std::cout << "** RAP_SIZE=" << jj_counter << " TOTAL=" << total_rap_size << '\n';
  Matrix new_matrix(nb_final_row, nb_final_row);
  new_matrix.setRowsSize(new_matrix_rows_size);
 
  //IntegerConstArrayView new_matrix_rows = new_matrix.rowsIndex();
  IntegerArrayView new_matrix_columns = new_matrix.columns();
  RealArrayView new_matrix_values = new_matrix.values();
 
  // Maintenant, remplit les coefficients de la matrice
  p_marker.fill(-1);
  a_marker.fill(-1);
  jj_counter = 0;
  for (Integer ic = 0; ic < nb_final_row; ++ic) {
    // Ajoute la diagonale
    p_marker[ic] = jj_counter;
    jj_row_begining = jj_counter;
    new_matrix_columns[jj_counter] = ic;
    new_matrix_values[jj_counter] = 0.0;
    ++jj_counter;
    // Boucle sur les colonnes de la ligne \a ic de \a matrix
    for (Integer jj1 = left_matrix_rows[ic]; jj1 < left_matrix_rows[ic + 1]; ++jj1) {
      Integer i1 = left_matrix_columns[jj1];
      Real r_entry = left_matrix_values[jj1];
 
      // Boucle sur les colonnes de la ligne \a i1 de \a matrix
      for (Integer jj2 = matrix_rows[i1]; jj2 < matrix_rows[i1 + 1]; ++jj2) {
        Integer i2 = matrix_columns[jj2];
        Real r_a_product = r_entry * matrix_values[jj2];
        /*--------------------------------------------------------------
         *  Vérifie A_marker pour voir si le point i2 a déjà été
         *  visité. Les nouvelles entrées dans RAP ne proviennent que
         *  des points non marqués.
         *--------------------------------------------------------------*/
        if (a_marker[i2] != ic) {
          a_marker[i2] = ic;
          /*-----------------------------------------------------------
           *  Boucle sur les entrées de la ligne i2 de P.
           *-----------------------------------------------------------*/
          for (Integer jj3 = right_matrix_rows[i2]; jj3 < right_matrix_rows[i2 + 1]; ++jj3) {
            Integer i3 = right_matrix_columns[jj3];
            Real r_a_p_product = r_a_product * right_matrix_values[jj3];
            /*--------------------------------------------------------
             *  Vérifie P_marker pour s'assurer que RAP_{ic,i3} n'a pas déjà
             *  été pris en compte. Si ce n'est pas le cas, crée une nouvelle entrée.
             *  S'il l'est, ajoute la nouvelle contribution.
             *--------------------------------------------------------*/
            if (p_marker[i3] < jj_row_begining) {
              p_marker[i3] = jj_counter;
              new_matrix_values[jj_counter] = r_a_p_product;
              new_matrix_columns[jj_counter] = i3;
              ++jj_counter;
            }
            else {
              new_matrix_values[p_marker[i3]] += r_a_p_product;
            }
          }
        }
        else {
          /*--------------------------------------------------------------
           *  Si i2 a déjà été visité ( A_marker[12]=ic ) cela ne produit
           *  aucune nouvelle entrée dans RAP et peut simplement ajouter de nouvelles contributions.
           *--------------------------------------------------------------*/
          for (Integer jj3 = right_matrix_rows[i2]; jj3 < right_matrix_rows[i2 + 1]; ++jj3) {
            Integer i3 = right_matrix_columns[jj3];
            Real r_a_p_product = r_a_product * right_matrix_values[jj3];
            new_matrix_values[p_marker[i3]] += r_a_p_product;
          }
        }
      }
    }
  }
  return new_matrix;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix MatrixOperation2::
applyGalerkinOperator2(const Matrix& left_matrix, const Matrix& matrix,
                       const Matrix& right_matrix)
{
  Integer nb_original_row = matrix.nbRow();
  Integer nb_final_row = left_matrix.nbRow();
  IntegerUniqueArray p_marker(nb_final_row);
  IntegerUniqueArray a_marker(nb_original_row);
  p_marker.fill(-1);
  a_marker.fill(-1);
 
  const Integer* left_matrix_rows = left_matrix.rowsIndex().data();
  const Integer* left_matrix_columns = left_matrix.columns().data();
  const Real* left_matrix_values = left_matrix.values().data();
 
  const Integer* right_matrix_rows = right_matrix.rowsIndex().data();
  const Integer* right_matrix_columns = right_matrix.columns().data();
  const Real* right_matrix_values = right_matrix.values().data();
 
  const Integer* matrix_rows = matrix.rowsIndex().data();
  const Integer* matrix_columns = matrix.columns().data();
  const Real* matrix_values = matrix.values().data();
 
  Integer jj_counter = 0;
  Integer jj_row_begining = 0;
 
  IntegerUniqueArray new_matrix_rows_size(nb_final_row);
 
  // D'abord, détermine le nombre de colonnes de chaque ligne de la
  // matrice finale
  for (Integer ic = 0; ic < nb_final_row; ++ic) {
    // Ajoute la diagonale
    p_marker[ic] = jj_counter;
    jj_row_begining = jj_counter;
    ++jj_counter;
 
    // Boucle sur les colonnes de la ligne \a ic de \a matrix
    for (Integer jj1 = left_matrix_rows[ic]; jj1 < left_matrix_rows[ic + 1]; ++jj1) {
      Integer i1 = left_matrix_columns[jj1];
 
      // Boucle sur les colonnes de la ligne \a i1 de \a matrix
      for (Integer jj2 = matrix_rows[i1]; jj2 < matrix_rows[i1 + 1]; ++jj2) {
        Integer i2 = matrix_columns[jj2];
        /*--------------------------------------------------------------
         *  Vérifie A_marker pour voir si le point i2 a déjà été
         *  visité. Les nouvelles entrées dans RAP ne proviennent que
         *  des points non marqués.
         *--------------------------------------------------------------*/
        if (a_marker[i2] != ic) {
          a_marker[i2] = ic;
          /*-----------------------------------------------------------
           *  Boucle sur les entrées de la ligne i2 de P.
           *-----------------------------------------------------------*/
          for (Integer jj3 = right_matrix_rows[i2]; jj3 < right_matrix_rows[i2 + 1]; ++jj3) {
            Integer i3 = right_matrix_columns[jj3];
            /*--------------------------------------------------------
             *  Vérifie P_marker pour s'assurer que RAP_{ic,i3} n'a pas déjà
             *  été pris en compte. Si ce n'est pas le cas, le marque et incrémente
             *  le compteur.
             *--------------------------------------------------------*/
            if (p_marker[i3] < jj_row_begining) {
              p_marker[i3] = jj_counter;
              ++jj_counter;
            }
          }
        }
      }
    }
    new_matrix_rows_size[ic] = jj_counter - jj_row_begining;
  }
 
  Matrix new_matrix(nb_final_row, nb_final_row);
  new_matrix.setRowsSize(new_matrix_rows_size);
 
  //Integer* new_matrix_rows = new_matrix.rowsIndex();
  Integer* ARCANE_RESTRICT new_matrix_columns = new_matrix.columns().data();
  Real* ARCANE_RESTRICT new_matrix_values = new_matrix.values().data();
 
  // Maintenant, remplit les coefficients de la matrice
  p_marker.fill(-1);
  a_marker.fill(-1);
  jj_counter = 0;
  for (Integer ic = 0; ic < nb_final_row; ++ic) {
    // Ajoute la diagonale
    p_marker[ic] = jj_counter;
    jj_row_begining = jj_counter;
    new_matrix_columns[jj_counter] = ic;
    new_matrix_values[jj_counter] = 0.0;
    ++jj_counter;
    // Boucle sur les colonnes de la ligne \a ic de \a matrix
    for (Integer jj1 = left_matrix_rows[ic]; jj1 < left_matrix_rows[ic + 1]; ++jj1) {
      Integer i1 = left_matrix_columns[jj1];
      Real r_entry = left_matrix_values[jj1];
 
      // Boucle sur les colonnes de la ligne \a i1 de \a matrix
      for (Integer jj2 = matrix_rows[i1]; jj2 < matrix_rows[i1 + 1]; ++jj2) {
        Integer i2 = matrix_columns[jj2];
        Real r_a_product = r_entry * matrix_values[jj2];
        /*--------------------------------------------------------------
         *  Vérifie A_marker pour voir si le point i2 a déjà été
         *  visité. Les nouvelles entrées dans RAP ne proviennent que
         *  des points non marqués.
         *--------------------------------------------------------------*/
        if (a_marker[i2] != ic) {
          a_marker[i2] = ic;
          /*-----------------------------------------------------------
           *  Boucle sur les entrées de la ligne i2 de P.
           *-----------------------------------------------------------*/
          for (Integer jj3 = right_matrix_rows[i2]; jj3 < right_matrix_rows[i2 + 1]; ++jj3) {
            Integer i3 = right_matrix_columns[jj3];
            Real r_a_p_product = r_a_product * right_matrix_values[jj3];
            /*--------------------------------------------------------
             *  Vérifie P_marker pour s'assurer que RAP_{ic,i3} n'a pas déjà
             *  été pris en compte. Si ce n'est pas le cas, crée une nouvelle entrée.
             *  S'il l'est, ajoute la nouvelle contribution.
             *--------------------------------------------------------*/
            if (p_marker[i3] < jj_row_begining) {
              p_marker[i3] = jj_counter;
              new_matrix_values[jj_counter] = r_a_p_product;
              new_matrix_columns[jj_counter] = i3;
              ++jj_counter;
            }
            else {
              new_matrix_values[p_marker[i3]] += r_a_p_product;
            }
          }
        }
        else {
          /*--------------------------------------------------------------
           *  Si i2 a déjà été visité ( A_marker[12]=ic ) cela ne produit
           *  aucune nouvelle entrée dans RAP et peut simplement ajouter de nouvelles contributions.
           *--------------------------------------------------------------*/
          for (Integer jj3 = right_matrix_rows[i2]; jj3 < right_matrix_rows[i2 + 1]; ++jj3) {
            Integer i3 = right_matrix_columns[jj3];
            Real r_a_p_product = r_a_product * right_matrix_values[jj3];
            new_matrix_values[p_marker[i3]] += r_a_p_product;
          }
        }
      }
    }
  }
  return new_matrix;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
enum
{
  TYPE_UNDEFINED = 0,
  TYPE_COARSE = 1,
  TYPE_FINE = 2,
  TYPE_SPECIAL_FINE = 3
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
class AMGLevel
: public TraceAccessor
{
 public:
 
  AMGLevel(ITraceMng* tm, Integer level)
  : TraceAccessor(tm)
  , m_level(level)
  , m_is_verbose(false)
  {}
  virtual ~AMGLevel() {}
 
 public:
 
  virtual void buildLevel(Matrix matrix, Real alpha);
 
 public:
 
  Matrix fineMatrix()
  {
    return m_fine_matrix;
  }
  Matrix coarseMatrix()
  {
    return m_coarse_matrix;
  }
  Matrix prolongationMatrix()
  {
    return m_prolongation_matrix;
  }
  Matrix restrictionMatrix()
  {
    return m_restriction_matrix;
  }
  Integer nbCoarsePoint() const
  {
    return m_coarse_matrix.nbRow();
  }
  Int32ConstArrayView pointsType() const
  {
    return m_points_type;
  }
  void printLevelInfo();
 
 private:
 
  Integer m_level;
  Matrix m_fine_matrix;
  Matrix m_coarse_matrix;
  Matrix m_prolongation_matrix;
  Matrix m_restriction_matrix;
  Int32UniqueArray m_points_type;
 
  bool m_is_verbose;
 
 private:
 
  void _buildCoarsePoints(Real alpha,
                          RealArray& rows_max_val,
                          UniqueArray<SharedArray<Integer>>& depends,
                          IntegerArray& weak_depends);
  void _buildInterpolationMatrix(RealConstArrayView rows_max_val,
                                 UniqueArray<SharedArray<Integer>>& depends,
                                 IntegerArray& weak_depends);
  void _printLevelInfo(Matrix matrix);
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
class AMG
: public TraceAccessor
{
 public:
 
  AMG(ITraceMng* tm)
  : TraceAccessor(tm)
  {}
  ~AMG();
 
 public:
 
  void build(Matrix matrix);
  void solve(const Vector& vector_b, Vector& vector_x);
 
 private:
 
  UniqueArray<AMGLevel*> m_levels;
  Matrix m_matrix;
 
 private:
 
  void _solve(const Vector& vector_b, Vector& vector_x, Integer level);
  void _relax(const Matrix& matrix, const Vector& vector_b, Vector& vector_x,
              Integer nb_relax);
  void _relax1(const Matrix& matrix, const Vector& vector_b, Vector& vector_x,
               Integer nb_relax);
  void _relaxJacobi(const Matrix& matrix, const Vector& vector_b, Vector& vector_x,
                    Real weight);
  void _relaxGaussSeidel(const Matrix& matrix, const Vector& vector_b, Vector& vector_x,
                         Integer point_type, Int32ConstArrayView points_type);
  void _relaxSymmetricGaussSeidel(const Matrix& matrix, const Vector& vector_b, Vector& vector_x);
  void _printResidualInfo(const Matrix& matrix, const Vector& vector_b,
                          const Vector& vector_x);
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
AMG::
~AMG()
{
  for (Integer i = 0; i < m_levels.size(); ++i)
    delete m_levels[i];
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
build(Matrix matrix)
{
  Matrix current_matrix = matrix;
  m_matrix = matrix;
  for (Integer i = 1; i < 100; ++i) {
    AMGLevel* level = new AMGLevel(traceMng(), i);
    level->buildLevel(current_matrix, 0.25);
    m_levels.add(level);
    Integer nb_coarse_point = level->nbCoarsePoint();
    if (nb_coarse_point < 20)
      break;
    current_matrix = level->coarseMatrix();
  }
  //for( Integer i=0; i<m_levels.size(); ++i )
  //m_levels[i]->printLevelInfo();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
solve(const Vector& vector_b, Vector& vector_x)
{
  //info() << "AMG::solve";
  if (0) {
    OStringStream ostr;
    RealConstArrayView v_values(vector_b.values());
    for (Integer i = 0; i < 20; ++i)
      ostr() << "VECTOR_F_" << i << " = " << v_values[i] << " X=" << vector_x.values()[i] << '\n';
    for (Integer i = 0; i < v_values.size(); ++i)
      if (math::abs(v_values[i]) > 1e-5)
        ostr() << "VECTOR_F_" << i << " = " << v_values[i] << '\n';
    info() << "VECTOR_F\n"
           << ostr.str();
  }
  //_printResidualInfo(m_matrix,vector_b,vector_x);
  _solve(vector_b, vector_x, 0);
  //info() << "END SOLVE";
  //_printResidualInfo(m_matrix,vector_b,vector_x);
  //info() << "END AMG::solve";
  /*{
    OStringStream ostr;
    ostr() << "\nVECTOR_X ";
    vector_x.dump(ostr());
    info() << ostr.str();
    }*/
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
_relax1(const Matrix& matrix, const Vector& vector_b, Vector& vector_x,
        Integer nb_relax)
{
  Vector r(vector_x.size());
  MatrixOperation mat_op;
  for (Integer i = 0; i < nb_relax; ++i) {
    // r = b - A * x
    mat_op.matrixVectorProduct(matrix, vector_x, r);
    mat_op.negateVector(r);
    mat_op.addVector(r, vector_b);
 
    // x = x + r
    mat_op.addVector(vector_x, r);
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
_relax(const Matrix& matrix, const Vector& vector_b, Vector& vector_x,
       Integer nb_relax)
{
  Real epsilon = 1.0e-10;
  DiagonalPreconditioner p(matrix);
  ConjugateGradientSolver solver;
  solver.setMaxIteration(nb_relax);
  //mat_op.matrixVectorProduct(restriction_matrix,vector_x,new_x);
  //new_x.values().fill(0.0);
 
  //OStringStream ostr;
  //vector_x.dump(ostr());
  //info() << " RELAX BEFORE VECTOR_X=" << ostr.str();
 
  solver.solve(matrix, vector_b, vector_x, epsilon, &p);
  //OStringStream ostr;
  //ostr() << " COARSE_B=";
  //new_b.dump(ostr());
  //ostr() << "\nCOARSE_X=";
  //new_x.dump(ostr());
  //info() << "SOLVE COARSE MATRIX nb_iter=" << solver.nbIteration()
  //       << ostr.str();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
_relaxJacobi(const Matrix& matrix, const Vector& vector_b, Vector& vector_x,
             Real weight)
{
  IntegerConstArrayView rows = matrix.rowsIndex();
  IntegerConstArrayView columns = matrix.columns();
  RealConstArrayView mat_values = matrix.values();
 
  RealArrayView x_values = vector_x.values();
  RealConstArrayView b_values = vector_b.values();
 
  Integer nb_row = matrix.nbRow();
  RealConstArrayView cx_values(vector_x.values());
  RealUniqueArray tmp_values(cx_values);
  if (0) {
    Integer v = math::min(40, nb_row);
    OStringStream ostr;
    for (Integer i = (nb_row - 1); i > (nb_row - v); --i)
      ostr() << "BEFORE_B=" << i << "=" << b_values[i] << " U=" << x_values[i] << " T=" << tmp_values[i] << '\n';
    info() << "B = X=" << x_values.data() << " T=" << tmp_values.data() << "\n"
           << ostr.str();
  }
  Real one_minus_weight = 1.0 - weight;
  for (Integer row = 0; row < nb_row; ++row) {
    Real diag = mat_values[rows[row]];
    if (math::isZero(diag))
      continue;
    Real res = b_values[row];
    for (Integer j = rows[row] + 1; j < rows[row + 1]; ++j) {
      Integer col = columns[j];
      res -= mat_values[j] * tmp_values[col];
    }
    x_values[row] *= one_minus_weight;
    x_values[row] += (weight * res) / diag;
  }
  if (0) {
    Integer v = math::min(40, nb_row);
    OStringStream ostr;
    for (Integer i = (nb_row - 1); i > (nb_row - v); --i)
      ostr() << "AFTER_B=" << i << "=" << b_values[i] << " U=" << x_values[i] << '\n';
    info() << "B\n"
           << ostr.str();
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
_relaxGaussSeidel(const Matrix& matrix, const Vector& vector_b, Vector& vector_x,
                  Integer point_type, Int32ConstArrayView points_type2)
{
#if 1
  IntegerConstArrayView rows = matrix.rowsIndex();
  IntegerConstArrayView columns = matrix.columns();
  RealConstArrayView mat_values = matrix.values();
 
  RealArrayView x_values = vector_x.values();
  RealConstArrayView b_values = vector_b.values();
  Int32ConstArrayView points_type = points_type2;
#else
  const Integer* rows = matrix.rowsIndex().data();
  const Integer* columns = matrix.columns().data();
  const Real* mat_values = matrix.values().data();
 
  Real* ARCANE_RESTRICT x_values = vector_x.values().data();
  const Real* b_values = vector_b.values().data();
  const Integer* points_type = points_type2.data();
#endif
 
  Integer nb_row = matrix.nbRow();
  if (0) {
    info() << " RELAX nb_relax=" << " nb_row=" << nb_row
           << " point_type=" << point_type;
    Integer v = math::min(40, nb_row);
    OStringStream ostr;
    for (Integer i = (nb_row - 1); i > (nb_row - v); --i)
      ostr() << "BEFORE_B=" << i << "=" << b_values[i] << " U=" << x_values[i] << '\n';
    info() << "B\n"
           << ostr.str();
  }
  for (Integer row = 0; row < nb_row; ++row) {
    Real diag = mat_values[rows[row]];
    if (points_type[row] != point_type || math::isZero(diag))
      continue;
    Real res = b_values[row];
    for (Integer j = rows[row] + 1; j < rows[row + 1]; ++j) {
      Integer col = columns[j];
      res -= mat_values[j] * x_values[col];
    }
    x_values[row] = res / diag;
  }
  if (0) {
    Integer v = math::min(40, nb_row);
    OStringStream ostr;
    for (Integer i = (nb_row - 1); i > (nb_row - v); --i)
      ostr() << "AFTER_B=" << i << "=" << b_values[i] << " U=" << x_values[i] << '\n';
    info() << "B\n"
           << ostr.str();
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
_relaxSymmetricGaussSeidel(const Matrix& matrix, const Vector& vector_b, Vector& vector_x)
{
  IntegerConstArrayView rows = matrix.rowsIndex();
  IntegerConstArrayView columns = matrix.columns();
  RealConstArrayView mat_values = matrix.values();
 
  RealArrayView x_values = vector_x.values();
  RealConstArrayView b_values = vector_b.values();
 
  Integer nb_row = matrix.nbRow();
  //info() << " RELAX nb_relax=" << nb_relax << " nb_row=" << nb_row
  //       << " point_type=" << point_type;
  for (Integer row = 0; row < nb_row; ++row) {
    Real diag = mat_values[rows[row]];
    if (math::isZero(diag))
      continue;
    Real res = b_values[row];
    for (Integer j = rows[row] + 1; j < rows[row + 1]; ++j) {
      Integer col = columns[j];
      res -= mat_values[j] * x_values[col];
    }
    x_values[row] = res / diag;
  }
 
  for (Integer row = nb_row - 1; row > -1; --row) {
    Real diag = mat_values[rows[row]];
    if (math::isZero(diag))
      continue;
    Real res = b_values[row];
    for (Integer j = rows[row] + 1; j < rows[row + 1]; ++j) {
      Integer col = columns[j];
      res -= mat_values[j] * x_values[col];
    }
    x_values[row] = res / diag;
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
_solve(const Vector& vector_b, Vector& vector_x, Integer level)
{
  AMGLevel* current_level = m_levels[level];
  Integer vector_size = vector_b.size();
  Integer nb_coarse = current_level->nbCoarsePoint();
  Matrix fine_matrix = current_level->fineMatrix();
  Matrix restriction_matrix = current_level->restrictionMatrix();
  Matrix coarse_matrix = current_level->coarseMatrix();
  Matrix prolongation_matrix = current_level->prolongationMatrix();
 
  Integer new_nb_row = nb_coarse;
  Vector new_b(new_nb_row);
  Vector new_x(new_nb_row);
  Vector tmp(vector_size);
 
  MatrixOperation mat_op;
 
  bool is_final_level = (level + 1) == m_levels.size();
 
  const bool use_gauss_seidel = false;
  Integer nb_relax1 = 2;
  Real jacobi_weight = 2.0 / 3.0;
  if (use_gauss_seidel) {
    for (Integer i = 0; i < nb_relax1; ++i)
      _relaxGaussSeidel(fine_matrix, vector_b, vector_x, TYPE_FINE, current_level->pointsType());
    for (Integer i = 0; i < nb_relax1; ++i)
      _relaxGaussSeidel(fine_matrix, vector_b, vector_x, TYPE_COARSE, current_level->pointsType());
  }
  else {
    //info() << "BEFORE SMOOTH";
    //_printResidualInfo(fine_matrix,vector_b,vector_x);
    for (Integer i = 0; i < nb_relax1; ++i) {
      //_relaxSymmetricGaussSeidel(fine_matrix,vector_b,vector_x);
      _relaxJacobi(fine_matrix, vector_b, vector_x, jacobi_weight);
    }
    //info() << "AFTER SMOOTH";
    //_printResidualInfo(fine_matrix,vector_b,vector_x);
    //_relax(fine_matrix,vector_b,vector_x,nb_relax1);
  }
 
  // Restreint le nouveau b à partir du b actuel
  // b(k+1) = I * (b(k) - A * x)
  {
    OStringStream ostr;
    mat_op.matrixVectorProduct(fine_matrix, vector_x, tmp);
    //ostr() << "\nCOARSE_B TMP(A*x) level=" << level << " ";
    //tmp.dump(ostr());
    mat_op.negateVector(tmp);
    mat_op.addVector(tmp, vector_b);
    //ostr() << "\nCOARSE_B TMP(b-A*x) level=" << level << " ";
    //tmp.dump(ostr());
    mat_op.matrixVectorProduct(restriction_matrix, tmp, new_b);
    //ostr() << "\nCOARSE_B level=" << level << " ";
    //new_b.dump(ostr());
    info() << ostr.str();
  }
 
  //mat_op.matrixVectorProduct(transpose_prolongation_matrix,vector_x,tmp);
 
  // Si le niveau final est atteint, résoudre la matrice.
  // Sinon, continuer en restreignant la matrice à nouveau
  if (is_final_level) {
 
    //info() << " SOLVE FINAL LEVEL";
 
    if (1) {
      DirectSolver ds;
      ds.solve(coarse_matrix, new_b, new_x);
      //_printResidualInfo(coarse_matrix,new_b,new_x);
    }
    else {
      Real epsilon = 1.0e-14;
      DiagonalPreconditioner p(coarse_matrix);
      ConjugateGradientSolver solver;
      //mat_op.matrixVectorProduct(restriction_matrix,vector_x,new_x);
      new_x.values().fill(0.0);
      solver.solve(coarse_matrix, new_b, new_x, epsilon, &p);
      OStringStream ostr;
      //ostr() << " COARSE_B=";
      //new_b.dump(ostr());
      //ostr() << "\nCOARSE_X=";
      //new_x.dump(ostr());
      //if (m_is_verbose)
      info() << "SOLVE COARSE MATRIX nb_iter=" << solver.nbIteration();
      //       << ostr.str();
      //_printResidualInfo(coarse_matrix,new_b,new_x);
    }
  }
  else {
    new_x.values().fill(0.0);
    _solve(new_b, new_x, level + 1);
  }
 
  // Interpole le nouveau x à partir de la solution trouvée
  // x(k) = x(k) + tI * x(k+1)
  mat_op.matrixVectorProduct(prolongation_matrix, new_x, tmp);
  mat_op.addVector(vector_x, tmp);
  /*{
    OStringStream ostr;
    vector_x.dump(ostr());
    info() << "NEW_X level=" << level << " X=" << ostr.str();
    _printResidualInfo(fine_matrix,vector_b,vector_x);
    }*/
 
  // Relaxation de Richardson
  if (use_gauss_seidel) {
    for (Integer i = 0; i < nb_relax1; ++i)
      _relaxGaussSeidel(fine_matrix, vector_b, vector_x, TYPE_FINE, current_level->pointsType());
    for (Integer i = 0; i < nb_relax1; ++i)
      _relaxGaussSeidel(fine_matrix, vector_b, vector_x, TYPE_COARSE, current_level->pointsType());
  }
  else {
    //info() << "BEFORE SMOOTH 2";
    //_printResidualInfo(fine_matrix,vector_b,vector_x);
    for (Integer i = 0; i < nb_relax1; ++i) {
      //_relaxSymmetricGaussSeidel(fine_matrix,vector_b,vector_x);
      _relaxJacobi(fine_matrix, vector_b, vector_x, jacobi_weight);
    }
    //info() << "AFTER SMOOTH 2";
    //_printResidualInfo(fine_matrix,vector_b,vector_x);
    //_relax(fine_matrix,vector_b,vector_x,nb_relax2);
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMG::
_printResidualInfo(const Matrix& a, const Vector& b, const Vector& x)
{
  OStringStream ostr;
  Vector tmp(b.size());
  // tmp = b - Ax
  MatrixOperation mat_op;
  mat_op.matrixVectorProduct(a, x, tmp);
  //ostr() << "\nAX=";
  //tmp.dump(ostr());
  mat_op.negateVector(tmp);
  mat_op.addVector(tmp, b);
  Real r = mat_op.dot(tmp);
  if (0) {
    Integer v = math::min(10, tmp.size());
    for (Integer i = 0; i < v; ++i)
      info() << "R_" << i << " = " << tmp.values()[i];
  }
  info() << " AMG_RESIDUAL_NORM=" << r << " sqrt=" << math::sqrt(r);
 
  //ostr() << "\nR=";
  //tmp.dump(ostr());
  //info() << " AMG_RESIDUAL_NORM="  << r <<  " AMG_RESIDUAL=" << ostr.str();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
class PointInfo
{
 public:
 
  PointInfo()
  : m_lambda(0)
  , m_index(0)
  {}
  PointInfo(Integer lambda, Integer index)
  : m_lambda(lambda)
  , m_index(index)
  {}
  Integer m_lambda;
  Integer m_index;
  bool operator<(const PointInfo& rhs) const
  {
    if (m_lambda == rhs.m_lambda)
      return m_index < rhs.m_index;
    return (m_lambda > rhs.m_lambda);
  }
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMGLevel::
printLevelInfo()
{
  _printLevelInfo(m_prolongation_matrix);
  _printLevelInfo(m_coarse_matrix);
}
 
void AMGLevel::
_printLevelInfo(Matrix matrix)
{
  OStringStream ostr;
  Integer nb_row = matrix.nbRow();
  Integer nb_column = matrix.nbColumn();
 
  IntegerConstArrayView rows = matrix.rowsIndex();
  //IntegerConstArrayView columns = matrix.columns();
  RealConstArrayView values = matrix.values();
  Integer nb_value = values.size();
 
  Real max_val = 0.0;
  Real min_val = 0.0;
  if (nb_value > 0) {
    max_val = values[0];
    min_val = values[0];
  }
 
  Real max_row_sum = 0.0;
  Real min_row_sum = 0.0;
  for (Integer row = 0; row < nb_row; ++row) {
    Real row_sum = 0.0;
    for (Integer z = rows[row], zs = rows[row + 1]; z < zs; ++z) {
      //Integer col = columns[z];
      Real v = values[z];
      if (v > max_val)
        max_val = v;
      if (v < min_val)
        min_val = v;
      row_sum += v;
    }
    if (row == 0) {
      max_row_sum = row_sum;
      min_row_sum = row_sum;
    }
    if (row_sum > max_row_sum)
      max_row_sum = row_sum;
    if (row_sum < max_row_sum)
      min_row_sum = row_sum;
  }
 
  Real sparsity = ((Real)nb_value) / ((Real)nb_row * (Real)nb_column);
 
  ostr() << "level=" << m_level
         << " nb_row=" << nb_row
         << " nb_col=" << nb_column
         << " nb_nonzero=" << nb_value
         << " sparsity=" << sparsity
         << " min=" << min_val
         << " max=" << max_val
         << " min_row=" << min_row_sum
         << " max_row=" << max_row_sum;
 
  info() << "INFO: " << ostr.str();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMGLevel::
_buildCoarsePoints(Real alpha,
                   RealArray& rows_max_val,
                   UniqueArray<SharedArray<Integer>>& depends,
                   IntegerArray& weak_depends)
{
  IntegerConstArrayView rows_index = m_fine_matrix.rowsIndex();
  IntegerConstArrayView columns = m_fine_matrix.columns();
  RealConstArrayView mat_values = m_fine_matrix.values();
  Integer nb_row = m_fine_matrix.nbRow();
 
  Int32UniqueArray lambdas(nb_row);
  lambdas.fill(0);
  UniqueArray<SharedArray<Integer>> influences(nb_row);
  depends.resize(nb_row);
  //UniqueArray<IntegerUniqueArray> weak_depends(nb_row);
  m_points_type.resize(nb_row);
  m_points_type.fill(TYPE_UNDEFINED);
 
  weak_depends.resize(mat_values.size());
  weak_depends.fill(0);
 
  const bool type_hypre = true;
  // Valeurs de chaque ligne qui influence
  rows_max_val.resize(nb_row);
  for (Integer row = 0; row < nb_row; ++row) {
    Real max_val = 0.0;
    Real min_val = 0.0;
    Real diag_val = mat_values[rows_index[row]];
    // Cherche le max (en valeur absolue) de la colonne, autre que la diagonale
    for (Integer z = rows_index[row] + 1, zs = rows_index[row + 1]; z < zs; ++z) {
      //Real mv = math::abs(mat_values[z]);
      Real mv = mat_values[z];
      if (!type_hypre)
        mv = math::abs(mv);
      if (mv > max_val)
        max_val = mv;
      if (mv < min_val)
        min_val = mv;
    }
    // Prend tous les éléments supérieurs à alpha * max_val
    //rows_max_val[row] = max_val * alpha;
    if (type_hypre) {
      if (diag_val < 0.0)
        rows_max_val[row] = max_val * alpha;
      else
        rows_max_val[row] = min_val * alpha;
    }
    else
      rows_max_val[row] = max_val * alpha;
  }
 
  for (Integer row = 0; row < nb_row; ++row) {
    // Prend tous les éléments supérieurs à max_val
    Real max_val = rows_max_val[row];
    Real diag_val = mat_values[rows_index[row]];
    for (Integer z = rows_index[row] + 1, zs = rows_index[row + 1]; z < zs; ++z) {
      //Real mv = math::abs(mat_values[z]);
      Real mv = mat_values[z];
      if (type_hypre) {
        if (diag_val < 0.0) {
          if (mv > max_val) {
            Integer column = columns[z];
            if (m_is_verbose)
              info() << " ADD INFLUENCE: ROW=" << row << " COL=" << column;
            ++lambdas[column];
            depends[row].add(column);
            influences[column].add(row);
            weak_depends[z] = 2;
          }
          else
            weak_depends[z] = 1;
        }
        else {
          if (mv < max_val) {
            Integer column = columns[z];
            if (m_is_verbose)
              info() << " ADD INFLUENCE: ROW=" << row << " COL=" << column;
            ++lambdas[column];
            depends[row].add(column);
            influences[column].add(row);
            weak_depends[z] = 2;
          }
          else
            weak_depends[z] = 1;
        }
      }
      else {
        if (math::abs(mv) > max_val) {
          Integer column = columns[z];
          if (m_is_verbose)
            info() << " ADD INFLUENCE: ROW=" << row << " COL=" << column;
          ++lambdas[column];
          depends[row].add(column);
          influences[column].add(row);
          weak_depends[z] = 2;
        }
        else {
          weak_depends[z] = 1;
        }
      }
      //else
      //weak_depends[row].add(column);
    }
  }
 
  if (0) {
    OStringStream ostr;
    Integer index = 0;
    int n = math::min(nb_row, 800);
    ostr() << "GRAPH\n";
    for (Integer i = 0; i < n; ++i) {
      ostr() << " GRAPH I=" << i << " ";
      for (Integer j = 0; j < depends[i].size(); ++j) {
        ++index;
        ostr() << " " << depends[i][j];
      }
      ostr() << " index=" << index << '\n';
    }
    ostr() << "\n MAXTRIX\n";
    index = 0;
    for (Integer i = 0; i < n; ++i) {
      ostr() << "MATRIX I=" << i << " ";
      for (Integer j = rows_index[i]; j < rows_index[i + 1]; ++j) {
        ++index;
        ostr() << " " << columns[j] << " " << mat_values[j];
      }
      ostr() << " index=" << index << '\n';
    }
    info() << ostr.str();
  }
 
  Integer nb_done = 0;
  Integer nb_iter = 0;
  Integer nb_fine = 0;
  Integer nb_coarse = 0;
  m_is_verbose = false;
  {
    // Marque comme point fin tous les points n'ayant aucune dépendance
    for (Integer row = 0; row < nb_row; ++row) {
      if (depends[row].size() == 0) {
        m_points_type[row] = TYPE_FINE;
        ++nb_done;
        if (m_is_verbose)
          info() << "FIRST MARK FINE point=" << row;
      }
    }
 
    // Les points qui n'influencent personne sont forcément fins.
    for (Integer row = 0; row < nb_row; ++row) {
      if (m_points_type[row] != TYPE_FINE && lambdas[row] <= 0) {
        m_points_type[row] = TYPE_FINE;
        ++nb_done;
        if (m_is_verbose)
          info() << "INIT MARK FINE NULL MEASURE point=" << row << " measure=" << lambdas[row];
        for (Integer j = 0, js = depends[row].size(); j < js; ++j) {
          Integer col = depends[row][j];
          if (m_points_type[col] != TYPE_FINE)
            if (col < row) {
              ++lambdas[col];
              if (m_is_verbose)
                printf("ADD MEASURE NULL point=%d measure=%d\n", (int)col, lambdas[col]);
            }
        }
      }
    }
 
    typedef std::set<PointInfo> PointSet;
    PointSet undefined_points;
    for (Integer i = 0; i < nb_row; ++i) {
      if (m_points_type[i] == TYPE_UNDEFINED)
        undefined_points.insert(PointInfo(lambdas[i], i));
    }
 
    while (nb_done < nb_row && nb_iter < 100000) {
      ++nb_iter;
      //for( PointSet::const_iterator i(undefined_points.data()); i!=undefined_points.end(); ++i ){
      //info() << " SET index=" << i->m_index << " value=" << i->m_lambda;
      //}
      // Prend le lambda max et note le point C
      //Integer max_value = -1;
      //Integer max_value_index = -1;
      //for( Integer i=0; i<nb_row; ++i ){
      //if (lambdas[i]>max_value && points_type[i]==TYPE_UNDEFINED){
      //  max_value = lambdas[i];
      //  max_value_index = i;
      //}
      //}
      if (undefined_points.empty())
        fatal() << "Undefined points is empty";
      PointSet::iterator max_point = undefined_points.begin();
      Integer max_value_index = max_point->m_index;
      Integer max_value = max_point->m_lambda;
      m_points_type[max_value_index] = TYPE_COARSE;
      ++nb_done;
      ++nb_coarse;
      undefined_points.erase(max_point);
      if (m_is_verbose)
        std::cout << "MARK COARSE point=" << max_value_index
                  << " measure=" << max_value
                  << " left=" << (nb_row - nb_done)
                  << "\n";
      IntegerConstArrayView point_influences = influences[max_value_index];
      for (Integer i = 0, is = point_influences.size(); i < is; ++i) {
        //for( Integer i=0, is=depends[max_value_index].size(); i<is; ++i ){
        Integer pt = point_influences[i];
        //Integer pt = depends[max_value_index][i];
        if (m_points_type[pt] == TYPE_UNDEFINED) {
          m_points_type[pt] = TYPE_FINE;
          ++nb_done;
          ++nb_fine;
          undefined_points.erase(PointInfo(lambdas[pt], pt));
          if (m_is_verbose)
            std::cout << "MARK FINE point=" << pt
                      << " measure=" << lambdas[pt]
                      << " left=" << (nb_row - nb_done)
                      << "\n";
          for (Integer z = 0, zs = depends[pt].size(); z < zs; ++z) {
            Integer pt2 = depends[pt][z];
            //for( Integer z=0, zs=point_influences.size(); z<zs; ++z ){
            //Integer pt2 = point_influences[z];
            if (m_points_type[pt2] == TYPE_UNDEFINED) {
              undefined_points.erase(PointInfo(lambdas[pt2], pt2));
              ++lambdas[pt2];
              undefined_points.insert(PointInfo(lambdas[pt2], pt2));
            }
          }
        }
      }
      for (Integer i = 0, is = depends[max_value_index].size(); i < is; ++i) {
        Integer pt3 = depends[max_value_index][i];
        if (m_points_type[pt3] == TYPE_UNDEFINED) {
          undefined_points.erase(PointInfo(lambdas[pt3], pt3));
          Integer n = lambdas[pt3];
          if (n < 0)
            info() << "N < 0";
          --lambdas[pt3];
          undefined_points.insert(PointInfo(lambdas[pt3], pt3));
        }
      }
      if (m_is_verbose)
        info() << "LAMBDA MAX = " << max_value << " index=" << max_value_index << " nb_done=" << nb_done;
    }
  }
 
  if (m_is_verbose)
    info() << "NB ROW=" << nb_row << " nb_done=" << nb_done << " nb_fine=" << nb_fine
           << " nb_coarse=" << nb_coarse << " nb_iter=" << nb_iter;
  if (nb_done != nb_row)
    fatal() << "Can not find all COARSE or FINE points nb_done=" << nb_done << " nb_point=" << nb_row;
 
  {
    //Maintenant, nous devons nous assurer que deux connexions F-F ont au moins un
    // point C commun. Sinon, le premier F est changé en C
    //info() << "SECOND PASS !!!";
    Int32UniqueArray points_marker(nb_row);
    points_marker.fill(-1);
    Integer ci_tilde_mark = -1;
    Integer ci_tilde = -1;
    bool C_i_nonempty = false;
    for (Integer row = 0; row < nb_row; ++row) {
      if ((ci_tilde_mark |= row))
        ci_tilde = -1;
      if (m_points_type[row] == TYPE_FINE) {
        for (Integer z = 0, zs = depends[row].size(); z < zs; ++z) {
          //for( Integer z=rows_index[row] ,zs=rows_index[row+1]; z<zs; ++z ){
          //Integer col = columns[z];
          Integer col = depends[row][z];
          if (m_points_type[col] == TYPE_COARSE)
            points_marker[col] = row;
        }
        for (Integer z = 0, zs = depends[row].size(); z < zs; ++z) {
          //for( Integer z=rows_index[row] ,zs=rows_index[row+1]; z<zs; ++z ){
          //Integer col = columns[z];
          Integer col = depends[row][z];
          if (m_points_type[col] == TYPE_FINE) {
            bool set_empty = true;
            for (Integer z2 = 0, zs2 = depends[row].size(); z2 < zs2; ++z2) {
              //for( Integer z2=rows_index[row] ,zs2=rows_index[row+1]; z2<zs2; ++z2 ){
              //Integer col2 = columns[z2];
              Integer col2 = depends[row][z2];
              if (points_marker[col2] == row) {
                set_empty = false;
                break;
              }
            }
            if (set_empty) {
              if (C_i_nonempty) {
                m_points_type[row] = TYPE_COARSE;
                //printf("SECOND PASS MARK COARSE1 point=%d\n",row);
                if (ci_tilde > -1) {
                  m_points_type[ci_tilde] = TYPE_FINE;
                  if (m_is_verbose)
                    printf("SECOND PASS MARK FINE point=%d\n", ci_tilde);
                  ci_tilde = -1;
                }
                C_i_nonempty = false;
              }
              else {
                ci_tilde = col;
                ci_tilde_mark = row;
                m_points_type[col] = TYPE_COARSE;
                if (m_is_verbose)
                  printf("SECOND PASS MARK COARSE2 point=%d\n", col);
                C_i_nonempty = true;
                --row;
                break;
              }
            }
          }
        }
      }
    }
  }
 
  if (0) {
    //Lecture depuis Hypre
    static int matrix_number = 0;
    ++matrix_number;
    info() << "READ HYPRE CF_marker n=" << matrix_number;
    StringBuilder fname("CF_marker-");
    fname += matrix_number;
    std::ifstream ifile(fname.toString().localstr());
    Integer nb_read_point = 0;
    ifile >> std::ws >> nb_read_point >> std::ws;
    if (nb_read_point != nb_row)
      fatal() << "Bad number of points for reading Hypre CF_marker read=" << nb_read_point
              << " expected=" << nb_row << " matrix_number=" << matrix_number;
    nb_coarse = 0;
    nb_fine = 0;
    for (Integer i = 0; i < nb_row; ++i) {
      int pt = 0;
      ifile >> pt;
      if (!ifile)
        fatal() << "Can not read marker point number=" << i;
      if (pt == (-1) || pt == (-3)) {
        m_points_type[i] = TYPE_FINE;
        ++nb_fine;
      }
      else if (pt == 1) {
        m_points_type[i] = TYPE_COARSE;
        ++nb_coarse;
      }
      else
        fatal() << "Bad value read=" << pt << " expected 1 or -1";
    }
  }
 
  // Vérifie que tous les points fins ont au moins un point d'influence
  nb_coarse = 0;
  for (Integer i = 0; i < nb_row; ++i) {
    if (m_points_type[i] == TYPE_UNDEFINED)
      fatal() << " Point " << i << " is undefined";
    if (m_points_type[i] != TYPE_FINE) {
      ++nb_coarse;
      continue;
    }
#if 0
    bool is_ok = false;
    //info() << "CHECK POINT point=" << i
    //       << " depend_size=" << depends[i].size();
    for( Integer z=0, zs=depends[i].size(); z<zs; ++z ){
      if (m_points_type[depends[i][z]]==TYPE_COARSE){
        is_ok = true;
        break;
      }
    }
    //if (!is_ok)
    //  fatal() << " Point " << i << " has no coarse point";
#endif
  }
 
  if (m_is_verbose) {
    OStringStream ostr;
    for (Integer i = 0; i < nb_row; ++i) {
      ostr() << " POINT i=" << i << " type=" << m_points_type[i] << " depends=";
      for (Integer j = 0, js = depends[i].size(); j < js; ++j)
        ostr() << depends[i][j] << ' ';
      ostr() << '\n';
    }
    info() << ostr.str();
  }
 
  nb_fine = nb_row - nb_coarse;
  Integer graph_size = 0;
  for (Integer i = 0; i < nb_row; ++i)
    graph_size += depends[i].size();
 
  info() << " NB COARSE=" << nb_coarse << " NB FINE=" << nb_fine
         << " MAXTRIX NON_ZEROS=" << m_fine_matrix.rowsIndex()[nb_row]
         << " GRAPH_SIZE=" << graph_size;
  bool dump_matrix = false;
  bool has_error = false;
  if (nb_fine == 0 || graph_size == 0) {
    has_error = true;
    dump_matrix = true;
  }
 
  if (dump_matrix) {
    OStringStream ostr;
    Integer index = 0;
    int n = math::min(nb_row, 40);
    if (0) {
      ostr() << "GRAPH\n";
      for (Integer i = 0; i < n; ++i) {
        ostr() << " GRAPH I=" << i << " ";
        for (Integer j = 0; j < depends[i].size(); ++j) {
          ++index;
          ostr() << " " << depends[i][j];
        }
        ostr() << " index=" << index << '\n';
      }
    }
    ostr() << "\n MAXTRIX\n";
    index = 0;
    for (Integer i = 0; i < n; ++i) {
      ostr() << "MATRIX I=" << i << " ";
      for (Integer j = rows_index[i]; j < rows_index[i + 1]; ++j) {
        ++index;
        ostr() << " " << columns[j] << " " << mat_values[j];
      }
      ostr() << " index=" << index << '\n';
    }
    info() << ostr.str();
  }
  if (has_error)
    throw FatalErrorException("AMGLevel::_buildCoarsePoints");
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMGLevel::
buildLevel(Matrix matrix, Real alpha)
{
  //Integer nb_row = matrix.nbRow();
  //if (nb_row<20)
  //return;
 
  m_fine_matrix = matrix;
 
  //bool is_verbose = false;
  matrix.sortDiagonale();
 
  IntegerUniqueArray points_type;
  RealUniqueArray rows_max_val;
  UniqueArray<SharedArray<Integer>> depends;
  IntegerUniqueArray weak_depends;
 
  _buildCoarsePoints(alpha, rows_max_val, depends, weak_depends);
  _buildInterpolationMatrix(rows_max_val, depends, weak_depends);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMGLevel::
_buildInterpolationMatrix(RealConstArrayView rows_max_val,
                          UniqueArray<SharedArray<Integer>>& depends,
                          IntegerArray& weak_depends)
{
  ARCANE_UNUSED(rows_max_val);
 
  IntegerConstArrayView rows_index = m_fine_matrix.rowsIndex();
  IntegerConstArrayView columns = m_fine_matrix.columns();
  RealConstArrayView mat_values = m_fine_matrix.values();
  Integer nb_row = m_fine_matrix.nbRow();
 
  IntegerUniqueArray points_in_coarse(nb_row);
  points_in_coarse.fill(-1);
  Integer nb_coarse = 0;
  {
    //Integer index = 0;
    nb_coarse = 0;
    for (Integer i = 0; i < nb_row; ++i) {
      if (m_points_type[i] == TYPE_COARSE) {
        points_in_coarse[i] = nb_coarse;
        ++nb_coarse;
      }
    }
  }
  bool type_hypre = true;
 
  // Maintenant,calcule les éléments de la matrice d'influence
  IntegerUniqueArray prolongation_matrix_columns;
  RealUniqueArray prolongation_matrix_values;
  IntegerUniqueArray prolongation_matrix_rows_size(nb_row);
 
  for (Integer row = 0; row < nb_row; ++row) {
    Integer nb_column = 0;
    if (m_points_type[row] == TYPE_FINE) {
      Real weak_connect_sum = 0.0;
      //Real max_value = rows_max_val[row];
      Real diag = mat_values[rows_index[row]];
      Real sign = 1.0;
      if (diag < 0.0)
        sign = -1.0;
      for (Integer z = rows_index[row] + 1, zs = rows_index[row + 1]; z < zs; ++z) {
        //Integer column = columns[z];
        if (weak_depends[z] == 1) {
          Real mv = mat_values[z];
          if (type_hypre)
            weak_connect_sum += mv;
          else {
            //weak_connect_sum += math::abs(mv);
            weak_connect_sum += mv;
          }
          //if (m_is_verbose || row<=5)
          //info() << "ADD WEAK_SUM mv=" << mv << " sum=" << weak_connect_sum
          //         << " row=" << row << " column=" << column;
        }
      }
      if (m_is_verbose)
        info() << "ROW row=" << row << " weak_connect_sum=" << weak_connect_sum;
      //for( Integer z=0, zs=depends[row].size(); z<zs; ++ z){
      //Integer j_column = dependcolumns[z];
      for (Integer z = rows_index[row] + 1, zs = rows_index[row + 1]; z < zs; ++z) {
        Integer j_column = columns[z];
        if (m_points_type[j_column] != TYPE_COARSE)
          continue;
        if (weak_depends[z] != 2)
          continue;
        Real num_add = 0.0;
        Real mv = mat_values[z];
        for (Integer z2 = 0, zs2 = depends[row].size(); z2 < zs2; ++z2) {
          Integer k_column = depends[row][z2];
          //if (m_is_verbose || row<=5)
          //info() << "CHECK K_COLUMN row=" << row << " col=" << k_column << " val=" << mv;
          if (m_points_type[k_column] != TYPE_FINE)
            continue;
          Real sum_coarse = 0.0;
          //if (m_is_verbose || row<=5)
          //info() << "CHECK COARSE row=" << row;
          //for( Integer z3=rows_index[row]+1 ,zs3=rows_index[row+1]; z3<zs3; ++z3 ){
          //Integer m_column = columns[z3];
          for (Integer z3 = 0, zs3 = depends[row].size(); z3 < zs3; ++z3) {
            Integer m_column = depends[row][z3];
            //if (m_is_verbose || row<=5)
            //info() << "CHECK COLUMN column=" << m_column;
            if (m_points_type[m_column] == TYPE_COARSE) {
              Real w = m_fine_matrix.value(k_column, m_column);
              //if (m_is_verbose || row<=5)
              //info() << "ADD SUM k=" << k_column << " m="<< m_column << " w=" << w;
              //if (math::isZero(w)){
              //fatal() << "WEIGHT is null k=" << k_column << " m=" << m_column << " row=" << row
              //        << " j=" << j_column;
              //}
              if (type_hypre) {
                if (w * sign < 0.0)
                  sum_coarse += w;
              }
              else {
                sum_coarse += math::abs(w);
                //if (w*sign<0.0)
                //sum_coarse += w;
              }
            }
          }
          Real to_add = 0.0;
          if (!math::isZero(sum_coarse)) {
            Real akj = m_fine_matrix.value(k_column, j_column);
            bool do_add = false;
            if (type_hypre) {
              if ((akj * sign) < 0.0)
                do_add = true;
            }
            else {
              //if ((akj*sign)<0.0)
              //do_add = true;
              akj = math::abs(akj);
              do_add = true;
            }
            if (do_add)
              //fatal() << "SUM_WEIGHT is null k=" << k_column << " row=" << row << " j=" << j_column;
              to_add = math::divide(m_fine_matrix.value(row, k_column) * akj, sum_coarse);
          }
          num_add += to_add;
        }
        Real weight = -(mv + num_add) / (diag + weak_connect_sum);
        Integer new_column = points_in_coarse[j_column];
        //if (m_is_verbose || row<=5)
        //info() << " ** WEIGHT row=" << row << " j_column=" << j_column
        //         << " weight=" << weight << " num_add=" << num_add << " mv=" << mv << " new_column=" << new_column
        //         << " diag=" << diag << " weak_sum=" << weak_connect_sum
        //         << " diag+wk=" << (diag+weak_connect_sum);
        if (new_column >= nb_coarse || new_column < 0)
          fatal() << " BAD COLUMN for fine point column=" << new_column << " nb=" << nb_coarse
                  << " jcolumn=" << j_column;
        prolongation_matrix_columns.add(new_column);
        prolongation_matrix_values.add(weight);
        ++nb_column;
      }
    }
    else {
      // Point grossier, met 1.0 dans la diagonale
      Integer column = points_in_coarse[row];
      if (column >= nb_coarse || column < 0)
        fatal() << " BAD COLUMN for coarse point j=" << column << " nb=" << nb_coarse
                << " row=" << row;
      prolongation_matrix_columns.add(column);
      prolongation_matrix_values.add(1.0);
      ++nb_column;
    }
    prolongation_matrix_rows_size[row] = nb_column;
  }
 
  m_prolongation_matrix = Matrix(nb_row, nb_coarse);
  m_prolongation_matrix.setRowsSize(prolongation_matrix_rows_size);
  //info() << "PROLONGATION_MATRIX_SIZE=" << m_prolongation_matrix.rowsIndex()[nb_row];
  m_prolongation_matrix.setValues(prolongation_matrix_columns, prolongation_matrix_values);
 
  if (0) {
    OStringStream ostr;
    Integer index = 0;
    int n = math::min(nb_row, 50);
    IntegerConstArrayView p_rows(m_prolongation_matrix.rowsIndex());
    IntegerConstArrayView p_columns(m_prolongation_matrix.columns());
    RealConstArrayView p_values(m_prolongation_matrix.values());
    for (Integer i = 0; i < n; ++i) {
      ostr() << "PROLONG I=" << i << " ";
      for (Integer j = p_rows[i]; j < p_rows[i + 1]; ++j) {
        ++index;
        ostr() << " " << p_columns[j] << " " << p_values[j];
      }
      ostr() << " index=" << index << '\n';
    }
    info() << "PROLONG\n"
           << ostr.str();
  }
 
  MatrixOperation2 mat_op2;
  if (1)
    m_restriction_matrix = mat_op2.transposeFast(m_prolongation_matrix);
  else
    m_restriction_matrix = mat_op2.transpose(m_prolongation_matrix);
  if (m_is_verbose) {
    OStringStream ostr;
    ostr() << "PROLONGATION_MATRIX ";
    m_prolongation_matrix.dump(ostr());
    ostr() << '\n';
    ostr() << "RESTRICTION_MATRIX ";
    m_restriction_matrix.dump(ostr());
    info() << ostr.str();
  }
  //info() << " ** TOTAL SUM=" << total_sum;
  // Calcule la matrice grossiere Ak+1 = I * Ak * tI
  //MatrixOperation mat_op;
  const bool old = false;
  if (old) {
    Matrix n1 = mat_op2.matrixMatrixProductFast(m_fine_matrix, m_prolongation_matrix);
    if (m_is_verbose) {
      OStringStream ostr;
      n1.dump(ostr());
      info() << "N1_MATRIX " << ostr.str();
    }
    m_coarse_matrix = mat_op2.matrixMatrixProductFast(m_restriction_matrix, n1);
  }
  else
    m_coarse_matrix = mat_op2.applyGalerkinOperator2(m_restriction_matrix, m_fine_matrix, m_prolongation_matrix);
  if (m_is_verbose) {
    OStringStream ostr;
    m_coarse_matrix.dump(ostr());
    info() << "level= " << m_level << " COARSE_MATRIX=" << ostr.str();
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
AMGPreconditioner::
~AMGPreconditioner()
{
  delete m_amg;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMGPreconditioner::
apply(Vector& out_vec, const Vector& vec)
{
  m_amg->solve(vec, out_vec);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMGPreconditioner::
build(const Matrix& matrix)
{
  delete m_amg;
  m_amg = new AMG(m_trace_mng);
  m_amg->build(matrix);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
AMGSolver::
~AMGSolver()
{
  delete m_amg;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMGSolver::
build(const Matrix& matrix)
{
  delete m_amg;
  m_amg = new AMG(m_trace_mng);
  m_amg->build(matrix);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void AMGSolver::
solve(const Vector& vector_b, Vector& vector_x)
{
  m_amg->solve(vector_b, vector_x);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
} // namespace Arcane::MatVec
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/