DistributedHaloVector.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   CEPS is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "linearAlgebra/DistributedHaloVector.hpp"
 
#ifdef CEPS_USE_PETSC
 
DistributedHaloVector::DistributedHaloVector()
{
  initialize();
}
 
DistributedHaloVector::DistributedHaloVector(DistributedHaloVector& v)
{
  initialize();
  // call assignment operator
  *this = v;
}
 
DistributedHaloVector::DistributedHaloVector(
  CepsIndex  M, 
  CepsIndex  m, 
  CepsIndex  nbHaloRows, 
  CepsIndex* haloRows
) :
  m_nbHaloRows(nbHaloRows)
{
  this->m_isAssembled = false;
  this->m_mayReadData = false;
 
  // create the vector with halo padding
  VecCreateGhost(ceps::getCommunicator(), m, M, nbHaloRows, haloRows, &(this->m_v));
  m_vecCreated = true;
 
  #ifdef PROFILING_VECTOR
    m_totalAssembly = m_tickAssembly = 0.0;
  #endif
 
  // update local copy of row ownership
  VecGetOwnershipRange(this->m_v, &(this->m_lo), &(this->m_hi));
  CepsIndex sum = 0;
  // mapping non-ghost values
  for (CepsIndex i = 0; i < m; i++, sum++)
    m_globalToLocalMapping[this->m_lo + i] = sum;
 
  // Mapping ghots values:
  // they are always located after the non-ghost values
  for (CepsIndex i = 0; i < nbHaloRows; i++, sum++)
    m_globalToLocalMapping[haloRows[i]] = sum;
 
  CEPS_ABORT_IF((this->m_hi - this->m_lo + nbHaloRows) != sum 
              or (CepsIndex) m_globalToLocalMapping.size () != sum,
    "DistributedHaloVector: fatal error, some DsOF were not mapped."
  );
 
  // Now ok to use this vector
  this->m_initialized = true;
  this->m_globalSize  = M;
  this->m_localSize   = m;
}
 
DistributedHaloVector::~DistributedHaloVector()
{
}
 
void
DistributedHaloVector::getLocalData()
{
  // Scatter the proper values before accessing the values.
  // Here, the scatter is forward, meaning the halo values
  // are set to the correct value from the owning processes:
  // i.e., any local modification in the halo rows is overwritten.
  VecGhostUpdateBegin(m_v,INSERT_VALUES,SCATTER_FORWARD);
  VecGhostUpdateEnd  (m_v,INSERT_VALUES,SCATTER_FORWARD);
 
  // Get the local representation of this vector
  VecGhostGetLocalForm(this->m_v,&m_localV);
  // Ok to get the local array
  VecGetArray(m_localV, &(this->m_localData));
 
  // ok to read data now
  this->m_mayReadData = true;
}
 
void
DistributedHaloVector::releaseLocalData()
{
  this->m_mayReadData = false;
  VecRestoreArray(m_localV, &(this->m_localData));
  VecGhostRestoreLocalForm(this->m_v, &m_localV);
}
 
void
DistributedHaloVector::getValues(CepsReal *values, CepsIndex n, const CepsIndex *indices)
{ 
 
  CEPS_ABORT_IF(not m_mayReadData,
    "Invalid access to local data.\n"
    "   Make sure you called DistributedVector::getLocalData() beforehand."
  );
 
  // Use the global to local mapping:
  // iterate on every global index parameter and check if
  // current process has it in local memory
  std::map<CepsIndex, CepsIndex>::iterator it;
  for (CepsIndex i = 0; i < n; i++)
  {
    it = m_globalToLocalMapping.find (indices[i]);
    // failure if current process does not own this row
    CEPS_ABORT_IF(it == m_globalToLocalMapping.end(),
      "Distributed vector: on process #" << ceps::getRank() << ", trying to access row #"
      << indices[i] << "\nwhich not local or halo."
    );
 
    // Now ok to read in local memory
    values[i] = this->m_localData[it->second];
  }
}
 
CepsInt
DistributedHaloVector::getNbHalo()
{
  return m_nbHaloRows;
}
 
void
DistributedHaloVector::initialize()
{
  m_mayReadData = false;
  m_initialized = false;
  m_isAssembled = false;
  m_globalSize = m_localSize = 0;
  #ifdef PROFILING_VECTOR
    m_totalAssembly = m_tickAssembly = 0.0;
  #endif
}
 
// OPERATORS
// non const version
CepsReal&
DistributedHaloVector::operator[](CepsIndex index)
{
  std::map<CepsIndex, CepsIndex>::iterator it = m_globalToLocalMapping.find(index);
 
  CEPS_ABORT_IF(not m_mayReadData,
    "access to memory not requested"
  );
  CEPS_ABORT_IF(it == m_globalToLocalMapping.end(),
    "out-of-bounds index " << index
  );
 
  return this->m_localData[it->second];
}
 
// const version
CepsReal
DistributedHaloVector::operator[](CepsIndex index) const
{
  std::map<CepsIndex, CepsIndex>::const_iterator it = m_globalToLocalMapping.find(index);
 
// NO : getReleaseLocalData not const compatible
//#ifdef CEPS_DEBUG_ENABLED
//  if (not m_mayReadData)
//    CEPS_ABORT ("DistributedHaloVector: access to memory not requested."); 
//#endif  // CEPS_DEBUG_ENABLED
 
  CEPS_ABORT_IF(it == m_globalToLocalMapping.end(),
    "out-of-bounds index " << index
  );
 
  return this->m_localData[it->second];
}
 
DistributedHaloVector&
DistributedHaloVector::operator=(DistributedHaloVector &v)
{
  CEPS_ABORT_IF(not v.m_initialized,
    "source vector is not initialized"
  );
 
  // duplicate
  if (not m_initialized)
    v.duplicate(*this,false);
  // copy
  CepsIndex size = v.m_hi-v.m_lo+v.m_nbHaloRows;
  v.getLocalData();
  getLocalData  ();
  CepsReal *source = v.localData();
  CepsReal *dest   = localData  ();
  memcpy(dest,source,size*sizeof(CepsReal));
  v.releaseLocalData();
  releaseLocalData  ();
 
  // copy the mapping
  m_globalToLocalMapping = v.m_globalToLocalMapping;
  // copy the number of halo rows
  m_nbHaloRows = v.m_nbHaloRows;
 
  // ok to use it now
  m_initialized = true;
 
  return *this;
}
 
DistributedHaloVector&
DistributedHaloVector::operator=(DistributedVector &v)
{
  CEPS_ABORT_IF(not m_initialized,
    "calling vector is not initialized"
  );
  CepsIndex hi,lo;
  v.getLocalRange(&lo,&hi);
  CEPS_ABORT_IF((lo != m_lo) or (hi != m_hi),
    "vectors are of different sizes"
  );
 
  v.getLocalData    ();
  this->getLocalData();
  CepsReal *source = v    .localData();
  CepsReal *dest   = this->localData();
  memcpy(dest,source,(hi-lo)*sizeof(CepsReal));
  v    .releaseLocalData();
  this->releaseLocalData();
 
  // ok to use it now
  m_initialized = true;
 
  return *this;
}
 
void
DistributedHaloVector::scale(CepsReal alpha)
{
  #ifdef CEPS_USE_PETSC
  VecScale(m_localV,alpha);
  #endif
  return;
}
 
void
DistributedHaloVector::abs()
{
  DistributedVector::abs();
  #ifdef CEPS_USE_PETSC
  VecAbs(m_localV);
  #endif
  return;
}
 
#endif  // CEPS_USE_PETSC