MD_Vector_tools.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 <Schema_Comm_Vecteurs.h>
#include <MD_Vector_composite.h>
#include <Echange_EV_Options.h>
#include <MD_Vector_tools.h>
#include <MD_Vector_seq.h>
#include <communications.h>
#include <Schema_Comm.h>
#include <TRUSTTrav.h>
#include <Perf_counters.h>
#include <vector>
#include <typeindex>
 
Schema_Comm_Vecteurs MD_Vector_tools::comm;
MD_Vector MD_Vector_tools::last_md;
static std::type_index last_type_idx(typeid(void *));  // love this one ;-)
int MD_Vector_tools::last_linesize = 0;
Echange_EV_Options MD_Vector_tools::last_opt;
 
 
class MD_Vector_renumber
{
public:
  // for each item in the source vector, should it be kept
  ArrOfBit items_to_keep_;
  // new C index of each item in the source vector to renumber.
  //  (for items not kept, must point to a geometrically equivalent item)
  // renum_ can be a global or a local number.
  IntVect renum_;
  // number of non-zero items in items_to_keep_:
  int  nb_items_to_keep_;
};
 
namespace
{
 
/** Should we use resize() or resize_dim0 to resize array ...
 */
template <typename _SIZE_, typename _TYPE_>
bool resize_tab_or_vect(TRUSTVect<_TYPE_,_SIZE_>& v, _SIZE_ sz, int sz_r, RESIZE_OPTIONS opt)
{
  using TAB = TRUSTTab<_TYPE_, _SIZE_>;
  TAB* vv = dynamic_cast<TAB*>(&v);
  if (vv)
    {
      TAB& t = *vv;
      const _SIZE_ n = t.dimension_tot(0);
      if (n == sz_r || n == 0) t.resize_dim0(sz, opt);
      else if (n == sz) { /* ok no resize */  }
      else return true;
    }
  else
    {
      const _SIZE_ n = v.size_totale();
      if (n == sz_r || n == 0) {  v.resize(sz, opt); }
      else if (n == sz) { /* ok no resize */ }
      else return true;
    }
  return false;
}
 
} // end anonymous NS
 
template <class VECT, class TAB>
static void creer_tableau_distribue_(const MD_Vector& md, VECT& v, RESIZE_OPTIONS opt)
{
  if (v.get_md_vector())
    {
      // If this message appears and you know what you are doing,
      // do this before calling creer_tableau_distribue:
      //  vect.set_md_vector(MD_Vector()); // Cancels the descriptor and keeps the vector.
      // or
      //  vect.reset(); // Completely destroys the vector contents
      Cerr << "Internal error in MD_Vector_tools::creer_tableau_distribue:\n"
           << " Vector already has a parallel vector structure" << finl;
      Process::exit();
    }
  int sz = md->get_nb_items_tot();
  // Note: sz_r may be -1 in some cases. Then the test n==sz_r will always be false,
  //  but that is the intended behaviour...
  int sz_r = md->get_nb_items_reels();
  bool err = ::resize_tab_or_vect(v, sz, sz_r, opt);
 
  if (err)
    {
      Cerr << "Internal error in MD_Vector_tools::creer_tableau_distribue:\n"
           << " Input vector has wrong size or dimension(0): expected 0 or size_reele=" << sz_r << finl;
      Process::exit();
    }
  v.set_md_vector(md);
}
 
template <typename _TYPE_, typename _SIZE_>
static void creer_tableau_seq_(const MD_Vector& md, TRUSTVect<_TYPE_,_SIZE_>& v, RESIZE_OPTIONS opt)
{
  const MD_Vector_seq * md_seq = dynamic_cast<const MD_Vector_seq *>(&md.valeur());
  assert(md_seq != nullptr); // sequential, cast should always be OK
 
  // _SIZE_ might be int, but md_seq->get_nb_items is always trustIdType:
  trustIdType sz0 = md_seq->get_nb_items();
  _SIZE_ sz = std::is_same<_SIZE_, int>::value ? Process::check_int_overflow(sz0) : sz0;
 
  bool err = ::resize_tab_or_vect(v, sz, -1, opt);
  if (err)
    {
      Cerr << "Internal error in MD_Vector_tools::creer_tableau_seq_:\n"
           << " Input vector has wrong size or dimension(0)" << finl;
      Process::exit();
    }
  // Important - keep the same MD_Vector object:
  v.set_md_vector(md);
}
 
 
/*! @brief Transforms v into a parallel array having the structure md.
 *
 * md must be non-null!
 *   The dimension(i>=1) of array v (if it is an IntTab or DoubleTab) are preserved,
 *   dimension(0) and dimension_tot(0) are modified according to the number of items
 *   specified in the MD_Vector.
 *   The initial size of the vector must be either 0, md.get_nb_items_reels(),
 *    or md.get_nb_items_tot(). If needed, the array size is modified and the
 *    array is initialized according to opt.
 *  WARNING, virtual_exchange() is NOT called. The virtual slots are not initialised...
 */
void MD_Vector_tools::creer_tableau_distribue(const MD_Vector& md, Array_base& v, RESIZE_OPTIONS opt)
{
  if (!md)
    {
      Cerr << "Error in MD_Vector_tools::creer_tableau_distribue(): MD_Vector is null" << finl;
      Process::exit();
    }
 
  IntVect* intV = dynamic_cast<IntVect*>(&v);
  DoubleVect* doubleV = dynamic_cast<DoubleVect*>(&v);
  FloatVect* floatV = dynamic_cast<FloatVect*>(&v);
#if INT_is_64_ == 2
  BigIntVect* bintV = dynamic_cast<BigIntVect*>(&v);
  TIDVect* tidV = dynamic_cast<TIDVect*>(&v);
  BigTIDVect* btidV = dynamic_cast<BigTIDVect*>(&v);
  BigDoubleVect* bdoubleV = dynamic_cast<BigDoubleVect*>(&v);
#endif
 
  if(sub_type(MD_Vector_seq, md.valeur()))
    {
      /* [ABN] Here we would like to do
              assert(Process::is_sequential());
       * but we can not: for example when reading a MED file we construct temporarily a dummy domain
       * which looks like a sequential one, even if we are actually running parallel ... too bad.
       */
      if (intV) creer_tableau_seq_<int, int>(md, *intV, opt);
      else if (doubleV) creer_tableau_seq_<double, int>(md, *doubleV, opt);
      else if (floatV) creer_tableau_seq_<float, int>(md, *floatV, opt);
#if INT_is_64_ == 2
      else if (bintV) creer_tableau_seq_<int, trustIdType>(md, *bintV, opt);
      else if (tidV) creer_tableau_seq_<trustIdType, int>(md, *tidV, opt);
      else if (btidV) creer_tableau_seq_<trustIdType, trustIdType>(md, *btidV, opt);
      else if (bdoubleV) creer_tableau_seq_<double, trustIdType>(md, *bdoubleV, opt);
#endif
      else
        {
          Cerr << "Internal error in MD_Vector_tools::creer_tableau_seq():\n"
               << " Invalid array type " << v.que_suis_je()
               << "\n Array must be a subtype of IntVect or DoubleVect or FloatVect" << finl;
          Process::exit();
        }
    }
  else  // parallel computation **or** composite descriptor (MD_Vector_composite)
    {
      if (intV) creer_tableau_distribue_<IntVect, IntTab>(md, *intV, opt);
      else if (doubleV) creer_tableau_distribue_<DoubleVect, DoubleTab>(md, *doubleV, opt);
      else if (floatV) creer_tableau_distribue_<FloatVect, FloatTab>(md, *floatV, opt);
#if INT_is_64_ == 2
      else if (tidV) creer_tableau_distribue_<TIDVect, TIDTab>(md, *tidV, opt);
      else if (bintV || btidV || bdoubleV)
        {
          Cerr << "Internal error in MD_Vector_tools::creer_tableau_distribue(const MD_Vector & md, Array_base & v):" << finl
               << "   -> Trying to create a distributed array for a 'big' array (greater than 32b) - shoud never be necessary." << finl;
          Process::exit();
        }
#endif
      else
        {
          Cerr << "Internal error in MD_Vector_tools::creer_tableau_distribue(const MD_Vector & md, Array_base & v):\n"
               << " Invalid array type " << v.que_suis_je()
               << "\n Array must be a subtype of IntVect or DoubleVect or FloatVect" << finl;
          Process::exit();
        }
    }
}
 
 
template <typename _TYPE_>
void MD_Vector_tools::perform_virtual_exchange(const MD_Vector& md, TRUSTVect<_TYPE_>& v, const Echange_EV_Options& opt, IsExchangeBlocking is_exchange_blocking, const std::string kernel_name)
{
  const MD_Vector_base& mdv = md.valeur();
  const std::type_index type_idx(typeid(_TYPE_));
 
  if ((is_exchange_blocking == IsExchangeBlocking::DefaultBlocking)||(is_exchange_blocking == IsExchangeBlocking::NonBlockingStart))
    {
 
      if (md == last_md && v.line_size() == last_linesize && last_type_idx == type_idx &&  last_opt == opt)
        {
          /* Do nothing if not the first pass, and nothing (including data type) has changed */
        }
      else
        {
          last_md = md;
          last_linesize = v.line_size();
          last_type_idx = type_idx;
          last_opt = opt;
          comm.begin_init();
          mdv.initialize_comm(opt, comm, v);
          comm.end_init();
        }
      bool bufferOnDevice = Process::is_parallel() && v.isDataOnDevice();
      comm.begin_comm(bufferOnDevice);     // buffer allocated on device
      mdv.prepare_send_data(opt, comm, v); // pack buffer on device (read_from_vect_items)
    }
  comm.exchange(is_exchange_blocking, kernel_name);                     // buffer d2h + MPI + buffer h2d
  if ((is_exchange_blocking == IsExchangeBlocking::DefaultBlocking)||(is_exchange_blocking == IsExchangeBlocking::NonBlockingFinish))
    {
      mdv.process_recv_data(opt, comm, v); // unpack buffer on device (write_to_vect_items + write_to_vect_blocs)
      comm.end_comm();
    }
}
 
template<typename _TYPE_>
void MD_Vector_tools::select_virtual_exchange_operation(const MD_Vector& md, TRUSTVect<_TYPE_>& v, MD_Vector_tools::Operations_echange opt, IsExchangeBlocking is_exchange_blocking, const std::string kernel_name)
{
  switch(opt)
    {
    case MD_Vector_tools::ECHANGE_EV:
      perform_virtual_exchange(md, v, echange_ev_opt_default,is_exchange_blocking, kernel_name);
      break;
    case MD_Vector_tools::EV_SOMME:
      perform_virtual_exchange(md, v, Echange_EV_Options(Echange_EV_Options::SUM));
      break;
    case MD_Vector_tools::EV_SOMME_ECHANGE:
      perform_virtual_exchange(md, v, Echange_EV_Options(Echange_EV_Options::SUM));
      perform_virtual_exchange(md, v, echange_ev_opt_default);
      break;
    case MD_Vector_tools::EV_MAX:
      perform_virtual_exchange(md, v, Echange_EV_Options(Echange_EV_Options::MAX));
      break;
    case MD_Vector_tools::EV_MINCOL1:
      perform_virtual_exchange(md, v, Echange_EV_Options(Echange_EV_Options::MINCOL1));
      break;
    default:
      Cerr << "select_virtual_exchange_operation operation not implemented" << finl;
      Process::exit();
    }
}
 
template<typename _TYPE_>
inline void MD_Vector_tools::call_virtual_exchange(TRUSTVect<_TYPE_>& v, MD_Vector_tools::Operations_echange opt, IsExchangeBlocking is_exchange_blocking, const std::string kernel_name)
{
  const MD_Vector& md = v.get_md_vector();
  if (md && Process::is_parallel())
    {
      // [ABN] in some weird cases (like building a temporary domain when reading a MED file)
      // we might end up calling the current method with a sequential MD_Vector:
      if (sub_type(MD_Vector_seq, md.valeur())) return;
      //Only time synchronous virtual exchange
      if (is_exchange_blocking == IsExchangeBlocking::DefaultBlocking) statistics().begin_count(STD_COUNTERS::virtual_swap);
      select_virtual_exchange_operation(v.get_md_vector(), v, opt,is_exchange_blocking, kernel_name);
      if (is_exchange_blocking == IsExchangeBlocking::DefaultBlocking) statistics().end_count(STD_COUNTERS::virtual_swap);
    }
  //else Cerr << "Warning: A call to ::virtual_exchange() is done on a non-distributed vector." << finl; /Process::exit();
}
 
void MD_Vector_tools::echange_espace_virtuel(IntVect& v, Operations_echange opt, IsExchangeBlocking is_exchange_blocking, const std::string kernel_name) { call_virtual_exchange<int>(v,opt, is_exchange_blocking, kernel_name); }
#if INT_is_64_ == 2
void MD_Vector_tools::echange_espace_virtuel(TIDVect& v, Operations_echange opt, IsExchangeBlocking is_exchange_blocking, const std::string kernel_name) { call_virtual_exchange<trustIdType>(v,opt,is_exchange_blocking, kernel_name); }
#endif
void MD_Vector_tools::echange_espace_virtuel(DoubleVect& v, Operations_echange opt, IsExchangeBlocking is_exchange_blocking, const std::string kernel_name) { call_virtual_exchange<double>(v,opt,is_exchange_blocking, kernel_name); }
void MD_Vector_tools::echange_espace_virtuel(FloatVect& v, Operations_echange opt, IsExchangeBlocking is_exchange_blocking, const std::string kernel_name) { call_virtual_exchange<float>(v,opt,is_exchange_blocking, kernel_name); }
 
void MD_Vector_tools::compute_sequential_items_index(const MD_Vector&, MD_Vector_renumber&, int line_size)
{
  Cerr << "MD_Vector_tools::compute_sequential_items_index" << finl;
  Process::exit();
}
 
/*! @brief Creates a descriptor for a subset of a vector.
 *
 * renum provides the structure and the descriptor of the source vector (must have line_size==1).
 *   renum must contain -1 for items of the source vector to discard and a
 *    non-negative value for items to keep.
 *   The value renum[i] gives the index of that item in the new array.
 *   Note: if an item of the source vector is received from another processor and must be kept,
 *    that same item must also be kept on the processor that sends it.
 *
 */
void MD_Vector_tools::creer_md_vect_renum(const IntVect& renum, MD_Vector& md_vect)
{
  const MD_Vector& src_md = renum.get_md_vector();
  if (!src_md)
    {
      Cerr << "Internal error in MD_Vector_tools::creer_md_vect_renum: null descriptor!" << finl;
      Process::exit();
    }
  if (renum.line_size() != 1)
    {
      Cerr << "Internal error in MD_Vector_tools::creer_md_vect_renum: line_size != 1" << finl;
      Process::exit();
    }
  src_md->fill_md_vect_renum(renum, md_vect);
}
 
/*! @brief Same as creer_md_vect_renum() but creates a default numbering.
 *
 * The flags_renum array must contain on input a NON-NEGATIVE value for the
 *   values to keep and a negative value for the others.
 *   On output, flags_renum contains the index of the item in the reduced array if it is
 *    kept, otherwise the original unmodified value.
 *   Items are placed in ascending order of their local index in the original descriptor.
 *
 */
void MD_Vector_tools::creer_md_vect_renum_auto(IntVect& flags_renum, MD_Vector& md_vect)
{
  const int n = flags_renum.size_array();
  int count = 0;
  for (int i = 0; i < n; i++)
    {
      int x = flags_renum[i];
      if (x >= 0)
        flags_renum[i] = count++;
    }
  creer_md_vect_renum(flags_renum, md_vect);
}
 
/*! @brief dumps vector v with its parallel descriptor to os.
 *
 * os must be a file type with one file per process (no shared file).
 *
 */
void MD_Vector_tools::dump_vector_with_md(const DoubleVect& v, Sortie& os)
{
  assert(v.line_size() == 1); // Only Vect! not tab.
 
  const MD_Vector_base& md = v.get_md_vector().valeur();
  os << md.que_suis_je() << finl;
  os << md << finl;
 
  os << v.size_array() << finl;
  os.put(v.addr(), v.size_array(), v.line_size());
  os << finl;
}
 
/*! @brief restores a vector saved by dump_vector_with_md.
 *
 * The code must run with the same number of processors.
 *   "v" will have a newly created copy of the descriptor (not equal to any
 *   other descriptor previously loaded). If "identical" descriptors are needed
 *   for several vectors, you might want to call set_md_vector() afterwards to
 *   attach another descriptor. See save_matrice_conditionnel()
 */
void MD_Vector_tools::restore_vector_with_md(DoubleVect& v, Entree& is)
{
 
  v.reset();
  OWN_PTR(MD_Vector_base) md_ptr;
  Nom md_type;
  is >> md_type;
  md_ptr.typer(md_type);
  is >> md_ptr.valeur();
 
  // Create the MD_Vector attached to the array
  MD_Vector md;
  md.copy(md_ptr.valeur());
 
  DoubleVect toto;
  toto.resize(md_ptr->get_nb_items_tot(), RESIZE_OPTIONS::NOCOPY_NOINIT);
  int size_tot;
  is >> size_tot;
  is.get(toto.addr(), size_tot);
  toto.set_md_vector(md);
 
  // Attach v to toto and forget toto...
  v.reset();
  v.ref(toto);
}