Matrice_Morse.cpp

/****************************************************************************
* Copyright (c) 2026, CEA
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
* 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
* 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
* IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
* OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*
*****************************************************************************/
 
#include <Matrice_Morse.h>
#include <Sparskit.h>
#include <unordered_map>
#include <Matrice_Morse_Sym.h>
#include <Check_espace_virtuel.h>
#include <SFichier.h>
#include <Noms.h>
#include <ArrOfBit.h>
#include <Array_tools.h>
#include <TRUSTTrav.h>
#include <TRUSTTrav.h>
 
 
Implemente_instanciable_sans_constructeur(Matrice_Morse,"Matrice_Morse",Matrice_Base);
 
/*! @brief Writes the three arrays of the Morse storage structure to an output stream.
 *
 * @param (Sortie& s) an output stream
 * @return (Sortie& s) the modified output stream
 */
Sortie& Matrice_Morse::printOn(Sortie& s) const
{
  s << tab1_;
  s << tab2_;
  s << coeff_;
  s << m_ << finl;
  return s;
}
 
/*! @brief NOT IMPLEMENTED
 *
 * @param (Entree& s) an input stream
 * @return (Entree& s) the input stream
 * @throws NOT IMPLEMENTED
 */
Entree& Matrice_Morse::readOn(Entree& s)
{
  s >> tab1_;
  s >> tab2_;
  s >> coeff_;
  s >> m_;
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  return s;
}
 
Sortie& Matrice_Morse::imprimer(Sortie& s) const
{
  int n=nb_lignes();
  for(int i=0; i<n; i++)
    {
      s <<i << ": " <<finl;
      s << "--------------------------------" << finl;
      for (auto k=tab1_(i)-1; k<tab1_(i+1)-1; k++)
        {
          s << "("<<(tab2_(k)-1) << "),(" <<coeff_(k)<< ") "
            << " k= " << k << finl;
        }
      s <<finl;
    }
  return s;
}
 
Sortie& Matrice_Morse::imprimer_formatte(Sortie& s) const
{
  return imprimer_formatte(s,0);
}
 
Sortie& Matrice_Morse::imprimer_formatte(Sortie& s, int symetrie) const
{
  int numerotation_fortran=(tab1_[0]==1);
  for (int proc=0; proc<Process::nproc(); proc++)
    {
      if (proc==Process::me())
        {
          s << "Matrix morse on the processor " << proc << " : " << finl;
          int n=nb_lignes();
          Noms tab_imp;
          tab_imp.dimensionner(nb_colonnes());
          for(int i=0; i<n; i++)
            {
              for (int k=0; k<nb_colonnes(); k++)
                tab_imp[k]="  .  ";
              if (i<10)
                s <<i << " :" ;
              else
                s <<i << ":" ;
 
              if (symetrie)
                {
                  for (int j=0; j<i; j++)
                    {
                      for (auto k=tab1_(j)-numerotation_fortran; k<tab1_(j+1)-numerotation_fortran; k++)
                        if (tab2_(k)-numerotation_fortran==i)
                          tab_imp[j] = coeff_(k);
                    }
                  int ligne=tab2_(tab1_(i)-numerotation_fortran)-numerotation_fortran;
                  if (i!=ligne)
                    {
                      Cerr << "Problem detected on this Matrice_Morse_Sym." << finl;
                      Cerr << "The diagonal of the line " << ligne << " must be stored even if it is null." << finl;
                      exit();
                    }
                }
              for (auto k=tab1_(i)-numerotation_fortran; k<tab1_(i+1)-numerotation_fortran; k++)
                if (tab2_(k)+!numerotation_fortran==0)
                  Cerr<<"Line " <<i<< " no coefficient "<<k<<finl;
                else
                  {
                    if (coeff_(k)>=0)
                      tab_imp[tab2_(k)-numerotation_fortran]=" ";
                    else
                      tab_imp[tab2_(k)-numerotation_fortran]="";
                    tab_imp[tab2_(k)-numerotation_fortran] += (Nom)coeff_(k);
                  }
              for(int k=0; k<nb_colonnes(); k++)
                s<<tab_imp[k];
              s<<finl;
            }
        }
      Process::barrier();
    }
  return s;
}
 
Sortie& Matrice_Morse::imprimer_image(Sortie& s) const
{
  return imprimer_image(s,0);
}
 
Sortie& Matrice_Morse::imprimer_image(Sortie& s, int symetrie) const
{
  int numerotation_fortran=(tab1_[0]==1);
  for (int proc=0; proc<Process::nproc(); proc++)
    {
      if (proc==Process::me())
        {
          s << "Matrix morse on the processor " << proc << " : " << finl;
          int n=nb_lignes();
          Noms tab_imp;
          tab_imp.dimensionner(nb_colonnes());
          for(int i=0; i<n; i++)
            {
              for (int k=0; k<nb_colonnes(); k++)
                tab_imp[k]="\u2588\u2588";
              if (i<10)
                s <<i << " :" ;
              else
                s <<i << ":" ;
 
              if (symetrie)
                {
                  for (int j=0; j<i; j++)
                    {
                      for (auto k=tab1_(j)-numerotation_fortran; k<tab1_(j+1)-numerotation_fortran; k++)
                        if (tab2_(k)-numerotation_fortran==i)
                          tab_imp[j] = (std::abs(coeff_(k)) < 1e-20) ? "  " : "\u2592\u2592";
                    }
                  int ligne=tab2_(tab1_(i)-numerotation_fortran)-numerotation_fortran;
                  if (i!=ligne)
                    {
                      Cerr << "Problem detected on this Matrice_Morse_Sym." << finl;
                      Cerr << "The diagonal of the line " << ligne << " must be stored even if it is null." << finl;
                      exit();
                    }
                }
              for (auto k=tab1_(i)-numerotation_fortran; k<tab1_(i+1)-numerotation_fortran; k++)
                if (tab2_(k)+!numerotation_fortran==0)
                  Cerr<<"Line " <<i<< " no coefficient "<<k<<finl;
                else
                  tab_imp[tab2_(k)-numerotation_fortran] = (std::abs(coeff_(k)) < 1e-20) ? "  " : "\u2592\u2592";
 
              for(int k=0; k<nb_colonnes(); k++)
                s<<tab_imp[k];
              s<<finl;
            }
        }
      Process::barrier();
    }
  return s;
}
 
void Matrice_Morse::WriteFileMTX(const Nom& name) const
{
  if (Process::is_parallel())
    {
      Cerr << "Warning, matrix market format is not available yet in parallel." << finl;
      return;
    }
  Nom filename(Objet_U::nom_du_cas());
  filename += "_";
  filename += name;
  filename += ".mtx";
  SFichier mtx(filename);
  mtx.precision(14);
  mtx.setf(ios::scientific);
  int rows = nb_lignes();
  Cerr << "Matrix (" << rows << " lines) written into file: " << filename << " ... " << finl;
  mtx << "%%MatrixMarket matrix coordinate real " << (sub_type(Matrice_Morse_Sym, *this) ? "symmetric" : "general") << finl;
  Cerr << "Matrix (" << rows << " lines) written into file: " << filename << finl;
  mtx << "%%matrix" << finl;
  mtx << rows << " " << rows << " " << get_tab1()[rows] << finl;
  for (int row=0; row<rows; row++)
    for (auto j=get_tab1()[row]; j<get_tab1()[row+1]; j++)
      mtx << row+1 << " " << get_tab2()[j-1] << " " << get_coeff()[j-1] << finl;
}
 
/*! @brief Copy constructor of a Matrice_Morse.
 *
 * Copies each member of the given parameter.
 *
 * @param (Matrice_Morse& acopier) the Morse matrix to copy
 */
Matrice_Morse::Matrice_Morse(const Matrice_Morse& acopier) :Matrice_Base(),
  tab1_(acopier.tab1_),
  tab2_(acopier.tab2_),
  coeff_(acopier.coeff_),
  m_(acopier.m_),
  symetrique_(0),
  zero_(0)
{
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  is_stencil_up_to_date_ = acopier.is_stencil_up_to_date_ ;
}
/*! @brief Constructor of a square Morse matrix of order n that can store at most nnz non-zero elements.
 *
 *     Also the default constructor as both parameters have a default value.
 *
 * @param (int n) the order of the square matrix to construct
 * @param (int nnz) the number of non-zero elements the matrix can store.
 */
template<typename _SIZE_>
Matrice_Morse::Matrice_Morse(int n, _SIZE_ nnz) :
  morse_matrix_structure_has_changed_(1), symetrique_(0), zero_(0)
{
  dimensionner(n,nnz), sorted_ = 0;
  is_stencil_up_to_date_ = false ;
}
 
Matrice_Morse::Matrice_Morse()
{
  dimensionner(0,0);
  morse_matrix_structure_has_changed_=1;
  symetrique_ = 0;
  sorted_ = 0;
  zero_ = 0;
  is_stencil_up_to_date_ = false ;
}
 
 
/*! @brief Constructor of a Morse matrix with n rows and m columns that can store at most nnz non-zero elements.
 *
 * @param (int n) the number of rows of the matrix
 * @param (int m) the number of columns of the matrix
 * @param (int nnz) the number of non-zero elements the matrix can store.
 */
template<typename _SIZE_>
Matrice_Morse::Matrice_Morse(int n, int m, _SIZE_ nnz):
  morse_matrix_structure_has_changed_(1), symetrique_(0), zero_(0)
{
  dimensionner(n,m,nnz);
  is_stencil_up_to_date_ = false, sorted_ = 0 ;
}
 
 
Matrice_Morse::Matrice_Morse(int n, int nnz, const IntLists& voisins,
                             const DoubleLists& valeurs,
                             const DoubleVect& terme_diag)
  :  morse_matrix_structure_has_changed_(1), symetrique_(0) , zero_(0)
{
  dimensionner(n,n,nnz);
  remplir(voisins, valeurs, terme_diag);
  is_stencil_up_to_date_ = false, sorted_ = 0;
}
 
void Matrice_Morse::set_nb_columns( const int nb_col )
{
  m_ = nb_col;
}
 
void Matrice_Morse::set_symmetric( const int symmetric )
{
  symetrique_ = symmetric ;
}
 
/*! @brief Size the matrix with n lines and n columns and nnz zero-values coefficients
 *
 */
template<typename _SIZE_>
void Matrice_Morse::dimensionner(int n, _SIZE_ nnz)
{
  dimensionner(n,n,nnz);
  return ;
}
 
 
/*! @brief Resizes the sparse matrix, optionally adding new non-zero coefficients
 *
 *
 *    Parameter: const IntTab &Ind
 *       Meaning: array of size nc * 2
 *                where nc is the number of pairs (i,j)
 *                for the indices of the new coefficients
 *
 *
 */
void Matrice_Morse::dimensionner(const IntTab& Ind)
{
  if (Ind.size()==0) return; // Do nothing if the structure is empty
  int n_ancien = nb_lignes(), m_ancien = nb_colonnes();
 
  assert(Ind.nb_dim() == 2);
  assert(Ind.dimension(1) == 2);
 
  // Compute the new number of rows
  //   = max (old, row indices of new coefficients)
  //
  // and the new number of columns
  //   = max (old, column indices of new coefficients)
 
  int nInd = Ind.dimension(0);
  int n = 0;
  int m = 0;
  for (int i=0; i<nInd; i++)
    {
      if (n < Ind(i,0)) n = Ind(i,0);
      if (m < Ind(i,1)) m = Ind(i,1);
    }
  n++;
  m++;
  if (n < n_ancien) n = n_ancien;
  if (m < m_ancien) m = m_ancien;
 
  // Copies of the old index arrays
 
  auto tab1_temp(tab1_);
 
  // Initialize to the number of coefficients already present at each row
 
  tab1_.resize(n+1);
  m_ = m;
 
  for (int i=1; i<=n_ancien; i++)
    tab1_[i] = tab1_temp[i] - tab1_temp[i-1];
  for (int i=n_ancien+1; i<=n; i++)
    tab1_[i] = 0;
 
  // Iterate over the indices of new coefficients to check if they are
  // already present
 
  int i_nouveaux = 0;
  for (int i=0; i<nInd; i++)
    {
      int i0 = Ind(i,0);
      int i1 = Ind(i,1) + 1;
 
      int test_present = 0;
 
      if (i0 < n_ancien)
        {
          auto kmin = tab1_temp[i0]-1;
          auto kmax = tab1_temp[i0+1]-1;
          for (auto k=kmin; k<kmax; k++)
            if (tab2_[k] == i1)
              {
                test_present = 1;
                break;
              }
        }
      if (!test_present)
        {
          i_nouveaux++;
          tab1_[i0+1] += 1;
        }
    }
  if (i_nouveaux == 0)
    {
      tab1_=tab1_temp;
      return;
    }
 
  // New array of positions of the first coefficients of each row
  tab1_[0] = 1;
  for (int i=1; i<=n; i++)
    tab1_[i] += tab1_[i-1];
 
  auto nnz_ancien = tab2_.size_array();
  auto nnz = nnz_ancien + i_nouveaux;
 
  auto tab2_temp(tab2_);
  auto coeff_temp(coeff_);
  tab2_.resize(nnz);
  coeff_.resize(nnz);
 
  // Copy old coefficients and their column indices
  // into the new arrays
 
  tab2_ = -1;
  for (int i=0; i<n_ancien; i++)
    {
      for (auto j1 = tab1_temp[i]-1, j2 = tab1_[i]-1;
           j1 < tab1_temp[i+1]-1;
           j1++, j2++)
        {
          tab2_[j2] = tab2_temp[j1];
          coeff_[j2] = coeff_temp[j1];
        }
    }
 
  for (int i=0; i<nInd; i++)
    {
      int j0 = Ind(i,0);
      int j1 = Ind(i,1) + 1;
      auto k(tab1_(0));
      for (k=tab1_[j0]-1; tab2_[k] >= 0; k++)
        if (tab2_[k] == j1)
          {
            break;
          }
      if (tab2_[k] < 0)
        {
          tab2_[k] = j1;
          coeff_[k] = 0.0;
        }
    }
  // reorder coefficients... not optimal but for checking..
  coeff_=0;
  //
  {
    int nbis=nb_lignes();
    for(int i=0; i<nbis; i++)
      {
        for (auto k=tab1_(i)-1; k<tab1_(i+1)-1; k++)
          {
            for (auto k2=k; k2<tab1_(i+1)-1; k2++)
              {
                int j1=tab2_(k);
                int j2=tab2_(k2);
                if (j1>j2)
                  {
                    tab2_(k)=j2;
                    tab2_(k2)=j1;
                  }
              }
          }
      }
  }
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
}
 
/*! @brief Size the matrix with n lines, m columns with nnz zero-values coefficients
 *
 */
template<typename _SIZE_>
void Matrice_Morse::dimensionner(int n, int m, _SIZE_ nnz)
{
  tab2_.resize(nnz);
  coeff_.resize(nnz);
  m_=m;
 
  // check if tab1 has the right size and if tab1[n1]==nnz.
  if ( tab1_.size_array()!=(n+1) || (tab1_[n]-1)!=nnz )
    {
      tab1_.resize(n+1);
      tab1_=1;
    }
  tab1_.resize(n+1);
  tab1_[n]=nnz+1;
 
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
}
 
/*! @brief Initialize to the identity matrix (modif MT)
 *
 */
void Matrice_Morse::unite()
{
  coeff_ = 0.0;
  int i,n  = ordre();
  for (i=0; i<n; i++)
    operator()(i,i) = 1.0;
}
 
/*! @brief Returns the order of the matrix: - the number of rows if the matrix is square
 *
 *      - 0 otherwise
 *
 * @return (int) the order of the matrix
 */
int Matrice_Morse::ordre() const
{
  if(nb_lignes()==nb_colonnes())
    return nb_lignes();
  else
    return 0;
}
 
/*! @brief Method to check/clean the Matrice_Morse matrix: -Suppress coefficient defined several times
 *
 *  -elim_coeff_nul=0, zero coefficients are not removed from the matrix
 *  -elim_coeff_nul=1, zero coefficients are removed from the matrix
 *  -elim_coeff_nul=2, zero and near-zero coefficients are removed from the matrix
 *
 */
void Matrice_Morse::compacte(int elim_coeff_nul)
{
  int n=nb_lignes();
  int coeff_nuls=0;
  int coeff_quasi_nuls=0;
  auto tab_elim_coeff(tab2_); // Possibly BigArrOfInt
  tab_elim_coeff = 0;
  if (elim_coeff_nul)
    {
      ArrOfDouble tab_coeff_max(n);
      tab_coeff_max = 0.;
      // Search for off-diagonal zero coefficients to remove from the Morse matrix
      {
        ArrOfInt tab_cnt(1);
        tab_cnt = 0;
        auto tab1 = tab1_.view_ro();
        CDoubleArrView coeff = coeff_.view_ro();
        DoubleArrView coeff_max = tab_coeff_max.view_rw();
        auto elim_coeff = tab_elim_coeff.view_rw();
        IntArrView cnt = tab_cnt.view_rw();
        Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n), KOKKOS_LAMBDA(const int i)
        {
          auto k1 = tab1(i)-1;
          auto k2 = tab1(i+1)-1;
          for (auto k = k1; k < k2; k++)
            {
              double abs_c = Kokkos::fabs(coeff(k));
              if (abs_c > coeff_max(i)) coeff_max(i) = abs_c;
              if (coeff(k) == 0)
                {
                  Kokkos::atomic_add(&cnt(0), 1);
                  elim_coeff(k) = 1;
                }
            }
        });
        end_gpu_timer(__KERNEL_NAME__);
        coeff_nuls = tab_cnt(0);
      }
 
      if (elim_coeff_nul==2)
        {
          // Search for near-zero off-diagonal coefficients (1.e-12 smaller than the largest coefficient in the row) to remove from the Morse matrix
          const double eps = Objet_U::precision_geom;
          ArrOfInt tab_cnt(1);
          tab_cnt = 0;
          auto tab1 = tab1_.view_ro();
          CDoubleArrView coeff = coeff_.view_ro();
          CDoubleArrView coeff_max = tab_coeff_max.view_ro();
          IntArrView elim_coeff = tab_elim_coeff.view_rw();
          IntArrView cnt = tab_cnt.view_rw();
          Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n), KOKKOS_LAMBDA(const int i)
          {
            double cm = coeff_max(i);
            if (!est_egal(cm, 0., eps) && cm < 1e10)
              {
                auto k1 = tab1(i) - 1;
                auto k2 = tab1(i + 1) - 1;
                for (auto k = k1; k < k2; k++)
                  if (coeff(k) != 0 && est_egal(Kokkos::fabs(coeff(k)) / cm, 0., eps))
                    {
                      Kokkos::atomic_add(&cnt(0), 1);
                      elim_coeff(k) = 1;
                    }
              }
          });
          end_gpu_timer(__KERNEL_NAME__);
          coeff_quasi_nuls = tab_cnt(0);
        }
    }
  // Recherche des coefficients doublons
  int nb_doublons=0;
  {
    auto tab1 = tab1_.view_ro();
    CIntArrView tab2 = tab2_.view_ro();
    CDoubleArrView coeff = coeff_.view_ro();
    IntArrView elim_coeff = tab_elim_coeff.view_rw();
    ArrOfInt tab_doublons(1);
    tab_doublons = 0;
    ArrOfInt tab_error(1);
    tab_error = 0;
    IntArrView doublons = tab_doublons.view_rw();
    IntArrView error = tab_error.view_rw();
    Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n), KOKKOS_LAMBDA(const int i)
    {
      auto k1 = tab1(i)-1;
      auto k2 = tab1(i+1)-1;
      int jmax = -1; // Highest column of a coefficient in the line i
      for (auto k = k1; k < k2; k++)
        {
          int j = tab2(k)-1;
          if (j > jmax)
            jmax = j;
          else
            {
              // Found a column j lower than jmax, check if not defined before:
              for (auto kk = k-1; kk >= k1; kk--)
                {
                  int jj = tab2(kk)-1;
                  if (jj == j)
                    {
                      // Already defined!
                      Kokkos::atomic_add(&doublons(0), 1);
                      elim_coeff(k) = 1;
                      // Check if same coefficients:
                      if (coeff(kk) != coeff(k))
                        Kokkos::atomic_add(&error(0), 1);
                      break;
                    }
                }
            }
        }
    });
    end_gpu_timer(__KERNEL_NAME__);
    nb_doublons = tab_doublons(0);
    if (tab_error(0))
      {
        Cerr << "Error in a Matrix Morse: duplicate entries with different values!" << finl;
        exit();
      }
  }
 
  auto nnz(tab1_(0));
  nnz=0;
  if (nb_doublons || coeff_nuls || coeff_quasi_nuls)
    {
      // Step 1: Count kept entries per row (parallel_for over rows)
      ArrOfInt tab_kept_per_row(n);
      {
        auto tab1 = tab1_.view_ro();
        CIntArrView elim_coeff = tab_elim_coeff.view_ro();
        IntArrView kept_per_row = tab_kept_per_row.view_wo();
        Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n), KOKKOS_LAMBDA(const int i)
        {
          int count = 0;
          auto k1 = tab1(i)-1;
          auto k2 = tab1(i+1)-1;
          for (auto k = k1; k < k2; k++)
            if (!elim_coeff(k)) count++;
          kept_per_row(i) = count;
        });
        end_gpu_timer(__KERNEL_NAME__);
      }
 
      // Step 2: Save old tab1_ (needed for source offsets in scatter step)
      auto old_tab1(tab1_);
 
      // Step 3: Update tab1_ via prefix scan (updates tab1_(1..n), tab1_(0)=1 unchanged)
      using tab1_scan_t = decltype(nnz);
      {
        auto tab1 = tab1_.view_rw();
        CIntArrView kept_per_row = tab_kept_per_row.view_ro();
        Kokkos::parallel_scan(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n), KOKKOS_LAMBDA(const int i, tab1_scan_t& update, const bool final)
        {
          update += kept_per_row(i);
          if (final) tab1(i+1) = update + 1;
        });
        end_gpu_timer(__KERNEL_NAME__);
      }
 
      // Step 4: Out-of-place scatter of coeff_ and tab2_ to new positions (parallel_for over rows)
      // Safe because new_pos(i) <= old_pos(i) always, and rows are processed independently
      nnz = tab1_[n] - 1;
      auto new_coeff(coeff_);
      auto new_tab2(tab2_);
      {
        auto tab1 = tab1_.view_ro();
        auto old_tab1_ro = old_tab1.view_ro();
        CDoubleArrView coeff_src = coeff_.view_ro();
        CIntArrView tab2_src = tab2_.view_ro();
        DoubleArrView coeff_dst = new_coeff.view_wo();
        IntArrView tab2_dst = new_tab2.view_wo();
        CIntArrView elim_coeff = tab_elim_coeff.view_ro();
        Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, n), KOKKOS_LAMBDA(const int i)
        {
          auto new_pos = tab1(i) - 1;
          auto k1 = old_tab1_ro(i)-1;
          auto k2 = old_tab1_ro(i+1)-1;
          for (auto k = k1; k < k2; k++)
            if (!elim_coeff(k))
              {
                coeff_dst(new_pos) = coeff_src(k);
                tab2_dst(new_pos) = tab2_src(k);
                new_pos++;
              }
        });
        end_gpu_timer(__KERNEL_NAME__);
      }
 
      // Step 5: Copy compacted data back
      {
        auto tab2 = tab2_.view_rw();
        auto coeff = coeff_.view_rw();
        CIntArrView new_tab2_ro = new_tab2.view_ro();
        CDoubleArrView new_coeff_ro = new_coeff.view_ro();
        Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__), range_1D(0, nnz), KOKKOS_LAMBDA(const int i)
        {
          tab2(i) = new_tab2_ro(i);
          coeff(i) = new_coeff_ro(i);
        });
        end_gpu_timer(__KERNEL_NAME__);
      }
    }
  else
    {
      nnz = tab1_[n] - 1;
    }
 
  // Resize the arrays
  tab2_.resize(nnz);
  coeff_.resize(nnz);
 
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  assert_check_morse_matrix_structure( );
}
 
/*! @brief Assignment operator from one Matrice_Morse to another.
 *
 * @param (Matrice_Morse& a) the right-hand side of the assignment
 */
Matrice_Morse& Matrice_Morse::operator=(const Matrice_Morse& a )
{
  tab1_ = a.get_tab1();
  tab2_ = a.get_tab2();
  coeff_ = a.get_coeff();
  m_=a.nb_colonnes();
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  is_stencil_up_to_date_=a.is_stencil_up_to_date();
  return(*this);
}
 
/*! @brief *this = transpose of a.
 *
 * @param (Matrice_Morse& a) the matrix to transpose
 */
Matrice_Morse& Matrice_Morse::transpose(const Matrice_Morse& a)
{
  int n=a.nb_lignes();
  int jk=nb_lignes();
  int job=1;
  int ipos=1;
  int m=a.nb_colonnes();
  int l=nb_lignes();
  if(m!=jk)
    {
      Cerr << "Matrice_Morse::transpose bad dimensions" << finl;
      exit();
    }
  m=a.nb_lignes();
  l=nb_colonnes();
  if(m!=l)
    {
      Cerr << "Matrice_Morse::transpose bad dimensions" << finl;
      exit();
    }
 
  for(int i=0; i<=jk; i++ ) tab1_[i] = 0 ;
  for(int i=0; i<n; i++)
    {
      for(auto k=a.tab1_[i]-1; k<a.tab1_[i+1]-1; k++)
        {
          int j = a.tab2_[k] ;
          tab1_[j] = tab1_[j] +1 ;
        }
    }
  tab1_[0] = ipos ;
  for(int i=0; i<jk; i++) tab1_[i+1] = tab1_[i] + tab1_[i+1] ;
  for(int i=0; i<n; i++)
    {
      for(auto k=a.tab1_[i]-1; k<a.tab1_[i+1]-1; k++)
        {
          int j = a.tab2_[k]-1 ;
          auto next = tab1_[j] ;
          if (job == 1) coeff_[next-1] = a.coeff_[k] ;
          tab2_[next-1] = i+1 ;
          tab1_[j]    = next+1 ;
        }
    }
  for(int i=jk-1; i>=0; i--) tab1_[i+1] = tab1_[i] ;
  tab1_[0] = ipos ;
 
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  return(*this);
}
 
 
//A=x*A with x a diagonal matrix stored as a vector
//the same method can be used to store the result in
//a different matrix than the initial one
Matrice_Morse& Matrice_Morse::diagmulmat(const DoubleVect& x)
{
  int m=nb_lignes();
  int l=0;
  int n=x.size_array();
  if(n!=m)
    {
      Cerr << "Matrice_Morse::diagmulmat bad dimensions" << finl;
      exit();
    }
  F77NAME(DIAMUA)(&m ,&l,
                  coeff_.addr(),tab2_.addr(),reinterpret_cast<const int*>(tab1_.addr()),x.addr(),
                  coeff_.addr(),tab2_.addr(),reinterpret_cast<int*>(tab1_.addr()));
  return(*this);
}
 
//extract the upper part of a Morse matrix
//the result matrix is the calling one
Matrice_Morse& Matrice_Morse::partie_sup(const Matrice_Morse& a)
{
  int m=nb_lignes();
  int n=a.nb_lignes();
  if(m!=n)
    {
      Cerr << "Matrice_Morse::partie_sup : bad dimensions m!=n." << finl;
      exit();
    }
  double t;
  auto ko(tab1_(0));
  auto kfirst (ko);
  auto kdiag(ko);
  ko = -1;
  for(int i=0; i< n; i++)
    {
      kfirst = ko + 1 ;
      kdiag = -1 ;
      for(auto k = a.tab1_[i]-1; k< a.tab1_[i+1]-1; k++)
        {
          if (a.tab2_[k]-1 >= i)
            {
              ko++ ;
              coeff_[ko] = a.coeff_[k] ;
              tab2_[ko] = a.tab2_[k] ;
              if (a.tab2_[k] == i) kdiag = ko ;
            }
        }
      if (kdiag != -1 && kdiag != kfirst)
        {
          t = coeff_[kdiag] ;
          coeff_[kdiag] = coeff_[kfirst] ;
          coeff_[kfirst] = t ;
          { int ktmp = tab2_[kdiag] ; tab2_[kdiag] = tab2_[kfirst] ; tab2_[kfirst] = ktmp ; }
        }
      tab1_[i] = kfirst+1 ;
    }
  auto nnz = (ko + 1) ;
 
  tab1_[n] = (nnz) + 1 ;
  tab2_.resize( nnz );
  coeff_.resize( nnz );
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  return(*this);
}
 
/*! @brief Operation de multiplication-accumulation (saxpy) matrice vecteur.
 *
 * Operation: resu = resu + A*x
 *
 */
DoubleVect& Matrice_Morse::ajouter_multvect_(const DoubleVect& tab_x,DoubleVect& tab_resu) const
{
  assert_check_morse_matrix_structure();
  const int n = tab1_.size_array() - 1;
  assert(tab_x.size_array() == nb_colonnes());
  // Test in this order as the size() attribute may be invalid:
  assert(tab_resu.size_array() == n || tab_resu.size() == n);
  // If matrix, x, resu are on device, we compute on the device to avoid expensive copy during TRUST GCP:
  if (tab_x.isDataOnDevice() && tab_resu.isDataOnDevice() && coeff_.isDataOnDevice())
    {
      //if (tab_x.line_size()>1) Process::exit("line_size>1 pour x dans Matrice_Morse::ajouter_multvect_");
      // Faster implementation on GPU (ToDo Kokkos: future, use Kokkos kernel?)
      auto tab1 = tab1_.view_ro();
      CIntArrView tab2 = tab2_.view_ro();
      CDoubleArrView coeff = coeff_.view_ro();
      CDoubleArrView x = tab_x.view_ro();
      DoubleArrView resu = tab_resu.view_rw();
      Kokkos::parallel_for(start_gpu_timer(__KERNEL_NAME__),
                           Kokkos::RangePolicy<>(0, n), KOKKOS_LAMBDA(
                             const int i)
      {
        auto start = tab1(i)-1;
        auto end = tab1(i + 1)-1;
        double tmp {};
 
        for (auto k = start; k < end; k++)
          {
            int j = tab2(k) - 1;
            tmp+= coeff(k) * x(j);
          }
        resu(i) += tmp;
      });
      end_gpu_timer(__KERNEL_NAME__);
    }
  else
    {
      tab_x.ensureDataOnHost();
      tab_resu.ensureDataOnHost();
      coeff_.ensureDataOnHost();
      // Fast CPU (old) implementation with pointer:
      const DoubleVect& x = tab_x;
      DoubleVect& resu = tab_resu;
      const auto *tab1_ptr = tab1_.addr() + 1;
      const int *tab2_ptr = tab2_.addr();
      const double *coeff_ptr = coeff_.addr();
      const double *x_fortran = x.addr() - 1; // to index x with a Fortran index
      auto k_fortran = 1; // Fortran index in tab2 and coeff
      for (int i = 0; i < n; i++, tab1_ptr++)
        {
          const auto kmax = *tab1_ptr; // tab1_[i+1] = Fortran index in tab2_
          assert(kmax >= k_fortran && kmax <= tab2_.size_array() + 1);
          double t = resu[i];
          assert(k_fortran == tab1_[i] && tab2_ptr == tab2_.addr() + (k_fortran - 1));
          for (; k_fortran < kmax; k_fortran++, tab2_ptr++, coeff_ptr++)
            {
              int colonne = *tab2_ptr; // Fortran index
              assert(colonne >= 1 && colonne <= nb_colonnes());
              t += (*coeff_ptr) * x_fortran[colonne];
            }
          resu[i] = t;
        }
    }
  return tab_resu;
}
 
// Multiplication of the matrix by a vector x, using only the non-shared real items for x
ArrOfDouble& Matrice_Morse::ajouter_multvect_(const ArrOfDouble& x,ArrOfDouble& resu,ArrOfInt& est_reel_pas_com) const
{
  ToDo_Kokkos("critical ?");
  assert_check_morse_matrix_structure( );
  int n = nb_lignes();
 
  assert(nb_colonnes()==x.size_array());
  assert(n==resu.size_array());
  for(int i=0; i<n; i++)
    {
      double t = 0.0;
      for (auto k=tab1_(i)-1; k<tab1_(i+1)-1; k++)
        {
          int j=tab2_(k)-1;
          if (est_reel_pas_com[j]) t += coeff_(k)*x[j];
        }
      resu[i] += t ;
    }
  return resu;
}
 
/*! @brief Matrix-matrix multiply-accumulate operation (saxpy) (matrix X represented by an array)
 *
 *     Operation: RESU = RESU + A*X
 *
 * @param (DoubleTab& x) the matrix to multiply
 * @param (DoubleTab& resu) the result matrix of the operation
 * @return (DoubleTab&) the result matrix of the operation
 */
DoubleTab& Matrice_Morse::ajouter_multTab_(const DoubleTab& x,DoubleTab& resu) const
{
 
  if ( (x.nb_dim() == 1) && (resu.nb_dim() == 1))
    {
      ajouter_multvect(x,resu);
      return resu;
    }
 
  assert_check_morse_matrix_structure( );
  int nb_comp = x.dimension(1);
 
  assert(resu.dimension(1) == nb_comp);
  double* t= new double[nb_comp];
  int ncomp;
  int n=nb_lignes();
  for(int i=0; i<n; i++)
    {
      for (ncomp=0; ncomp<nb_comp; ncomp++)
        t[ncomp] = 0.0;
      for (auto k=tab1_(i)-1; k<tab1_(i+1)-1; k++)
        for (ncomp=0; ncomp<nb_comp; ncomp++)
          t[ncomp] += coeff_(k)*x(tab2_(k)-1,ncomp);
      for (ncomp=0; ncomp<nb_comp; ncomp++)
        resu(i,ncomp) += t[ncomp] ;
    }
  delete [] t;
  return resu;
}
 
 
/*! @brief Matrix-vector multiply-accumulate operation (saxpy), by the transposed matrix.
 *
 *     Operation: resu = resu + A^{T}*x
 *
 * @param (DoubleVect& x) the vector to multiply
 * @param (DoubleVect& resu) the result vector of the operation
 * @return (DoubleVect&) the result vector of the operation
 */
DoubleVect& Matrice_Morse::ajouter_multvectT_(const DoubleVect& x,DoubleVect& resu) const
{
  assert_check_morse_matrix_structure( );
 
  int n=nb_lignes();
  for(int i=0; i<n; i++)
    {
      double xi = x(i);
      for (auto k=tab1_(i)-1; k<tab1_(i+1)-1; k++)
        resu(tab2_(k)-1) += coeff_(k) * xi;
    }
  return resu;
}
 
// Multiplication of the transposed matrix by a vector x, using only the non-shared real items
ArrOfDouble& Matrice_Morse::ajouter_multvectT_(const ArrOfDouble& x,ArrOfDouble& resu,ArrOfInt& est_reel_pas_com) const
{
  assert_check_morse_matrix_structure( );
  int n=nb_lignes();
 
  assert(n==x.size_array());
  assert(nb_colonnes()==resu.size_array());
  for(int i=0; i<n; i++)
    {
      if (est_reel_pas_com[i])
        {
          double xi = x[i];
          for (auto k=tab1_(i)-1; k<tab1_(i+1)-1; k++)
            resu[tab2_(k)-1] += coeff_(k) * xi;
        }
    }
  return resu;
}
 
/*! @brief Friend function (outside the class) of the Matrice_Morse class. Addition of 2 Morse-format matrices.
 *
 *     Operation: returns (A+B)
 *
 * @param (Matrice_Morse& A) a Morse-format matrix
 * @param (Matrice_Morse& B) a Morse-format matrix
 * @return (Matrice_Morse) the result of the operation
 */
Matrice_Morse operator+(const Matrice_Morse& A , const Matrice_Morse& B )
{
  int nrow=A.nb_lignes();
  int ncol=A.nb_colonnes();
  Matrice_Morse C;
  // PL: before sizing to nzmax we check if A and B don't happen to have the same structure...
  // This avoids a memory spike caused by matrix addition in Equation_base::dimensionner_matrice
  auto nzmax = A.has_same_morse_matrix_structure(B) ? A.nb_coeff() : A.nb_coeff() + B.nb_coeff();
  C.dimensionner(nrow, ncol, nzmax);
#ifndef TRUST_USE_GPU
  // Fortran call cause faster on serail version on some Baltik:
  int job = 1;
  int ierr = -1;
  IntVect iw(ncol);
  F77NAME(APLB)(&nrow, &ncol, &job, A.get_coeff().addr(), A.get_tab2().addr(), A.get_tab1().addr(),
                B.get_coeff().addr(), B.get_tab2().addr(), B.get_tab1().addr(), C.get_set_coeff().addr(),
                C.get_set_tab2().addr(), C.get_set_tab1().addr(),
                &nzmax, iw.addr(), &ierr);
#else
  // Algorithm (per row i):
  //   1. Collect entries from row i of A and B into a small temporary buffer
  //   2. Sort by column index
  //   3. Merge duplicate columns (accumulate coefficients)
  //   4. Write result into C and advance c_tab1
  //
  // Time: O((nnz_A + nnz_B) * log(max_nnz_per_row))  [sort dominates]
  // Space: O(max_nnz_per_row_A + max_nnz_per_row_B)   [reused buffer]
  //
  // ToDo: Kokkos parallel version for GPU once CPU version is validated
  const auto& a_tab1 = A.get_tab1();
  const auto& a_tab2 = A.get_tab2();
  const auto& a_coeff = A.get_coeff();
  const auto& b_tab1 = B.get_tab1();
  const auto& b_tab2 = B.get_tab2();
  const auto& b_coeff = B.get_coeff();
  auto& c_tab1 = C.get_set_tab1();
  auto& c_tab2 = C.get_set_tab2();
  auto& c_coeff = C.get_set_coeff();
 
  using idx_t = std::remove_reference_t<decltype(c_tab1[0])>;
  idx_t nnz_c = 0; // running count of non-zeros written into C (0-based offset)
  c_tab1[0] = 1;    // 1-based (Morse/Fortran convention)
 
  std::unordered_map<int, idx_t> col_to_pos;
  col_to_pos.reserve(256);
  for (int i = 0; i < nrow; ++i)
    {
      col_to_pos.clear();
 
      // Step 1: copy A row i into C, recording each column's position
      for (auto k = a_tab1[i] - 1; k < a_tab1[i + 1] - 1; ++k)
        {
          c_tab2[nnz_c] = (int) a_tab2[k];
          c_coeff[nnz_c] = a_coeff[k];
          col_to_pos[(int) a_tab2[k]] = nnz_c;
          ++nnz_c;
        }
 
      // Step 2: merge B row i — accumulate if column already in C, else append
      for (auto k = b_tab1[i] - 1; k < b_tab1[i + 1] - 1; ++k)
        {
          const int col = (int) b_tab2[k];
          auto it = col_to_pos.find(col);
          if (it != col_to_pos.end())
            c_coeff[it->second] += b_coeff[k]; // column shared with A: accumulate
          else
            {
              c_tab2[nnz_c] = col;
              c_coeff[nnz_c] = b_coeff[k];
              col_to_pos[col] = nnz_c;
              ++nnz_c;
            }
        }
 
      c_tab1[i + 1] = nnz_c + 1; // 1-based pointer to start of next row
    }
#endif
  const auto nnz = C.tab1_[nrow] - 1;
  C.get_set_tab2().resize(nnz);
  C.get_set_coeff().resize(nnz);
  C.morse_matrix_structure_has_changed_ = 1, C.sorted_ = 0;
  return(C);
}
 
bool Matrice_Morse::has_same_morse_matrix_structure(const Matrice_Morse& A) const
{
  int nrow = A.nb_lignes();
  for (int i = 0; i < nrow; i++)
    if (tab1_(i) != A.tab1_(i))
      return false;
  auto ncoeff = tab2_.size_array(), ncoeff_A = A.tab2_.size_array();
  if (ncoeff != ncoeff_A) return false;
 
  for (auto i = 0; i < ncoeff; i++)
    if (tab2_(i) != A.tab2_(i))
      return false;
  return true;
}
 
/*! @brief Computes the solution of the linear system: A * solution = secmem.
 *
 * The method used is GMRES preconditioned with ILUT.
 *   WARNING: this method has probably never been tested in parallel
 *
 * @param (DoubleVect& secmem) the right-hand side of the linear system
 * @param (DoubleVect& solution) the solution of the system
 * @param (double coeff_seuil)
 * @return (int) always returns 1
 * @throws Error in ilut 'matrix may be wrong' dixit SAAD
 * @throws Error in ilut: overflow in L
 * @throws Error in ilut: overflow in U
 * @throws Illegal value for lfil: probably a memory corruption
 * @throws Empty row encountered
 * @throws Null pivot encountered! at step
 * @throws Something abnormal has happened: it is preferable to stop.
 */
// Delegates to the 4-arg version with max_iter=-1 (retry-on-failure mode, maxits=ordre())
int Matrice_Morse::inverse(const DoubleVect& secmem, DoubleVect& solution,
                           double coeff_seuil) const
{
  return inverse(secmem, solution, coeff_seuil, -1);
}
 
// Solves A*solution=secmem using ILUT-preconditioned PGMRES.
// max_iter<0: use ordre() as iteration limit and retry with stronger preconditioner on failure (returns 1).
// max_iter>=0: use max_iter as iteration limit and return 0 on failure (used by hyperbolic implicit).
int Matrice_Morse::inverse(const DoubleVect& secmem, DoubleVect& solution,
                           double coeff_seuil, int max_iter) const
{
  if (Process::is_parallel())
    {
      Cerr << "Matrice_Morse::inverse has never been tested in parallel" << finl;
      Cerr << "Try 'Solveur Gmres { diag }' or 'Solveur Petsc Gmres { precond diag { } }'" << finl;
      Cerr << "instead of 'Solveur Gmres { }' which is not parallelized yet." << finl;
      exit();
    }
 
  const bool retry_on_failure = (max_iter < 0);
 
  DoubleVect toto(secmem);
  int prems=1;                         // recompute L and U only when prems=1
  int lf_min = 10;
  int lf = std::min(lf_min, ordre()/2); // fill level for ILUT
  int nn = ordre();
  int ima = std::min(lf_min, ordre()/2); // Krylov space dimension
  IntVect ju, jlu;
  DoubleVect alu, vv;
  DoubleVect Resini(toto);
 
  int ie=1;
  auto n2 = nb_coeff()+(2*lf*nn); // number of non-zeros in LU
 
  double r, coeff_seuilr;
 
precond:
  if (prems)
    {
      int iw = (int)(n2 + 2);
      ju.resize(nn);
      jlu.resize(iw);
      alu.resize(iw);
      double to = 1.e-10; // drop tolerance for ILUT
      DoubleVect w(nn+1);
      IntVect jw(2*nn);
      set_tab1_int32();
      F77NAME(ILUT)(&nn, coeff_.addr(), tab2_.addr(), get_tab1_int32().addr(), &lf,
                    &to, alu.addr(), jlu.addr(), ju.addr(),
                    &iw, w.addr(), jw.addr(), &ie);
      switch(ie)
        {
        case  0:
          break;
        case -1:
          Cerr << "Error in ilut 'matrix may be wrong' dixit SAAD" << finl;
          exit();
          break;
        case -2:
          Cerr << "Error in ilut : overflow in L" << finl;
          exit();
          break;
        case -3:
          Cerr << "Error in ilut : overflow in U" << finl;
          exit();
          break;
        case -4:
          Cerr << "Illegal value for lfil : it may be a memory trouble" << finl;
          exit();
          break;
        case -5:
          Cerr << "Empty line met" << finl;
          exit();
          break;
        default:
          Cerr << "Pivot null met ! at step " << ie << finl;
          exit();
          break;
        }
      prems=0;
    }
 
  vv.resize(nn*(ima+1));
  assert_espace_virtuel_vect(solution);
  multvect(solution, Resini);
  Resini -= toto;
  r = mp_prodscal(Resini, Resini);
  r = sqrt(r);
  Cout << " Initial residu : " << r << finl;
  coeff_seuilr = (r == 0) ? DMAXFLOAT : coeff_seuil/r;
  Resini = toto;
  int minits = 10;
  int maxits = std::max(minits, retry_on_failure ? nn : max_iter);
  int io = 0;
  F77NAME(PGMRES)(&nn, &ima, toto.addr(), solution.addr(), vv.addr(), &coeff_seuilr,
                  &maxits, &io, coeff_.addr(), tab2_.addr(), get_tab1_int32().addr(),
                  alu.addr(), jlu.addr(), ju.addr(), &ie);
  switch(ie)
    {
    case 0:
      Cout << "     ** PGMRES has converged **" << finl;
      break;
    case 1:
      Cout << "     ** No convergence after " << maxits << " iterations **" << finl;
      if (retry_on_failure)
        {
          toto = Resini;
          if (lf < 50)
            {
              lf += 5;
              Cerr << "  The degree of the preconditioning matrix LU is increased: " << lf << finl;
              n2 = (int)tab2_.size_array()+(2*lf*nn);
              prems = 1;
              goto precond;
            }
        }
      else
        return 0;
      break;
    case -1:
      Cerr << "Convergence after 0 iterations !! 'stationnary state may be obtained'" << finl;
      break;
    default:
      Cerr << "Something abnormal has happened : it is preferable to stop." << finl;
      exit();
    }
  return 1;
}
 
 
/*! @brief Matrix-vector multiplication operator: scales rows of the matrix by the corresponding
 *
 *     coefficients of the vector passed as a parameter.
 *     A *= x performs the following scalings:
 *       A(i,:) = A(i,:) * x(i), for all rows i of A
 *
 * @param (DoubleVect& x) scaling vector
 * @return (Matrice_Morse&) the result of the operation (*this)
 */
Matrice_Morse& Matrice_Morse::operator *=(const DoubleVect& x)
{
  for(int i = 0; i<nb_lignes(); i++)
    for(auto k = tab1_(i)-1; k<tab1_(i+1)-1; k++)
      coeff_(k) *= x(i);
  return *this;
}
 
 
/*! @brief Assigns the product of 2 Morse matrices A and B to this object.
 *
 * Operation: this = A * B
 *
 * @param (Matrice_Morse& A) a Morse-format matrix
 * @param (Matrice_Morse& B) a Morse-format matrix
 * @return (Matrice_Morse&) the result of the operation (*this)
 */
Matrice_Morse& Matrice_Morse::affecte_prod(const Matrice_Morse& a, const Matrice_Morse& b)
{
  int nrow= a.nb_lignes();                // number of rows of A
  int ncol= b.nb_colonnes();                // number of columns of B
  //assert(nrow==ncol);
  // Jloi non?
  assert(a.nb_colonnes()==b.nb_lignes());
  tab1_.resize(nrow+1);
  m_ = ncol;
  int job = 1 ;                      // retrieve everything (tab1, tab2, coeff of matrice_resu)
  auto  nzmax = nb_coeff();                // max number of values of the result matrix
  if(nzmax==0)
    {
      nzmax=a.nb_coeff();
      tab2_.resize(nzmax);
      coeff_.resize(nzmax);
      assert(nzmax==nb_coeff());
    }
  IntVect iw(ncol+1);                        // work array
  double scal=0. ;
  int ii, jj ;
  int values = 0;
  if (job != 0) values = 1 ;
  auto len = -1 ;
  tab1_[0] = 1  ;
  iw = -1 ;
  for(ii=0; ii< nrow; ii++)
    {
      for(auto ka=a.tab1_[ii]-1; ka < a.tab1_[ii+1]-1; ka++)
        {
          if (values == 1) scal = a.coeff_[ka] ;
          jj   = a.tab2_[ka] - 1 ;
          for (auto kb=b.tab1_[jj]-1; kb < b.tab1_[jj+1]-1; kb++)
            {
              int jcol = b.tab2_[kb] -1 ;
              int jpos = iw[jcol]  ;
              if (jpos == -1)
                {
                  len++ ;
                  if (len > nzmax-1)
                    {
                      // Cerr << "Matrice_Morse::affect_prod len > nzmax -1 " << nzmax << finl;
                      nzmax *= 2;
                      coeff_.resize(nzmax);
                      tab2_.resize(nzmax);
                    }
                  tab2_[len] = jcol + 1 ;
                  iw[jcol]= (int)len ;
                  if (values == 1) coeff_[len]  = scal*b.coeff_[kb] ;
                }
              else
                {
                  if (values == 1) coeff_[jpos] += scal*b.coeff_[kb] ;
                }
            }
        }
 
      for (auto k=tab1_[ii]-1; k < len+1 ; k++) iw[tab2_[k]-1] = -1 ;
      tab1_[ii+1] = (len+1) + 1 ;
    }
 
  coeff_.resize(tab1_[nrow]-1);
  tab2_.resize(tab1_[nrow]-1);
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  return *this;
}
 
 
 
/*! @brief Friend function (outside the class) of the Matrice_Morse class. Scaling of the matrix by a scalar: multiplies all
 *
 *     elements of the matrix by a scalar.
 *     Operation: returns x*A
 *
 * @param (double x) scaling value
 * @param (Matrice_Morse& B) a Morse-format matrix
 * @return (Matrice_Morse) the result of the operation
 */
Matrice_Morse operator *(double x , const Matrice_Morse& A)
{
  Matrice_Morse mat_res(A);
  mat_res.coeff_*=x;
  return(mat_res);
}
 
/*! @brief Unary negation operator, returns the opposite of the matrix: - A. Calls operator*(double,const Matrice_Morse&)
 *
 * @return (Matrice_Morse) the result of the call to operator*(double,const Matrice_Morse&)
 */
Matrice_Morse Matrice_Morse::operator -() const
{
  return((-1)*(*this));
}
 
 
/*! @brief DOES NOTHING
 *
 * @param (Matrice_Morse&) a Morse matrix
 * @return (Matrice_Morse&) always returns *this
 */
Matrice_Morse& Matrice_Morse::operator +=(const Matrice_Morse& A)
{
  // PL: Before doing expensive RAM operations, we check
  // if it is not the same structure:
  if (has_same_morse_matrix_structure(A))
    {
      auto size = A.nb_coeff();
      const auto& coeff = A.get_coeff();
      for (auto i=0; i<size; i++)
        coeff_(i)+=coeff(i);
    }
  else
    {
      *this = *this + A;
      morse_matrix_structure_has_changed_ = 1, sorted_ = 0;
    }
  return *this;
}
 
 
/*! @brief Operator multiplying all elements of a matrix by a scalar.
 *
 *     Operation: A = x * A
 *
 * @param (double x) the scaling parameter
 * @return (Matrice_Morse&) the result of the operation (*this)
 */
Matrice_Morse& Matrice_Morse::operator *=(double x )
{
  scale( x );
  return(*this);
}
 
void Matrice_Morse::scale( const double x )
{
  coeff_ *= x;
}
 
void Matrice_Morse::get_stencil( Stencil& stencil ) const
{
  assert_check_morse_matrix_structure( );
 
  if( is_stencil_up_to_date_ )
    {
      stencil = stencil_;
      return;
    }
 
 
  stencil.resize( 0, 2 );
  auto nnz = tab2_.size_array();
  stencil.resize(nnz, 2);
 
 
  ArrOfInt tmp;
 
 
  decltype(nnz) compteur = 0;
 
  const int nb_lines = nb_lignes( );
  for ( int i=0; i<nb_lines; ++i )
    {
      auto k0   = tab1_( i ) - 1;
      auto k1   = tab1_( i + 1 ) - 1;
      const auto size = k1 - k0;
      const int  size_int = (int)size;
 
      tmp.resize_array( 0 );
      tmp.resize_array( size_int );
 
      for ( int k=0; k<size_int; ++k )
        {
          tmp[ k ] = tab2_( k + k0 ) - 1;
        }
 
      tmp.ordonne_array( );
 
      for ( int k=0; k<size_int; ++k )
        {
          stencil( k+compteur , 0 ) = i;
          stencil( k+compteur , 1 ) =  tmp[ k ];
        }
      compteur += size;
    }
 
 
}
 
// Local template method : copy either value or ptr to value!
namespace
{
template<typename _T_> static inline void _fill_slot(_T_& dest, const double& src);
template<> inline void _fill_slot<double>(double& dest, const double& src)
{
  dest = src;
}
template<> inline void _fill_slot<const double *>(const double*& dest, const double& src)
{
  dest = &src;
}
}
 
 
template<typename _TAB_T_, typename _VALUE_T_>
inline void Matrice_Morse::get_stencil_coeff_templ( Stencil& stencil, _TAB_T_& coeffs_span) const
{
  auto nnz = tab2_.size_array();
  coeffs_span.resize(nnz);
  stencil.resize(nnz, 2);
  decltype(nnz) compteur = 0;
  const int nb_lines = nb_lignes( );
  for ( int i=0; i<nb_lines; ++i )
    {
      const auto k0      = tab1_( i ) - 1;
      const auto k1      = tab1_( i + 1 ) - 1;
      const int  size_int = (int)(k1 - k0);
      for ( int k=0; k<size_int; ++k )
        {
          stencil( compteur + k , 0 ) = i;
          stencil( compteur + k , 1 ) = tab2_( k + k0 ) - 1;
          ::_fill_slot<_VALUE_T_>(coeffs_span[ compteur + k ], coeff_(k+k0));
        }
      compteur += size_int;
    }
 
 
}
 
 
void Matrice_Morse::get_stencil_and_coeff_ptrs(Stencil& stencil,
                                               std::vector<const double *>& coeff_ptr) const
{
  assert_check_morse_matrix_structure( );
 
  if( is_stencil_up_to_date_ )
    {
      Cerr << "Error in Matrice_Morse::get_symmetric_stencil_and_coeff_ptrs( )"<<finl;
      Cerr << "  stencil up to date - function not impl. in this case."<<finl;
      Cerr << "  Aborting..." << finl;
      Process::abort( );
      return;
    }
 
  get_stencil_coeff_templ< std::vector<const double *>, const double *>(stencil, coeff_ptr);
  assert( (trustIdType)coeff_ptr.size( ) == stencil.dimension( 0 ));
}
 
 
void Matrice_Morse::get_stencil_and_coefficients( Stencil&       stencil,
                                                  StencilCoeffs& coefficients ) const
{
  if( is_stencil_up_to_date_ )
    {
      if( coeff_.size( ) == 0 )
        {
          Cerr << "Error in Matrice_Morse::get_stencil_and_coefficients( )"<<finl;
          Cerr << "  The coefficients are not filled."<<finl;
          Cerr << "  Aborting..." << finl;
          Process::abort( );
        }
      stencil = stencil_ ;
      { const auto sz = coeff_.size_array(); coefficients.resize(sz); for (auto k=sz-sz; k<sz; k++) coefficients[k] = coeff_[k]; }
      return;
    }
 
 
 
  get_stencil_coeff_templ<StencilCoeffs, double>(stencil, coefficients);
  assert( coefficients.size_array( ) == stencil.dimension( 0 ));
}
 
 
/*! @brief Operator dividing all elements of a matrix by a scalar.
 *
 *     Operation: A =  A / x
 *
 * @param (double x) the scaling parameter
 * @return (Matrice_Morse&) the result of the operation (*this)
 * @throws division by zero not possible
 */
Matrice_Morse& Matrice_Morse::operator /=(double x )
{
  coeff_/=x;
  return(*this);
}
 
void Matrice_Morse::remplir(const IntLists& voisins,
                            const DoubleLists& valeurs,
                            const DoubleVect& terme_diag)
{
 
  int num_elem;
  int compteur,rang =0;
 
  // Fill arrays tab1, tab2 and coeff_:
  auto* p_tab1 = tab1_.addr();
  int* p_tab2 = tab2_.addr();
  double* p_coeff = coeff_.addr();
 
  int* tab2_ptr = p_tab2;
  int n=nb_lignes();
 
  for (num_elem=0; num_elem<n; num_elem++)
    {
 
      IntList_Curseur liste_vois(voisins[num_elem]);
      DoubleList_Curseur liste_val(valeurs[num_elem]);
      compteur =0;
      *p_tab1++ = rang;
 
      *tab2_ptr++=num_elem;
      *p_coeff++ = terme_diag[num_elem];
 
      while  (liste_vois)
        {
          *tab2_ptr++ = liste_vois.valeur();
          *p_coeff++ = liste_val.valeur();
          ++liste_vois;
          ++liste_val;
          compteur++;
        }
      //       tab2[rang] = compteur;
      rang += (compteur + 1);
    }
  tab1_(num_elem)=rang;
  formeF();
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  is_stencil_up_to_date_=false;
}
 
void Matrice_Morse::remplir(const IntLists& voisins,
                            const DoubleLists& valeurs)
{
 
  int num_elem;
  int compteur,rang =0;
 
  // Fill arrays tab1, tab2 and coeff_:
  auto* p_tab1 = tab1_.addr();
  int* p_tab2 = tab2_.addr();
  double* p_coeff = coeff_.addr();
 
  int* tab2_ptr = p_tab2;
  int n=nb_lignes();
 
  for (num_elem=0; num_elem<n; num_elem++)
    {
 
      IntList_Curseur liste_vois(voisins[num_elem]);
      DoubleList_Curseur liste_val(valeurs[num_elem]);
      compteur =0;
      *p_tab1++ = rang;
 
      while  (liste_vois)
        {
          *tab2_ptr++ = liste_vois.valeur();
          *p_coeff++ = liste_val.valeur();
          ++liste_vois;
          ++liste_val;
          compteur++;
        }
      //       tab2[rang] = compteur;
      rang += (compteur);
    }
  tab1_(num_elem)=rang;
  formeF();
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  is_stencil_up_to_date_=false;
}
 
/*! @brief Fill a Morse matrix from a smaller Morse matrix
 *
 */
void Matrice_Morse::remplir(const int ideb, const int jdeb, const int n, const int m, const Matrice_Morse& mat)
{
  // Verification
  assert(ideb<=n);
  assert(jdeb<=m);
 
  // Build a local matrix
  Matrice_Morse matrice_locale(mat);
  // Case where the local matrix is symmetric
  if (sub_type(Matrice_Morse_Sym,mat))
    {
      // Create the lower part L
      Matrice_Morse L(matrice_locale);
      L.transpose(matrice_locale);
      int lordre = L.ordre();
      for (int i=0; i<lordre; i++)
        L(i, i) = 0.;
      // Add M=U+L
      matrice_locale += L;
    }
 
  // Size the global matrix
  auto nnz=matrice_locale.nb_coeff();
  dimensionner(n,m,(int)nnz);
 
  // Fill the global matrix from the local matrix:
  // Fill tab1_ with offset ideb:
  int mon_nb_lignes=matrice_locale.nb_lignes();
  assert(mon_nb_lignes+ideb<=n);
  for (int i=0; i<ideb; i++)
    tab1_(i)=1;
  for (int i=0; i<mon_nb_lignes; i++)
    tab1_(i+ideb)=matrice_locale.tab1_(i);
  for (int i=mon_nb_lignes+ideb; i<n+1; i++)
    tab1_(i)=matrice_locale.tab1_(mon_nb_lignes);
 
  // Fill tab2_ with offset jdeb:
  for (auto i=0; i<nnz; i++)
    tab2_(i)=matrice_locale.tab2_(i)+jdeb;
 
  // Fill coeff_:
  for (auto i=0; i<nnz; i++)
    coeff_(i)=matrice_locale.coeff_(i);
 
  morse_matrix_structure_has_changed_=1, sorted_ = 0;
  is_stencil_up_to_date_=false;
}
 
void Matrice_Morse::formeC()
{
  int n=nb_lignes();
  for(int ii=0; ii<=n; ii++)
    tab1_(ii)--;
  for(int ii=0; ii<n; ii++)
    tab2_(tab1_(ii))=nb_vois(ii);
  for(auto k=0; k<nb_coeff(); k++)
    tab2_(k)--;
  morse_matrix_structure_has_changed_=1;
  is_stencil_up_to_date_=false;
}
 
void Matrice_Morse::formeF()
{
  int n=nb_lignes();
  for(int ii=0; ii<=n; ii++)
    tab1_(ii)++;
  for(auto k=0; k<nb_coeff(); k++)
    tab2_(k)++;
  morse_matrix_structure_has_changed_=1;
  is_stencil_up_to_date_=false;
}
 
 
 
 
/*! @brief DOES NOTHING
 *
 * @return (int) always returns 1
 */
int Matrice_Morse_test()
{
  return 1;
}
 
/*! @brief Remplit la matrice avec des zeros.
 *
 */
 
void Matrice_Morse::clean()
{
  coeff_ = 0;
}
 
/*! @brief Computes the bandwidth of a Morse matrix
 *
 */
int Matrice_Morse::largeur_de_bande() const
{
  int ldist,min = 0;
  const auto* p_tab1_ = get_tab1().addr();
  const int* p_tab2_ = get_tab2().addr();
  int N=ordre();
 
  for(int i=0; i<N; i++)
    for(auto k = p_tab1_[i]; k < p_tab1_[i+1]; k++)
      {
        if (p_tab2_[k-1]-1<N)
          {
            ldist = p_tab2_[k-1] - i;
            if( min < ldist ) min = ldist;
          }
      };
  return min;
}
 
bool Matrice_Morse::check_morse_matrix_structure() const
{
  const int nb_lines   = nb_lignes( );
  const int nb_columns = nb_colonnes( );
  const auto nb_coefficients = tab1_( nb_lines ) - 1;
 
  if ( tab2_.size_array( ) != nb_coefficients )
    {
      Cerr << "invalid tab2 size" << finl;
      return false;
    }
 
  if ( coeff_.size_array( ) != nb_coefficients )
    {
      Cerr << "invalid coeff size" << finl;
      return false;
    }
 
  ArrOfBit flags( nb_columns );
 
  for ( int i=0; i<nb_lines; ++i )
    {
      flags = 0;
 
      auto k0 = tab1_( i ) - 1;
      auto k1 = tab1_( i + 1 ) - 1;
 
      for ( auto k=k0; k<k1; ++k )
        {
          int j = tab2_( k ) - 1;
 
          if ( j < 0 )
            {
              Cerr << "invalid column index (<0): " << j << finl;
              return false;
            }
 
          if ( j >= nb_columns )
            {
              Cerr << "invalid column index (>nb_cols): " << j << " > " << nb_columns << finl;
              return false;
            }
 
          if ( flags[ j ] )
            {
              Cerr << "invalid coefficient ( " << i << ", " << j << " ): already defined ( " << k << " )" << finl;
              return false;
            }
 
          flags.setbit( j );
        }
    }
 
  return true;
}
 
bool Matrice_Morse::check_sorted_morse_matrix_structure() const
{
  const int nb_lines   = nb_lignes( );
  const int nb_columns = nb_colonnes( );
  const auto nb_coefficients = tab1_( nb_lines ) - 1;
 
  if ( tab2_.size_array( ) != nb_coefficients )
    {
      Cerr << "invalid tab2 size" << finl;
      return false;
    }
 
  if ( coeff_.size_array( ) != nb_coefficients )
    {
      Cerr << "invalid coeff size" << finl;
      return false;
    }
 
  ArrOfBit flags( nb_columns );
 
  for ( int i=0; i<nb_lines; ++i )
    {
      flags = 0;
 
      auto k0 = tab1_( i ) - 1;
      auto k1 = tab1_( i + 1 ) - 1;
 
      int j0 = tab2_( k0 ) - 1 - 1;
 
      for ( auto k=k0; k<k1; ++k )
        {
          int j = tab2_( k ) - 1;
 
          if ( j < 0 )
            {
              Cerr << "invalid column index (<0): " << j << finl;
              return false;
            }
 
          if ( j >= nb_columns )
            {
              Cerr << "invalid column index (>nb_cols): " << j << " > " << nb_columns << finl;
              return false;
            }
 
          if ( flags[ j ] )
            {
              Cerr << "invalid coefficient ( " << i << ", " << j << " ): already defined ( " << k << " )" << finl;
              return false;
            }
 
          if ( j <= j0 )
            {
              Cerr << "unsorted coefficient: ( " << i << ", " << j << " ) after ( " << i << ", " << j0 << " ) " << finl;;
              return false;
            }
 
          j0 = j;
          flags.setbit( j );
        }
    }
 
  return true;
}
 
 
void Matrice_Morse::assert_check_morse_matrix_structure() const
{
  if (!morse_matrix_structure_has_changed_) return;
#ifndef NDEBUG
  if ( ! ( check_morse_matrix_structure( ) ) )
    {
      Cerr << "Error in 'Matrice_Morse::assert_check_morse_matrix_structure( )':" << finl;
      Cerr << "  Exiting..." << finl;
      Process::exit( );
    }
  else
    morse_matrix_structure_has_changed_=0;
#endif
}
 
void Matrice_Morse::assert_check_sorted_morse_matrix_structure() const
{
  if (!morse_matrix_structure_has_changed_) return;
#ifndef NDEBUG
  if ( ! ( check_sorted_morse_matrix_structure( ) ) )
    {
      Cerr << "Error in 'Matrice_Morse::assert_check_sorted_morse_matrix_structure( )':" << finl;
      Cerr << "  Exiting..." << finl;
      Process::exit( );
    }
  else
    morse_matrix_structure_has_changed_=0;
#endif
}
 
// Build a new Morse matrix spanning the rectangular area defined by the two points (nl0, nc0) and (nl1, nc1)
// in the original matrix.
// Indices are provided in C mode (0-based indexing).
void Matrice_Morse::construire_sous_bloc(int nl0, int nc0, int nl1, int nc1, Matrice_Morse& result) const
{
  // count non-zero entries:
  assert(nl0 >= 0);
  assert(nc0 >= 0);
  assert(nl0 <= nl1);
  assert(nc0 <= nc1);
 
  auto max_nnz = tab1_(nl1+1) - tab1_(nl0); // maximum number of zeros that we will find
  int tot=0;
  IntTab loca((int)max_nnz, 2);
  DoubleTab sub_coeffs((int)max_nnz);
  for (int li=nl0; li <= nl1; li++)
    {
      auto idx_coeff = tab1_(li)-1;
      int nb_coeff_on_line = (int)(tab1_(li+1)-tab1_(li));
      for (int j=0; j < nb_coeff_on_line; j++)
        {
          int col_idx = tab2_(j+idx_coeff)-1;
          if (col_idx >= nc0 && col_idx <= nc1) // is the coeff in the window?
            {
              loca(tot, 0) = li - nl0;
              loca(tot, 1) = col_idx - nc0;
              sub_coeffs(tot) = coeff_(j+idx_coeff);
              tot++;
            }
        }
    }
  loca.resize(tot,2);
  sub_coeffs.resize(tot);
 
  result.dimensionner(loca);
  // Set coefficient values:
  for (int i =0 ; i < tot; i++)
    {
      int il = loca(i, 0);
      int ic = loca(i, 1);
      result.coef(il, ic) = sub_coeffs(i);
    }
}
 
void Matrice_Morse::sort_stencil()
{
  if (sorted_) return; //already done
  for (int i = 0; i + 1 < tab1_.size_array(); i++) //indice de ligne
    std::sort(tab2_.addr() + tab1_(i) - 1, tab2_.addr() + tab1_(i + 1) - 1);
  morse_matrix_structure_has_changed_ = sorted_ = 1;
}
 
// Check if the matrix is sorted based on a stencil condition
bool Matrice_Morse::is_sorted_stencil() const
{
  if (!sorted_)
    {
      const int n = nb_lignes();
      for (int i = 0; i < n; i++)
        {
          const auto k0 = tab1_( i ) - 1;
          const auto k1 = tab1_( i + 1 ) - 1;
          for (auto k=k0; k<k1-1; k++)
            if (tab2_(k)>tab2_(k+1))
              return sorted_; // not sorted
        }
      sorted_ = true;
    }
  return sorted_;
}
 
// Check the matrix is diagonal:
// Faster than using:
// Stencil stencil;
// A.get_stencil(stencil);
// Matrix_tools::is_diagonal_stencil(A.nb_lignes(), A.nb_colonnes(), stencil);
bool Matrice_Morse::is_diagonal()
{
  bool is_diagonal = true;
  const int n = nb_lignes();
  for (int i = 0; i < n; i++)
    {
      const auto k1 = get_tab1()(i) - 1;
      const auto k2 = get_tab1()(i + 1) - 1;
      for (auto k = k1; k < k2; k++)
        {
          if (k2-k1>1 || get_tab2()(k)-1!=i)
            {
              is_diagonal = false;
              break;
            }
        }
    }
  return is_diagonal;
}
 
// Explicit instantiations for 'auto nnz' abbreviated function templates
template Matrice_Morse::Matrice_Morse(int, int);
template Matrice_Morse::Matrice_Morse(int, int, int);
template void Matrice_Morse::dimensionner(int, int);
template void Matrice_Morse::dimensionner(int, int, int);
#ifdef TRUST_USE_GPU
template Matrice_Morse::Matrice_Morse(int, trustIdType);
template Matrice_Morse::Matrice_Morse(int, int, trustIdType);
template void Matrice_Morse::dimensionner(int, trustIdType);
template void Matrice_Morse::dimensionner(int, int, trustIdType);
#endif