Matrix.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
//-----------------------------------------------------------------------------
/*---------------------------------------------------------------------------*/
/* Matrix.cc                                                   (C) 2000-2026 */
/*                                                                           */
/* Linear algebra matrix.                                                    */
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
#include "arcane/utils/Array.h"
#include "arcane/utils/FatalErrorException.h"
#include "arcane/utils/Iostream.h"
#include "arcane/utils/Numeric.h"
#include "arcane/utils/ArgumentException.h"
#include "arcane/utils/NotImplementedException.h"
#include "arcane/utils/TraceInfo.h"
#include "arcane/utils/StringBuilder.h"
#include "arcane/utils/ValueConvert.h"
 
#include "arcane/core/matvec/Matrix.h"
#include "arcane/core/matvec/Vector.h"
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
namespace Arcane::MatVec
{
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

class MatrixImpl
{
 public:
 
  MatrixImpl()
  : m_nb_reference(0)
  , m_nb_row(0)
  , m_nb_column(0)
  , m_nb_element(0)
  {}
  MatrixImpl(Integer nb_row, Integer nb_column);
 
 private:
 
  void operator=(const Matrix& rhs);
 
 public:
 
  MatrixImpl* clone()
  {
    MatrixImpl* m2 = new MatrixImpl();
    m2->m_nb_row = m_nb_row;
    m2->m_nb_column = m_nb_column;
    m2->m_nb_element = m_nb_element;
    m2->m_values = m_values.clone();
    m2->m_rows_index = m_rows_index.clone();
    m2->m_columns = m_columns.clone();
    return m2;
  }
 
 private:
 public:
 
  Integer nbRow() const { return m_nb_row; }
  Integer nbColumn() const { return m_nb_column; }
  void setRowsSize(IntegerConstArrayView rows_size);
  void setValues(IntegerConstArrayView columns, RealConstArrayView values);
  void dump(std::ostream& o);
  RealConstArrayView values() const { return m_values; }
  RealArrayView values() { return m_values; }
  Real value(Integer row, Integer column) const;
  IntegerConstArrayView rowsIndex() const { return m_rows_index; }
  IntegerConstArrayView columns() const { return m_columns; }
  IntegerArrayView rowsIndex() { return m_rows_index; }
  IntegerArrayView columns() { return m_columns; }
  void setValue(Integer row, Integer column, Real value);
  void sortDiagonale();
  void assemble();
 
 public:
 
  void checkValid();
 
 public:
 
  void addReference()
  {
    ++m_nb_reference;
  }
  void removeReference()
  {
    --m_nb_reference;
    if (m_nb_reference == 0)
      delete this;
  }
 
 private:
 
  Integer m_nb_reference;
  Integer m_nb_row;
  Integer m_nb_column;
  Integer m_nb_element;
  SharedArray<Real> m_values;
  SharedArray<Integer> m_rows_index;
  SharedArray<Integer> m_columns;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
MatrixImpl::
MatrixImpl(Integer nb_row, Integer nb_column)
: m_nb_reference(0)
, m_nb_row(nb_row)
, m_nb_column(nb_column)
, m_nb_element(0)
, m_rows_index(m_nb_row + 1)
{
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix::
Matrix(Integer nb_row, Integer nb_column)
: m_impl(new MatrixImpl(nb_row, nb_column))
{
  m_impl->addReference();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix::
Matrix(MatrixImpl* impl)
: m_impl(impl)
{
  if (m_impl)
    m_impl->addReference();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix::
Matrix(const Matrix& rhs)
: m_impl(rhs.m_impl)
{
  if (m_impl)
    m_impl->addReference();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void Matrix::
operator=(const Matrix& rhs)
{
  if (rhs.m_impl)
    rhs.m_impl->addReference();
  if (m_impl)
    m_impl->removeReference();
  m_impl = rhs.m_impl;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix::
~Matrix()
{
  if (m_impl)
    m_impl->removeReference();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Integer Matrix::
nbRow() const
{
  return m_impl->nbRow();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void Matrix::
sortDiagonale()
{
  return m_impl->sortDiagonale();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Integer Matrix::
nbColumn() const
{
  return m_impl->nbColumn();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix Matrix::
clone() const
{
  MatrixImpl* m = m_impl->clone();
  return Matrix(m);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
IntegerConstArrayView Matrix::
rowsIndex() const
{
  return m_impl->rowsIndex();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void Matrix::
dump(std::ostream& o) const
{
  m_impl->dump(o);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void Matrix::
setValue(Integer row, Integer column, Real value)
{
  m_impl->setValue(row, column, value);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Real Matrix::
value(Integer row, Integer column) const
{
  return m_impl->value(row, column);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void Matrix::
setValues(IntegerConstArrayView columns, RealConstArrayView values)
{
  m_impl->setValues(columns, values);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void Matrix::
setRowsSize(IntegerConstArrayView rows_size)
{
  m_impl->setRowsSize(rows_size);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
RealConstArrayView Matrix::
values() const
{
  return m_impl->values();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
RealArrayView Matrix::
values()
{
  return m_impl->values();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
IntegerArrayView Matrix::
columns()
{
  return m_impl->columns();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
IntegerArrayView Matrix::
rowsIndex()
{
  return m_impl->rowsIndex();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
IntegerConstArrayView Matrix::
columns() const
{
  return m_impl->columns();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
DiagonalPreconditioner::
DiagonalPreconditioner(const Matrix& matrix)
: m_inverse_diagonal(matrix.nbRow())
{
  Integer size = m_inverse_diagonal.size();
  IntegerConstArrayView rows_index = matrix.rowsIndex();
  IntegerConstArrayView columns = matrix.columns();
  RealConstArrayView mat_values = matrix.values();
  RealArrayView vec_values = m_inverse_diagonal.values();
  for (Integer i = 0, is = size; i < is; ++i) {
    for (Integer z = rows_index[i], zs = rows_index[i + 1]; z < zs; ++z) {
      Integer mj = columns[z];
      if (mj == i)
        vec_values[i] = mat_values[z];
    }
  }
#if 0
  cout << "DIAG1=";
  m_inverse_diagonal.dump(cout);
  cout << "\n";
#endif
  const double epsilon = std::numeric_limits<Real>::min();
 
  // Invert the diagonal if the value is not less than the Real epsilon
  // (otherwise it generates an FPE)
  for (Integer i = 0; i < size; ++i) {
    Real v = vec_values[i];
    bool is_zero = v > -epsilon && v < epsilon;
    vec_values[i] = (is_zero) ? 1.0 : (1.0 / v);
  }
 
#if 0
  cout << "DIAG2=";
  m_inverse_diagonal.dump(cout);
  cout << "\n";
#endif
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void DiagonalPreconditioner::
apply(Vector& out_vec, const Vector& vec)
{
  Integer size = m_inverse_diagonal.size();
  RealConstArrayView inverse_diagonal_values = m_inverse_diagonal.values();
  RealConstArrayView vec_values = vec.values();
  RealArrayView out_vec_values = out_vec.values();
  for (Integer i = 0; i < size; ++i)
    out_vec_values[i] = vec_values[i] * inverse_diagonal_values[i];
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixOperation::
matrixVectorProduct(const Matrix& mat, const Vector& vec, Vector& out_vec)
{
  Integer nb_row = mat.nbRow();
  Integer nb_column = mat.nbColumn();
  //cout << "** MATRIX ROW=" << nb_row << " COL=" << nb_column
  //     << " intput_size=" << vec.size()
  //     << " output_size=" << out_vec.size()
  //     << "\n";
  if (nb_column != vec.size())
    throw ArgumentException("MatrixVectorProduct", "Bad size for input vector");
  if (nb_row != out_vec.size())
    throw ArgumentException("MatrixVectorProduct", "Bad size for output_vector");
  IntegerConstArrayView rows_index = mat.rowsIndex();
  IntegerConstArrayView columns = mat.columns();
  RealConstArrayView mat_values = mat.values();
  RealConstArrayView vec_values = vec.values();
  RealArrayView out_vec_values = out_vec.values();
 
  for (Integer i = 0, is = nb_row; i < is; ++i) {
    Real sum = 0.0;
    for (Integer z = rows_index[i], zs = rows_index[i + 1]; z < zs; ++z) {
      Integer mj = columns[z];
      Real mv = mat_values[z];
      //cout << "ADD vec=" << vec_values[mj] << " mv=" << mv << " sum=" << sum << '\n';
      sum += vec_values[mj] * mv;
    }
    out_vec_values[i] = sum;
  }
}
 
Real MatrixOperation::
dot(const Vector& vec)
{
  Integer size = vec.size();
  RealConstArrayView vec_values = vec.values();
  Real v = 0.0;
  for (Integer i = 0; i < size; ++i)
    v += vec_values[i] * vec_values[i];
  return v;
}
 
Real MatrixOperation::
dot(const Vector& vec1, const Vector& vec2)
{
  Real v = 0.0;
  Integer size = vec1.size();
  RealConstArrayView vec1_values = vec1.values();
  RealConstArrayView vec2_values = vec2.values();
  for (Integer i = 0; i < size; ++i) {
    v += vec1_values[i] * vec2_values[i];
    //cout << " i=" << i << " v=" << v << '\n';
  }
  return v;
}
 
void MatrixOperation::
negateVector(Vector& vec)
{
  Integer size = vec.size();
  RealArrayView vec_values = vec.values();
  for (Integer i = 0; i < size; ++i)
    vec_values[i] = -vec_values[i];
}
 
void MatrixOperation::
scaleVector(Vector& vec, Real mul)
{
  Integer size = vec.size();
  RealArrayView vec_values = vec.values();
  for (Integer i = 0; i < size; ++i)
    vec_values[i] *= mul;
}
 
void MatrixOperation::
addVector(Vector& out_vec, const Vector& vec)
{
  Integer size = vec.size();
  RealConstArrayView vec_values = vec.values();
  RealArrayView out_vec_values = out_vec.values();
  for (Integer i = 0; i < size; ++i)
    out_vec_values[i] += vec_values[i];
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void ConjugateGradientSolver::
_applySolver(const Matrix& a, const Vector& b, Vector& x, Real epsilon, IPreconditioner* p)
{
  MatrixOperation mat_op;
  Integer vec_size = a.nbRow();
  Vector r(vec_size);
 
  // r = b - Ax
  mat_op.matrixVectorProduct(a, x, r);
  mat_op.negateVector(r);
  mat_op.addVector(r, b);
 
  m_nb_iteration = 0;
  m_residual_norm = 0.0;
 
  /*cout << " R=";
  r.dump(cout);
  cout << '\n';*/
 
  Vector d(r.size(), 0.0);
  if (p)
    p->apply(d, r);
  else
    d.copy(r);
 
  Vector q(r.size());
  Vector t(r.size());
  Vector s(r.size(), 0.0);
  Real delta_new = 0.0;
  //Real r0=mat_op.dot(r);
  if (p) {
    delta_new = mat_op.dot(r, d);
  }
  else
    delta_new = mat_op.dot(r);
#if 0
  cout << "R=";
  r.dump(cout);
  cout << "\n";
  cout << "D=";
  d.dump(cout);
  cout << "\n";
#endif
  //Real norm0 = r.normInf();
  const Real delta0 = delta_new;
 
  //cout << " TOL=" << epsilon << " delta0=" << delta0 << '\n';
  //cout << " deltanew=" << delta_new << '\n';
  Integer nb_iter = 0;
  for (nb_iter = 0; nb_iter < m_max_iteration; ++nb_iter) {
    if (delta_new < epsilon * epsilon * delta0)
      break;
#if 0
    // If we use the infinity norm
    {
      Real norm = r.normInf();
      cout << " norm=" << norm << " norm0=" << norm0 << '\n';
      if (norm < epsilon * norm0)
        break;
    }
#endif
 
    //cout << "delta_new=" << delta_new << '\n';
    // q = A * d
    mat_op.matrixVectorProduct(a, d, q);
    Real alpha = delta_new / mat_op.dot(d, q);
 
    // x = x + alpha * d
    t.copy(d);
    mat_op.scaleVector(t, alpha);
    mat_op.addVector(x, t);
 
#if 0
    // r = b - Ax
    mat_op.matrixVectorProduct(a,x,r);
    mat_op.negateVector(r);
    mat_op.addVector(r,b);
#endif
    // r = r - alpha * q
    mat_op.scaleVector(q, -alpha);
    mat_op.addVector(r, q);
 
    if (p)
      p->apply(s, r);
    Real delta_old = delta_new;
    if (p)
      delta_new = mat_op.dot(r, s);
    else
      delta_new = mat_op.dot(r);
    Real beta = delta_new / delta_old;
    //cout << " alpha=" << alpha << " beta=" << beta << " delta_new=" << delta_new << '\n';
    // d = beta * d + s
    mat_op.scaleVector(d, beta);
    if (p) {
      //mat_op.addVector(s,r);
      mat_op.addVector(d, s);
    }
    else
      mat_op.addVector(d, r);
    //cout << '\n';
  }
  //cout << " X=";
  //x.dump(cout);
  //cout << '\n';
  //cout << "NB ITER=" << nb_iter << " epsilon=" << epsilon
  //     << " delta0=" << delta0
  //     << " delta_new=" << delta_new << " r=" << mat_op.dot(r) << " r0=" << r0 << '\n';
  m_nb_iteration = nb_iter;
  m_residual_norm = delta_new;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void ConjugateGradientSolver::
_applySolverAsHypre(const Matrix& a, const Vector& b, Vector& x, Real tol,
                    IPreconditioner* precond)
{
  //cout << "** APPLY ARCANE_SOLVER PCG using hypre method\n";
  ARCANE_CHECK_POINTER(precond);
 
  MatrixOperation mat_op;
  Integer vec_size = a.nbRow();
  Vector r(vec_size);
  Vector p(vec_size);
 
  const bool is_two_norm = true;
  Real bi_prod = 0.0;
  if (is_two_norm) {
    bi_prod = mat_op.dot(b, b);
  }
  else {
    precond->apply(p, b);
    bi_prod = mat_op.dot(p, b);
  }
  Real eps = tol * tol;
  //TODO: regarder stop_crit
 
  // r = b - Ax
  mat_op.matrixVectorProduct(a, x, r);
  mat_op.negateVector(r);
  mat_op.addVector(r, b);
 
  m_nb_iteration = 0;
  m_residual_norm = 0.0;
 
  /*cout << " R=";
  r.dump(cout);
  cout << '\n';*/
 
  Vector tmp_x(x.size());
 
  Vector d(r.size());
  Vector s(r.size());
  d.copy(r);
  // p = C*r
  precond->apply(p, r);
 
  /* gamma = <r,p> */
  Real gamma = mat_op.dot(r, p);
 
  Real i_prod = 0.0;
 
  //Real delta_new = 0.0;
  //Real r0=mat_op.dot(r);
  //if (p){
  //  delta_new = mat_op.dot(r,d);
  //}
  //else
  // delta_new = mat_op.dot(r);
#if 0
  cout << "R=";
  r.dump(cout);
  cout << "\n";
  cout << "D=";
  d.dump(cout);
  cout << "\n";
#endif
  //Real delta0 = delta_new;
  //cout << "** TOL=" << tol << " bi_prod=" << bi_prod << '\n';
  //cout << " deltanew=" << delta_new << '\n';
  Integer nb_iter = 0;
  for (nb_iter = 0; nb_iter < m_max_iteration; ++nb_iter) {
 
    // s = A * p
    mat_op.matrixVectorProduct(a, p, s);
 
    /* alpha = gamma / <s,p> */
    Real sdotp = mat_op.dot(s, p);
 
    Real alpha = gamma / sdotp;
 
    Real gamma_old = gamma;
 
    /* x = x + alpha*p */
    mat_op.matrixVectorProduct(a, p, tmp_x);
    mat_op.addVector(x, tmp_x);
 
    /* r = r - alpha*s */
    tmp_x.copy(s);
    mat_op.scaleVector(tmp_x, -alpha);
    mat_op.addVector(r, tmp_x);
 
    /* s = C*r */
    precond->apply(s, r);
 
    /* gamma = <r,s> */
    gamma = mat_op.dot(r, s);
 
    /* set i_prod for convergence test */
    if (is_two_norm)
      i_prod = mat_op.dot(r, r);
    else
      i_prod = gamma;
 
    //cout << "** sdotp=" << sdotp << " i_prod=" << i_prod << " bi_prod=" << bi_prod << '\n';
 
    /* check for convergence */
    if (i_prod / bi_prod < eps) {
#if 0
      if (rel_change && i_prod > guard_zero_residual)
      {
        pi_prod = (*(pcg_functions->InnerProd))(p,p); 
        xi_prod = (*(pcg_functions->InnerProd))(x,x);
        ratio = alpha*alpha*pi_prod/xi_prod;
        if (ratio < eps)
        {
          (pcg_data -> converged) = 1;
          break;
        }
      }
      else
      {
#endif
      //(pcg_data -> converged) = 1;
      break;
#if 0
      }
#endif
    }
 
    /* beta = gamma / gamma_old */
    Real beta = gamma / gamma_old;
 
    //cout << " alpha=" << alpha << " beta=" << beta << " delta_new=" << gamma << '\n';
 
    /* p = s + beta p */
    tmp_x.copy(p);
    mat_op.scaleVector(tmp_x, beta);
    mat_op.addVector(tmp_x, s);
    p.copy(tmp_x);
  }
  //cout << " X=";
  //x.dump(cout);
  //cout << '\n';
  //cout << "NB ITER=" << nb_iter << " epsilon=" << epsilon
  //     << " delta0=" << delta0
  //     << " delta_new=" << delta_new << " r=" << mat_op.dot(r) << " r0=" << r0 << '\n';
  m_nb_iteration = nb_iter;
  m_residual_norm = gamma;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
bool ConjugateGradientSolver::
solve(const Matrix& a, const Vector& b, Vector& x, Real epsilon, IPreconditioner* p)
{
  //epsilon = 1e-12;
  //x.values().fill(0.0);
  m_nb_iteration = 0;
  m_residual_norm = 0.0;
  //_doConjugateGradient(a,b,x,epsilon,&p);
  cout.precision(20);
  const bool do_print = false;
  if (do_print) {
    cout << "A=";
    a.dump(cout);
    cout << '\n';
    cout << "b=";
    b.dump(cout);
    cout << '\n';
  }
 
  DiagonalPreconditioner p2(a);
  if (!p)
    p = &p2;
  _applySolver(a, b, x, epsilon, p);
 
  if (do_print) {
    cout << "\n\nSOLUTION\nx=";
    x.dump(cout);
    cout << '\n';
  }
 
  return false;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void
_testArcaneMatrix1()
{
  Integer s = 5;
  Matrix m(s, s);
  Matrix m1(Matrix::read("test.mat"));
  cout << "** o=" << m.nbRow() << '\n';
  cout << "** M1=";
  m1.dump(cout);
  cout << '\n';
  Vector v1(10);
  for (Integer i = 0; i < 10; ++i) {
    v1.values()[i] = (Real)(i + 1);
  }
  Real epsilon = 1.0e-10;
  {
    Vector v2(10);
    MatrixOperation mat_op;
    Vector r(5);
    mat_op.matrixVectorProduct(m1, v1, v2);
    cout << "** V1=";
    v1.dump(cout);
    cout << "\n";
    cout << "** V2=";
    v2.dump(cout);
    cout << "\n";
    Vector b3(10);
    for (Integer i = 0; i < 10; ++i) {
      b3.values()[i] = (Real)(i + 1);
    }
    Vector x3(10);
    DiagonalPreconditioner p(m1);
    ConjugateGradientSolver solver;
    solver.solve(m1, b3, x3, epsilon, &p);
  }
 
  IntegerUniqueArray rows_size(s);
  rows_size[0] = 5;
  rows_size[1] = 2;
  rows_size[2] = 2;
  rows_size[3] = 2;
  rows_size[4] = 2;
  m.setRowsSize(rows_size);
  RealUniqueArray values(13);
  IntegerUniqueArray columns(13);
  values[0] = 9.0;
  values[1] = 1.5;
  values[2] = 6.0;
  values[3] = 0.75;
  values[4] = 3.0;
 
  values[5] = 1.5;
  values[6] = 0.5;
 
  values[7] = 6.0;
  values[8] = 0.5;
 
  values[9] = 0.75;
  values[10] = 5.0 / 8.0;
 
  values[11] = 3.0;
  values[12] = 16.0;
 
  columns[0] = 0;
  columns[1] = 1;
  columns[2] = 2;
  columns[3] = 3;
  columns[4] = 4;
  columns[5] = 0;
  columns[6] = 1;
  columns[7] = 0;
  columns[8] = 2;
  columns[9] = 0;
  columns[10] = 3;
  columns[11] = 0;
  columns[12] = 4;
  m.setValues(columns, values);
  m.dump(cout);
  cout << '\n';
  Vector b(5);
  RealArrayView rav(b.values());
  rav[0] = 1.0;
  rav[1] = 1.0;
  rav[2] = 1.0;
  rav[3] = 3.0;
  rav[4] = 1.0;
  Vector x(5);
  ConjugateGradientSolver solver;
  solver.solve(m, b, x, epsilon);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void _testArcaneMatrix2(Integer matrix_number)
{
  cout.precision(16);
  cout.flags(std::ios::scientific);
  String dir_name = "test-matrix";
  StringBuilder ext_name;
  ext_name += matrix_number;
  ext_name += ".00000";
  Matrix m(Matrix::readHypre(dir_name + "/MATRIX_matrix" + ext_name.toString()));
  Vector b(Vector::readHypre(dir_name + "/MATRIX_b" + ext_name.toString()));
  Vector xref(Vector::readHypre(dir_name + "/MATRIX_x" + ext_name.toString()));
  Vector x(m.nbRow());
  cout << "** XREF=" << xref.size() << '\n';
  //m.dump(cout);
  cout << '\n';
  {
    Real epsilon = 1.0e-15;
    {
      ConjugateGradientSolver solver;
      solver.solve(m, b, x, epsilon);
    }
    MatrixOperation mat_op;
    mat_op.negateVector(xref);
    mat_op.addVector(xref, x);
    //xref.dump(cout);
    Real x_norm = mat_op.dot(x);
    Real diff_norm = mat_op.dot(xref);
    cout << "** MATRIX_NUMBER = " << matrix_number;
    if (!math::isNearlyZero(x_norm)) {
      cout << "** norm=" << x_norm << " REL=" << diff_norm / x_norm << '\n';
    }
    else
      cout << "** norm=" << x_norm << " DIFF=" << diff_norm << '\n';
    cout << '\n';
  }
  {
    Real epsilon = 1.0e-15;
    x.values().fill(0.0);
    {
      DiagonalPreconditioner p(m);
      ConjugateGradientSolver solver;
      solver.solve(m, b, x, epsilon, &p);
    }
    MatrixOperation mat_op;
    mat_op.negateVector(xref);
    mat_op.addVector(xref, x);
    //xref.dump(cout);
    Real x_norm = mat_op.dot(x);
    Real diff_norm = mat_op.dot(xref);
    cout << "** PRECOND MATRIX_NUMBER = " << matrix_number;
    if (!math::isNearlyZero(x_norm)) {
      cout << "** norm=" << x_norm << " REL=" << diff_norm / x_norm << '\n';
    }
    else
      cout << "** norm=" << x_norm << " DIFF=" << diff_norm << '\n';
    cout << '\n';
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
extern "C" void _testArcaneMatrix()
{
  _testArcaneMatrix1();
  //for( Integer i=1; i<15; ++i )
  //_testArcaneMatrix2(i);
  exit(0);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixImpl::
setRowsSize(IntegerConstArrayView rows_size)
{
  Integer index = 0;
  for (Integer i = 0, is = m_nb_row; i < is; ++i) {
    m_rows_index[i] = index;
    index += rows_size[i];
  }
  m_rows_index[m_nb_row] = index;
  m_nb_element = index;
  m_columns.resize(index);
  m_columns.fill(-1);
  m_values.resize(index);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixImpl::
setValues(IntegerConstArrayView columns, RealConstArrayView values)
{
  m_columns.copy(columns);
  m_values.copy(values);
  if (arcaneIsCheck())
    checkValid();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixImpl::
checkValid()
{
  IntegerConstArrayView columns = m_columns;
  IntegerConstArrayView rows_index = m_rows_index;
 
  Integer nb_column = m_nb_column;
  for (Integer row = 0, nb_row = m_nb_row; row < nb_row; ++row) {
    for (Integer j = rows_index[row], js = rows_index[row + 1]; j < js; ++j) {
      if (columns[j] >= nb_column || columns[j] < 0) {
        cout << "BAD COLUMN VALUE for row=" << row << " column=" << columns[j]
             << " column_index=" << j
             << " nb_column=" << nb_column << " nb_row=" << nb_row << '\n';
        throw FatalErrorException("MatrixImpl::checkValid");
      }
    }
  }
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Real MatrixImpl::
value(Integer row, Integer column) const
{
  IntegerConstArrayView rows_index = m_rows_index;
  IntegerConstArrayView columns = m_columns;
  RealConstArrayView values = m_values;
 
  for (Integer z = rows_index[row], zs = rows_index[row + 1]; z < zs; ++z) {
    if (columns[z] == column)
      return values[z];
  }
  return 0.0;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixImpl::
setValue(Integer row, Integer column, Real value)
{
  IntegerConstArrayView rows_index = m_rows_index;
  IntegerArrayView columns = m_columns;
#ifdef ARCANE_CHECK
  if (row >= m_nb_row)
    throw ArgumentException("MatrixImpl::setValue", "Invalid row");
  if (column >= m_nb_column)
    throw ArgumentException("MatrixImpl::setValue", "Invalid column");
  if (row < 0)
    throw ArgumentException("MatrixImpl::setValue", "Invalid row");
  if (column < 0)
    throw ArgumentException("MatrixImpl::setValue", "Invalid column");
#endif
  for (Integer j = rows_index[row], js = rows_index[row + 1]; j < js; ++j) {
    if (columns[j] == (-1) || columns[j] == column) {
      columns[j] = column;
      m_values[j] = value;
      return;
    }
  }
  throw ArgumentException("Matrix::setValue", "column not found");
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixImpl::
sortDiagonale()
{
  assemble();
  if (arcaneIsCheck()) {
    checkValid();
  }
  // Sort the diagonal so that it is the first element of the row
  IntegerConstArrayView rows_index = m_rows_index;
  IntegerArrayView columns = m_columns;
  RealArrayView values = m_values;
  for (Integer row = 0, nb_row = m_nb_row; row < nb_row; ++row) {
    Integer first_col = rows_index[row];
    for (Integer j = first_col, js = rows_index[row + 1]; j < js; ++j) {
      if (columns[j] == row) {
        Integer c = columns[first_col];
        Real v = values[first_col];
        columns[first_col] = columns[j];
        values[first_col] = values[j];
        columns[j] = c;
        values[j] = v;
      }
    }
  }
  if (arcaneIsCheck())
    checkValid();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixImpl::
assemble()
{
  IntegerConstArrayView rows_index = m_rows_index;
  IntegerArrayView columns = m_columns;
  RealConstArrayView values = m_values;
 
  SharedArray<Integer> new_rows_index(m_nb_row + 1);
  SharedArray<Integer> new_columns;
  SharedArray<Real> new_values;
 
  new_rows_index[0] = 0;
  for (Integer row = 0, nb_row = m_nb_row; row < nb_row; ++row) {
    Integer first_col = rows_index[row];
    for (Integer j = first_col, js = rows_index[row + 1]; j < js; ++j) {
      if (columns[j] >= 0) {
        new_columns.add(columns[j]);
        new_values.add(values[j]);
      }
    }
    new_rows_index[row + 1] = new_columns.size();
  }
 
  m_rows_index = new_rows_index;
  m_columns = new_columns;
  m_values = new_values;
  m_nb_element = new_values.size();
  if (arcaneIsCheck())
    checkValid();
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
void MatrixImpl::
dump(std::ostream& o)
{
  o << "(Matrix nb_row=" << m_nb_row << " nb_col=" << m_nb_column << ")\n";
  //Integer index = 0;
  for (Integer i = 0, is = m_nb_row; i < is; ++i) {
    o << "[" << i << "] ";
    Real sum = 0.0;
    for (Integer z = m_rows_index[i], zs = m_rows_index[i + 1]; z < zs; ++z) {
      Integer j = m_columns[z];
      Real v = m_values[z];
      sum += v;
      o << " [" << j << "]=" << v;
    }
    o << " S=" << sum << '\n';
  }
  //o << ")";
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix Matrix::
read(const String& filename)
{
  std::ifstream ifile(filename.localstr());
  Integer nb_x = 0;
  Integer nb_y = 0;
  ifile >> ws >> nb_x >> ws >> nb_y;
  //cout << "** ** N=" << nb_x << ' ' << nb_y << '\n';
  if (nb_x != nb_y)
    ARCANE_FATAL("not a square matrix nb_x={0} nb_y={1}", nb_x, nb_y);
  Matrix m(nb_x, nb_y);
  UniqueArray<Int32> rows_size(nb_x);
  UniqueArray<Int32> columns;
  RealUniqueArray values;
  for (Int32 x = 0; x < nb_x; ++x) {
    Int32 nb_column = 0;
    for (Int32 y = 0; y < nb_y; ++y) {
      Real v = 0.0;
      ifile >> v;
      if (!math::isZero(v)) {
        columns.add(y);
        values.add(v);
        ++nb_column;
      }
    }
    rows_size[x] = nb_column;
  }
  m.setRowsSize(rows_size);
  m.setValues(columns, values);
  return m;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
Matrix Matrix::
readHypre(const String& filename)
{
  const bool is_verbose = false;
  std::ifstream ifile(filename.localstr());
  if (!ifile.good())
    ARCANE_FATAL("Can not read Hypre matrix file='{0}'", filename);
  Int32 lower_i = 0;
  Int32 upper_i = 0;
  Int32 lower_j = 0;
  Int32 upper_j = 0;
  ifile >> lower_i >> upper_i >> lower_j >> upper_j;
  if (is_verbose) {
    std::cout << "** ** I=" << lower_i << ' ' << upper_i << '\n';
    std::cout << "** ** IJ" << lower_j << ' ' << upper_j << '\n';
  }
  if (lower_i != 0 || lower_j != 0)
    ARCANE_FATAL("Only lower_i==0 or lower_j==0 is supported (lower_i={0} lower_j={1})", lower_i, lower_j);
  Int32 nb_x = 1 + upper_i - lower_i;
  Int32 nb_y = 1 + upper_j - lower_j;
  if (is_verbose)
    std::cout << "** ** N=" << nb_x << ' ' << nb_y << '\n';
  if (nb_x != nb_y)
    ARCANE_FATAL("Not a square matrix nb_x={0} nb_y={1}", nb_x, nb_y);
 
  Matrix m(nb_x, nb_y);
  UniqueArray<Int32> rows_size(nb_x);
  UniqueArray<Int32> columns;
  UniqueArray<Real> values;
 
  Int32 nb_column = 0;
  Int32 last_row = 0;
  while (!ifile.eof()) {
    Int32 x = 0;
    Int32 y = 0;
    Real v = 0.0;
    ifile >> ws >> x >> ws >> y >> ws >> v >> ws;
    if (is_verbose)
      std::cout << "X=" << x << " Y=" << y << " V=" << v << "\n";
    if (!ifile.good())
      break;
    if (x != last_row) {
      // New row.
      if (is_verbose)
        std::cout << "New row last_row_size=" << nb_column << "\n";
      rows_size[last_row] = nb_column;
      nb_column = 0;
      last_row = x;
    }
    columns.add(y);
    values.add(v);
    ++nb_column;
  }
  rows_size[last_row] = nb_column;
  m.setRowsSize(rows_size);
  m.setValues(columns, values);
  return m;
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
} // namespace Arcane::MatVec
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/