NormalizeOpt.h

// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
//-----------------------------------------------------------------------------
// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: Apache-2.0
//-----------------------------------------------------------------------------
#pragma once
 
#include <alien/data/CompositeMatrix.h>
#include <alien/data/CompositeVector.h>
#include <alien/data/IMatrix.h>
#include <alien/data/IVector.h>
#include <alien/kernels/composite/CompositeMatrix.h>
#include <alien/kernels/composite/CompositeVector.h>
#include <alien/kernels/simple_csr/SimpleCSRMatrix.h>
#include <alien/kernels/simple_csr/SimpleCSRVector.h>
#include <alien/utils/Precomp.h>
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
namespace Alien
{
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

class ALIEN_EXPORT NormalizeOpt
{
 public:
  enum eAlgoType
  {
    StdLU,
    EigenLU
  };

  enum eOptType
  {
    SumFirstEq
  };

  enum eErrorType
  {
    NoError,
    PivotError,
    WithErrors
  };

  typedef SimpleCSRMatrix<Arccore::Real> MatrixImpl;
  typedef SimpleCSRVector<Arccore::Real> VectorImpl;

  NormalizeOpt();

  virtual ~NormalizeOpt() {}

  void setAlgo(eAlgoType algo);

  void setOpt(eOptType opt, bool flag);

  eErrorType normalize(IMatrix& matrix, IVector& vector) const;

  eErrorType normalize(CompositeMatrix& matrix, CompositeVector& vector,
                       ConstArrayView<Integer> eq_ids) const;
 
  template <int N, bool check_null_pivot = false>
  class LU;
  class Op;
  class Op2;
 
 private:
  eErrorType _normalize(MatrixImpl& matrix, VectorImpl& vector) const;

  eErrorType _normalize(MatrixImpl& matrix, MatrixImpl& matrix2, bool trans,
                        ConstArrayView<Integer> eq_ids, VectorImpl& vector) const;

  eErrorType _normalize(MatrixImpl& matrix,
                        Arccore::ConstArrayView<Arccore::Integer> eq_ids, MatrixImpl& matrix2, bool trans,
                        VectorImpl& vector) const;

  eAlgoType m_algo;
  bool m_sum_first_eq;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

template <int N, bool check_null_pivot>
class NormalizeOpt::LU
{
 public:
  typedef Arccore::Real Block2DType[N][N];
  /*
   * \brief Constructor
   * \param[in] Ap The matrix
   * \param[in] factorize Whether or not the matrix should be factorized
   */
  LU(Arccore::Real* Ap, bool factorize = true)
  : m_A(*(Block2DType*)Ap)
  {
    if (factorize) {
      Block2DType& A = *(Block2DType*)Ap;
 
      for (int k = 0; k < N; ++k) {
        if (check_null_pivot) {
          assert(A[k][k] != 0);
          A[k][k] = 1 / A[k][k];
        }
        else {
          if (A[k][k] == 0)
            A[k][k] = 0;
          else
            A[k][k] = 1 / A[k][k];
        }
        for (int i = k + 1; i < N; ++i) {
          A[i][k] *= A[k][k];
          // TOCHECK : to be removed ?
          // A[i][k] /= A[k][k];
        }
        for (int i = k + 1; i < N; ++i) {
          for (int j = k + 1; j < N; ++j) {
            A[i][j] -= A[i][k] * A[k][j];
          }
        }
      }
    }
  }

  void setIdentity()
  {
    for (int i = 0; i < N; ++i) {
      for (int j = 0; j < i; ++j) {
        m_A[i][j] = 0;
      }
      m_A[i][i] = 1;
      for (int j = i + 1; j < N; ++j) {
        m_A[i][j] = 0;
      }
    }
  }

  void setZero()
  {
    Real* ptr = (Real*)m_A;
    for (int i = 0; i < N * N; ++i) {
      ptr[i] = 0;
    }
  }

  template <int NRhs, bool TransRhs>
  void LXSolve(Arccore::Real* Xp) const
  {
    if (TransRhs) {
      typedef Real Rhs2DType[NRhs][N];
      Rhs2DType& X = *(Rhs2DType*)Xp;
 
      for (int k = 0; k < NRhs; ++k) {
        for (int i = 1; i < N; ++i) {
          for (int j = 0; j < i; ++j) {
            X[k][i] -= m_A[i][j] * X[k][j];
          }
        }
      }
    }
    else {
      typedef Real Rhs2DType[N][NRhs];
      Rhs2DType& X = *(Rhs2DType*)Xp;
 
      for (int i = 1; i < N; ++i) {
        for (int j = 0; j < i; ++j) {
          for (int k = 0; k < NRhs; ++k) {
            X[i][k] -= m_A[i][j] * X[j][k];
          }
        }
      }
    }
  }

  template <int NRhs, bool TransRhs>
  void UXSolve(Arccore::Real* Xp) const
  {
    if (TransRhs) {
      typedef Real Rhs2DType[NRhs][N];
      Rhs2DType& X = *(Rhs2DType*)Xp;
      for (int i = N - 1; i >= 0; --i) {
        for (int k = NRhs - 1; k >= 0; --k) {
          for (int j = N - 1; j > i; --j) {
            X[k][i] -= m_A[i][j] * X[k][j];
          }
        }
        for (int k = NRhs - 1; k >= 0; --k) {
          // TOCHECK : to be removed ?
          // X[k][i] /= m_A[i][i];
          X[k][i] *= m_A[i][i];
        }
      }
    }
    else {
      typedef Real Rhs2DType[N][NRhs];
      Rhs2DType& X = *(Rhs2DType*)Xp;
 
      for (int i = N - 1; i >= 0; --i) {
        for (int j = N - 1; j > i; --j) {
          for (int k = NRhs - 1; k >= 0; --k) {
            X[i][k] -= m_A[i][j] * X[j][k];
          }
        }
        for (int k = NRhs - 1; k >= 0; --k) {
          X[i][k] *= m_A[i][i];
          // TOCHECK: to be removed ?
          // X[i][k] /= m_A[i][i];
        }
      }
    }
  }

  template <int NRhs, bool TransRhs>
  void solve(Arccore::Real* Xp) const
  {
    LXSolve<NRhs, TransRhs>(Xp);
    UXSolve<NRhs, TransRhs>(Xp);
  }
 
 private:
  Block2DType& m_A;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

class NormalizeOpt::Op
{
 public:
  Op(MatrixImpl& m, VectorImpl& x, eAlgoType algo, bool sum_first_eq)
  : m_local_size(m.getCSRProfile().getNRow())
  , m_local_offset(m.getLocalOffset())
  , m_cols(m.getCSRProfile().getCols())
  , m_row_offset(m.getCSRProfile().getRowOffset())
  , m_matrix(m.internal()->getDataPtr())
  , m_rhs(x.getDataPtr())
  , m_sum_first_eq(sum_first_eq)
  , m_algo(algo)
  {
    // TOCHECK : to be removed ?
    // m_equations_num = m.space().structInfo().size() ;
    m_equations_num = 1;
    if (m.block() && m.block()->size())
      m_equations_num = m.block()->size();
    m_unknowns_num = m_equations_num;
    m_block_size = m_equations_num * m_unknowns_num;
    m_diag_first = m.getCSRProfile().getDiagFirstOpt();
    if (!m_diag_first)
    {
      if(m.isParallel())
        m_upper_diag_offset = m.getDistStructInfo().getUpperDiagOffset(m.getCSRProfile());
      else
        m_upper_diag_offset = m.getCSRProfile().getUpperDiagOffset().constView();
    }
    m_keep_diag = true;
    if (m_sum_first_eq) {
      if (m_diag_first)
        sumBlockEq<true>();
      else
        sumBlockEq<false>();
    }
  }

  void checkNullEq(Arccore::Integer irow, Arccore::Integer diag_offset)
  {
    // search line with only zero
    for (Integer ieq = 0; ieq < m_equations_num; ++ieq) {
      bool ok = false;
      Integer col = m_row_offset[irow] + diag_offset;
      for (Integer ui = 0; ui < m_unknowns_num; ++ui)
        if (m_matrix[col * m_block_size + ij(ieq, ui)] != 0) {
          ok = true;
          break;
        }
      if (!ok) {
        // search null column
        for (Arccore::Integer ui = 0; ui < m_unknowns_num; ++ui) {
          bool ok2 = false;
          for (Integer jeq = 0; jeq < m_equations_num; ++jeq)
            if (m_matrix[col * m_block_size + ij(jeq, ui)] != 0) {
              ok2 = true;
              break;
            }
          if (!ok2) {
            // put 1 on ui column
            m_matrix[col * m_block_size + ij(ieq, ui)] = 1.;
            break;
          }
        }
      }
    }
  }

  template <bool diag_first>
  void sumBlockEq()
  {
    for (Integer irow = 0; irow < m_local_size; ++irow) {
      checkNullEq(irow, (diag_first ? 0 : m_upper_diag_offset[irow] - m_row_offset[irow]));
      for (Integer col = m_row_offset[irow]; col < m_row_offset[irow + 1]; ++col) {
        for (Integer ieq = 1; ieq < m_equations_num; ++ieq)
          for (Integer ui = 0; ui < m_unknowns_num; ++ui)
            m_matrix[col * m_block_size + ij(0, ui)] +=
            m_matrix[col * m_block_size + ij(ieq, ui)];
      }
      for (Integer ieq = 1; ieq < m_equations_num; ++ieq)
        m_rhs[irow * m_equations_num] += m_rhs[irow * m_equations_num + ieq];
    }
  }

  template <int N>
  void multInvDiag();
 
 protected:
  Integer ijk(Integer i, Integer j, Integer k) const
  {
    return k * m_block_size + i * m_unknowns_num + j;
  }
 
  Arccore::Integer ij(Arccore::Integer i, Arccore::Integer j) const
  {
    return i * m_unknowns_num + j;
  }
 
  Arccore::Integer ik(Arccore::Integer i, Arccore::Integer k) const
  {
    return k * m_equations_num + i;
  }
 
  Arccore::Integer jk(Arccore::Integer j, Arccore::Integer k) const
  {
    return k * m_unknowns_num + j;
  }

  Arccore::Integer m_local_size;
  Arccore::Integer m_local_offset;
  Arccore::ConstArrayView<Arccore::Integer> m_cols;
  Arccore::ConstArrayView<Arccore::Integer> m_row_offset;
  Arccore::ConstArrayView<Arccore::Integer> m_upper_diag_offset;
  Arccore::Integer m_equations_num;
  Arccore::Integer m_unknowns_num;
  Arccore::Integer m_block_size;
  Arccore::Real* m_matrix;
  Arccore::Real* m_rhs;
  bool m_keep_diag;
  bool m_diag_first;
  bool m_sum_first_eq;
  eAlgoType m_algo;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
// TOCHECK : to be removed ?
/*
  class NormalizeOpt::Op
  {
  public :
  Op(MatrixImpl& m, VectorImpl& x,eAlgoType algo,bool sum_first_eq)
  : m_local_size(m.getCSRProfile().getNRow())
  , m_local_offset(m.getLocalOffset())
  , m_cols(m.getCSRProfile().getCols())
  , m_row_offset(m.getCSRProfile().getRowOffset())
  , m_matrix(m.internal().getDataPtr())
  , m_rhs(x.getDataPtr())
  , m_sum_first_eq(sum_first_eq)
  , m_algo(algo)
  {
  m_equations_num = m.space().structInfo().size() ;
  m_unknowns_num = m_equations_num ;
  m_block_size = m_equations_num*m_unknowns_num ;
  m_diag_first = m.getCSRProfile().getDiagFirstOpt() ;
  if(!m_diag_first)
  m_upper_diag_offset = m.getCSRProfile().getUpperDiagOffset() ;
  m_keep_diag = true ;
  if(m_sum_first_eq)
  {
  if(m_diag_first)
  sumBlockEq<true>() ;
  else
  sumBlockEq<false>() ;
  }
  }
 
  void checkNullEq(Integer irow,Integer diag_offset)
  {
  //search line with only zero
  for(Integer ieq=0;ieq<m_equations_num;++ieq)
  {
  bool ok = false ;
  Integer col = m_row_offset[irow]+diag_offset ;
  for(Integer ui=0;ui<m_unknowns_num;++ui)
  if(m_matrix[col*m_block_size+ij(ieq,ui)]!=0)
  {
  ok = true ;
  break ;
  }
  if(!ok)
  {
  //search null colonne
  for(Integer ui=0;ui<m_unknowns_num;++ui)
  {
  bool ok2 = false ;
  for(Integer jeq=0;jeq<m_equations_num;++jeq)
  if(m_matrix[col*m_block_size+ij(jeq,ui)]!=0)
  {
  ok2 = true ;
  break ;
  }
  if(!ok2)
  {
  //put 1 on ui colonne
  m_matrix[col*m_block_size+ij(ieq,ui)] = 1. ;
  break ;
  }
  }
  }
  }
  }
 
  template<bool diag_first>
  void
  sumBlockEq()
  {
  for(Integer irow=0;irow<m_local_size;++irow)
  {
  checkNullEq(irow,(diag_first?0:m_upper_diag_offset[irow]-m_row_offset[irow])) ;
  for(Integer col= m_row_offset[irow];col<m_row_offset[irow+1];++col)
  {
  for(Integer ieq=1;ieq<m_equations_num;++ieq)
  for(Integer ui=0;ui<m_unknowns_num;++ui)
  m_matrix[col*m_block_size+ij(0,ui)] += m_matrix[col*m_block_size+ij(ieq,ui)] ;
  }
  for(Integer ieq=1;ieq<m_equations_num;++ieq)
  m_rhs[irow*m_equations_num] += m_rhs[irow*m_equations_num+ieq] ;
  }
  }
 
  template<int N>
  void multInvDiag() ;
  protected :
  Integer ijk(Integer i, Integer j, Integer k) const {
  return k*m_block_size+i*m_unknowns_num+j ;
  }
  Integer ij(Integer i, Integer j) const {
  return i*m_unknowns_num+j ;
  }
  Integer ik(Integer i, Integer k) const {
  return k*m_equations_num+i ;
  }
  Integer jk(Integer j, Integer k) const {
  return k*m_unknowns_num+j ;
  }
 
  Integer                 m_local_size ;
  Integer                 m_local_offset ;
  ConstArrayView<Integer> m_cols  ;
  ConstArrayView<Integer> m_row_offset ;
  ConstArrayView<Integer> m_upper_diag_offset ;
  Integer                 m_equations_num ;
  Integer                 m_unknowns_num ;
  Integer                 m_block_size ;
  Real*                   m_matrix ;
  Real*                   m_rhs ;
  bool                    m_keep_diag ;
  bool                    m_diag_first ;
  bool                    m_sum_first_eq;
  eAlgoType               m_algo ;
  } ;
*/

class NormalizeOpt::Op2 : public NormalizeOpt::Op
{
 public:
  Op2(MatrixImpl& m, MatrixImpl& m2, bool trans,
      Arccore::ConstArrayView<Arccore::Integer> eq_ids, VectorImpl& x, eAlgoType algo,
      bool sum_first_eq)
  : Op(m, x, algo, sum_first_eq)
  , m_nb_extra_eq(m2.getCSRProfile().getNRow())
  , m_eq_ids(eq_ids)
  , m_submatrix1(false)
  , m_submatrix2(true)
  , m_trans(trans)
  , m_nuk2(0)
  , m_extra_eq_cols(m2.getCSRProfile().getCols())
  , m_extra_eq_row_offset(m2.getCSRProfile().getRowOffset())
  , m_extra_eq_matrix(m2.internal()->getDataPtr())
  {
    if (m_sum_first_eq) {
      if (m_trans) {
        // TOCHECK : to be removed ?
        // Integer neq = m2.space().structInfo().size() ;
        Integer neq = 1;
        if (m2.block() && m2.block()->size())
          neq = m2.block()->size();
        for (Integer i = 0; i < m_eq_ids.size(); ++i) {
          Integer ieq = m_eq_ids[i];
          for (Integer k = m_extra_eq_row_offset[ieq]; k < m_extra_eq_row_offset[ieq + 1];
               ++k) {
            Real val = m_extra_eq_matrix[k * neq];
            for (Integer ieq = 1; ieq < neq; ++ieq)
              val += m_extra_eq_matrix[k * neq + ieq];
            m_extra_eq_matrix[k * neq] = val;
          }
        }
      }
      else {
        // TOCHECK : to be removed ?
        // Integer neq = m2.eqSpace().structInfo().size() ;
        Integer neq = 1;
        if (m2.rowBlock() && m2.rowBlock()->maxBlockSize())
          neq = m2.rowBlock()->maxBlockSize();
        for (Integer i = 0; i < m_eq_ids.size(); ++i) {
          Integer irow = m_eq_ids[i];
          for (Integer k = m_extra_eq_row_offset[irow];
               k < m_extra_eq_row_offset[irow + 1]; ++k) {
            Real val = m_extra_eq_matrix[k * neq];
            for (Integer ieq = 1; ieq < neq; ++ieq)
              val += m_extra_eq_matrix[k * neq + ieq];
            m_extra_eq_matrix[k * neq] = val;
          }
        }
      }
    }
  }

  Op2(MatrixImpl& m, Arccore::ConstArrayView<Arccore::Integer> eq_ids, MatrixImpl& m2,
      bool trans, VectorImpl& x, eAlgoType algo, bool sum_first_eq)
  : Op(m, x, algo, sum_first_eq)
  , m_nb_extra_eq(m2.getCSRProfile().getNRow())
  , m_eq_ids(eq_ids)
  , m_submatrix1(true)
  , m_submatrix2(false)
  , m_trans(trans)
  , m_nuk2(0)
  , m_extra_eq_cols(m2.getCSRProfile().getCols())
  , m_extra_eq_row_offset(m2.getCSRProfile().getRowOffset())
  , m_extra_eq_matrix(m2.internal()->getDataPtr())
  {
    if (m_trans)
      // TOCHECK : to be removed
      // m_nuk2 = m2.eqSpace().structInfo().size() ;
      m_nuk2 = m2.rowBlock()->maxBlockSize();
    else
      // TOCHECK : to be removed
      // m_nuk2 = m2.uSpace().structInfo().size() ;
      m_nuk2 = m2.colBlock()->maxBlockSize();
  }

  template <int N>
  void multInvDiag();
 
 private:
  Arccore::Integer m_nb_extra_eq;
  Arccore::ConstArrayView<Arccore::Integer> m_eq_ids;
  bool m_submatrix1;
  bool m_submatrix2;
  bool m_trans;
  Arccore::Integer m_nuk2;
  Arccore::ConstArrayView<Arccore::Integer> m_extra_eq_cols;
  Arccore::ConstArrayView<Arccore::Integer> m_extra_eq_row_offset;
  Arccore::Real* m_extra_eq_matrix;
};
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
} // namespace Alien
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/