PETScVector.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/DistributedVector.hpp"
 
DistributedVector::DistributedVector():
  m_unit        ("None" ),
  m_lo          (0      ),
  m_hi          (0      ), 
  m_globalSize  (0      ),
  m_localSize   (0      ),
  m_localData   (nullptr),
  m_mayReadData (false  ),
  m_isAssembled (false  ),
  #ifdef PROFILING_VECTOR
  m_totalAssembly(0),
  m_tickAssembly (0),
  #endif
  m_initialized (false  ),
  m_vecCreated  (false  )
{
}
 
DistributedVector::DistributedVector(const DistributedVector &v):
  DistributedVector()
{
  createEmptyVector();
  *this = v;
}
 
DistributedVector::DistributedVector(CepsInt M, CepsInt m):
  DistributedVector()
{
  createEmptyVector();
  setSize(M,m);
}
 
DistributedVector::DistributedVector (
  std::shared_ptr<DistributedVector> templateVector,
  CepsBool                           copyValues,
  CepsBool                           isSequential
):
  DistributedVector()
{
  if (not isSequential)
  {
    createEmptyVector();
    templateVector->duplicate(*this,copyValues);
  }
  else
  {
    VecScatter context;
    // VecScatterCreateToZero(templateVector->getVector(), &context, &m_v);
    VecScatterCreateToAll(templateVector->getVector(),&context,&m_v);
    VecScatterBegin      (context,templateVector->getVector(),m_v,INSERT_VALUES,SCATTER_FORWARD);
    VecScatterEnd        (context,templateVector->getVector(),m_v,INSERT_VALUES,SCATTER_FORWARD);
    VecScatterDestroy    (&context);
 
    #ifdef PROFILING_VECTOR
      m_totalAssembly = m_tickAssembly = 0.0;
    #endif
 
    VecGetOwnershipRange(m_v,&m_lo,&m_hi);
    m_globalSize  = m_hi;
    m_localSize   = m_hi;
    m_initialized = true;
  }
}
 
DistributedVector::DistributedVector(
  DistributedVector& templateVector,
  CepsBool           copyValues,
  CepsBool           isSequential
):
  DistributedVector()
{
  if (not isSequential)
  {
    createEmptyVector();
    templateVector.duplicate(*this,copyValues);
  }
  else
  {
    VecScatter context;
    // VecScatterCreateToZero(templateVector->getVector(), &context, &m_v);
    VecScatterCreateToAll(templateVector.getVector(),&context,&m_v);
    VecScatterBegin      (context,templateVector.getVector(),m_v,INSERT_VALUES,SCATTER_FORWARD);
    VecScatterEnd        (context,templateVector.getVector(),m_v,INSERT_VALUES,SCATTER_FORWARD);
    VecScatterDestroy    (&context);
 
    #ifdef PROFILING_VECTOR
      m_totalAssembly = m_tickAssembly = 0.0;
    #endif
 
    VecGetOwnershipRange(m_v,&m_lo,&m_hi);
    m_globalSize  = m_hi;
    m_localSize   = m_hi;
    m_initialized = true;
  }
}
 
DistributedVector::DistributedVector(DMatPtr mat, CepsGlobalIndex row):
  DistributedVector()
{
  createEmptyVector();
  CepsInt         M,N;
  CepsGlobalIndex lo,hi;
  mat->getLocalRange(&lo,&hi);
  mat->getSize(&M,&N);
  setSize(M,hi-lo);
  CepsBool isOwner = lo <= row and row < hi;
 
  CepsInt            nbEntries;
  const CepsInt     *columnIndices;
  const PetscScalar *values;
 
  // as we are in constructor, we are sure that data will be exchanged
  if (isOwner)
  {
    MatGetRow(mat->getMatrix(),row,&nbEntries,&columnIndices,&values);
    VecSetValues(m_v,nbEntries,columnIndices,values,INSERT_VALUES);
  }
 
  // assemble
  finalize();
 
  if(isOwner)
    MatRestoreRow(mat->getMatrix(),row,&nbEntries,&columnIndices,&values);
 
}
 
DistributedVector::~DistributedVector()
{
  #ifdef PROFILING_VECTOR
    if (ceps::isMaster())
      std::cout << "/* PROFILING: Distributed vector - time spent assembling = " << m_totalAssembly
                << " seconds. *\\" << std::endl;
  #endif
  if (m_vecCreated)
    VecDestroy(&(this->m_v)); 
}
 
const PetscVector&
DistributedVector::getVector() const
{
  return m_v;
}
 
CepsMathScalar *
DistributedVector::localData()
{
  CEPS_ABORT_IF(not m_mayReadData, "call getLocalData() before accessing local data.\n");
  return m_localData;
}
 
void
DistributedVector::setSize(CepsInt M, CepsInt m)
{
  CEPS_ABORT_IF(M < 0 or m < 0, "(M = " << M << ", n = " << m << "): cannot set negative size!");
  VecSetSizes(m_v,m,M);
  m_globalSize = M;
  m_localSize  = m;
  VecGetOwnershipRange(m_v,&m_lo,&m_hi);
  m_initialized = true;
}
 
void
DistributedVector::getSize(CepsInt *globalSize) const
{
  CepsInt s;
  VecGetSize(m_v, &s);
  *globalSize = s;
  return;
}
 
CepsInt
DistributedVector::getSize() const
{
  CepsInt s;
  VecGetSize(m_v, &s);
  return s;
}
 
void
DistributedVector::getLocalSize(CepsInt *localSize) const
{
  CepsInt s;
  VecGetLocalSize(m_v, &s);
  *localSize = s;
  return;
}
 
CepsInt
DistributedVector::getLocalSize() const
{
  CepsInt s;
  VecGetLocalSize(m_v,&s);
  return s;
}
 
void
DistributedVector::getLocalRange(CepsGlobalIndex* lo, CepsGlobalIndex* hi) const
{
  *lo = m_lo;
  *hi = m_hi;
  return;
}
 
void
DistributedVector::getLocalData()
{
  VecGetArray(m_v,&m_localData);
  m_mayReadData = true;
  return;
}
 
void
DistributedVector::releaseLocalData()
{
  m_mayReadData = false;
  VecRestoreArray(m_v,&m_localData);
  return;
}
 
void
DistributedVector::setValue(CepsMathScalar value, CepsGlobalIndex row)
{
  VecSetValue(m_v,row,value,INSERT_VALUES);
  return;
}
 
void
DistributedVector::addValue(CepsMathScalar value, CepsGlobalIndex row)
{
  VecSetValue(m_v,row,value,ADD_VALUES);
  return;
}
 
void
DistributedVector::ignoreOffProcEntries(CepsBool flag)
{
  PetscBool value = flag ? PETSC_TRUE : PETSC_FALSE;
  VecSetOption(m_v,VEC_IGNORE_OFF_PROC_ENTRIES,value);
  return;
}
 
void
DistributedVector::setValues(
  const CepsMathScalar  *values,
  CepsInt                n,
  const CepsGlobalIndex *indices
)
{
  VecSetValues(m_v,n,indices,values,INSERT_VALUES);
  return;
}
 
void
DistributedVector::setValuesLocal(
  DistributedVector     &values,
  CepsInt                n,
  const CepsGlobalIndex *indices
)
{
  // request local access on both vectors
  getLocalData();
  values.getLocalData();
 
  if (ceps::isValidPtr(indices))
    for (CepsInt i=0; i<n; i++)
      m_localData[indices[i]-m_lo] = values[indices[i]];
  else
    // contiguous, starting from first owned row
    for (CepsInt i=0; i<n; i++)
      m_localData[i] = values.m_localData[i];
 
  // done with local data
  releaseLocalData();
  values.releaseLocalData();
  return;
}
 
void
DistributedVector::setValuesLocal(
  const CepsMathScalar  *values,
  CepsInt                n,
  const CepsGlobalIndex *indices
)
{
  // request access to local data
  getLocalData();
 
  if (ceps::isValidPtr(indices))
    for (CepsInt i=0; i<n; i++)
      m_localData[indices[i]-m_lo] = values[i];
  else
    // contiguous, starting from first owned row
    for (CepsInt i=0; i<n; i++)
      m_localData[i] = values[i];
  // done with local data
  releaseLocalData();
  return;
}
 
void
DistributedVector::setValuesLocal(
  const CepsVector<CepsMathScalar> &values,
  const CepsVector<CepsInt>        &indices
)
{
  CepsUInt N = values.size();
  CEPS_ABORT_IF(indices.size() != N,
    "Given arrays of values and row indices have different size."
  );
  getLocalData();
  for (CepsUInt i=0U; i<N; i++)
    m_localData[indices[i]-m_lo] = values[i];
  releaseLocalData();
  return;
}
 
void
DistributedVector::addValues(
  const CepsMathScalar  *values,
  CepsInt                n,
  const CepsGlobalIndex *indices
)
{
  VecSetValues(m_v,n,indices,values,ADD_VALUES);
  return;
}
 
void
DistributedVector::addValuesLocal(
  DistributedVector     &values,
  CepsInt                n,
  const CepsGlobalIndex *indices
)
{
  // request local access on both vectors
  getLocalData();
  values.getLocalData();
 
  if (ceps::isValidPtr(indices))
    for (CepsInt i=0; i<n; i++)
      m_localData[indices[i]-m_lo] += values[indices[i]];
  else
  {
    // contiguous, starting from first owned row
    for (CepsInt i=0; i<n; i++) {
      m_localData[i] += values.m_localData[i];
    }
  }
 
  // done with local data
  releaseLocalData();
  values.releaseLocalData();
  return;
}
 
void
DistributedVector::addValuesLocal(
  const CepsMathScalar  *values,
  CepsInt                n,
  const CepsGlobalIndex *indices
)
{
  // request access to local data
  getLocalData();
  CepsInt i;
 
  if (ceps::isValidPtr(indices))
    for (i=0; i<n; i++)
      m_localData[indices[i]-m_lo] += values[i];
  else
    // contiguous, starting from first owned row
    for (i=0; i<n; i++)
      m_localData[i] += values[i];
  // done with local data
  releaseLocalData();
  return;
}
 
void
DistributedVector::getValues(CepsMathScalar *values, CepsInt n, const CepsInt *indices)
{
  CEPS_ABORT_IF(not m_mayReadData,
    "invalid access to local data.\nMake sure you called "
    "getLocalData() beforehand."
  );
 
  // standard non-ghost version
  for (CepsInt i=0; i<n; i++)
    values[i] = m_localData[indices[i]-m_lo];
  return;
}
 
CepsVector<CepsMathScalar>
DistributedVector::getSubVector(const CepsVector<CepsInt> &indices)
{
  CepsVector<CepsMathScalar> res;
  getLocalData();
  // extract from local data
  for (CepsInt row : indices)
  {
    CEPS_ABORT_IF(row<m_lo or row>=m_hi,
      "trying to access row " << row << "which is out of local range ("
                              << m_lo << "," << m_hi << ")"
    );
    res.push_back(m_localData[row-m_lo]);
  }
 
  // release array
  releaseLocalData();
  return res;
}
 
void
DistributedVector::zero()
{
  VecZeroEntries(m_v);
  return;
}
 
void
DistributedVector::fill(CepsMathScalar value)
{
  VecSet(m_v,value);
  return;
}
 
void
DistributedVector::duplicate(DistributedVector& v, CepsBool copyValues) const
{
  if (not v.m_vecCreated)
    v.createEmptyVector();
 
  if (copyValues)
  {
    v.setSize(m_globalSize,m_localSize);
    VecCopy(m_v,v.m_v);
  }
  else
  {
    VecDuplicate(m_v,&v.m_v);
    v.setSize(m_globalSize,m_localSize);
  }
  return;
}
 
void
DistributedVector::beginAssembly()
{
#ifdef PROFILING_VECTOR
  m_tickAssembly = ceps::tick();
#endif
  m_isAssembled = false;
  VecAssemblyBegin(m_v);
  return;
}
 
void
DistributedVector::endAssembly()
{
  VecAssemblyEnd(m_v);
  m_isAssembled = true;
#ifdef PROFILING_VECTOR
  m_totalAssembly += ceps::tock(m_tickAssembly);
#endif
  return;
}
 
void
DistributedVector::finalize()
{
#ifdef PROFILING_VECTOR
  m_tickAssembly = ceps::tick();
#endif
  beginAssembly();
  endAssembly();
#ifdef PROFILING_VECTOR
  m_totalAssembly += ceps::tock(m_tickAssembly);
#endif
  return;
}
 
CepsBool
DistributedVector::sameSize(const DistributedVector &vec) const
{
  return (m_globalSize == vec.m_globalSize);
}
 
CepsBool
DistributedVector::sameLocalRange(const DistributedVector &vec) const
{
  if (! this->sameSize(vec))
    return false;
 
  return ((m_lo == vec.m_lo) and (m_hi == vec.m_hi));
}
 
CepsBool
DistributedVector::equals(const DistributedVector &vec, CepsMathScalar errorFactor) const
{
  CepsInt lo, hi;
  vec.getLocalRange(&lo, &hi);
 
  // check size
  CEPS_ABORT_IF(not sameSize(vec),
    "DistributedVector::equals: trying to compare vectors of sizes "
    << m_globalSize << " and " << vec.m_globalSize
  );
  // check ownership range
  CEPS_ABORT_IF(not sameLocalRange(vec),
    "DistributedVector:equals: trying to compare two vectors with different data repartition."
  );
 
  // check the values
  const_cast<DistributedVector *> (this)->getLocalData();
  const_cast<DistributedVector &> (vec) . getLocalData();
 
  CepsUInt myflag = 1U, flag = 0U;
  for (CepsInt i=0; i<(hi-lo); i++)
  {
    if (! ceps::equals(vec.m_localData[i], m_localData[i], errorFactor))
    {
      myflag = 0U;
      break;
    }
  }
 
  // Share results
  MPI_Allreduce(&myflag,&flag,1,CEPS_MPI_UINT,MPI_MIN,ceps::getCommunicator());
 
  const_cast<DistributedVector *> (this)->releaseLocalData();
  const_cast<DistributedVector &> (vec)  .releaseLocalData();
 
  return (CepsBool) flag;
}
 
CepsBool
DistributedVector::approxEquals(const DistributedVector &vec, CepsMathScalar epsilon) const
{
  CepsInt lo, hi;
  vec.getLocalRange (&lo, &hi);
 
  // check size
  CEPS_ABORT_IF(not sameSize(vec),
    "trying to compare vectors of sizes " << m_globalSize << " and " << vec.m_globalSize
  );
  // check ownership range
  CEPS_ABORT_IF(not sameLocalRange(vec),
    "trying to compare two vectors with different data repartition."
  );
 
  // check the values
  const_cast<DistributedVector *> (this)->getLocalData();
  const_cast<DistributedVector &> (vec)  .getLocalData();
 
  CepsUInt myflag = 1U, flag = 0U;
  for (CepsInt i=0; i<(hi-lo); i++)
  {
    if (not ceps::approxEquals(vec.m_localData[i],m_localData[i],epsilon))
    {
      myflag = false;
      //      break;
    }
  }
  // Share results
  MPI_Allreduce(&myflag, &flag, 1, CEPS_MPI_UINT, MPI_MIN, ceps::getCommunicator());
 
  const_cast<DistributedVector *> (this)->releaseLocalData();
  const_cast<DistributedVector &> (vec)  .releaseLocalData();
 
  return (CepsBool) flag;
}
 
void
DistributedVector::add(const DistributedVector &x)
{
  this->addScaled(x,1.0);
  return;
}
 
void
DistributedVector::addScaled(const DistributedVector &x, CepsMathScalar alpha)
{
  CEPS_ABORT_IF(ceps::isSamePtr(this,&x),
    "cannot compute x = x + a.x when using PETSc."
  );
 
  // check sizes and allocate if not initialized
  CEPS_ABORT_IF(not x.m_initialized,
    "added vector is not properly initialized."
  );
  checkAndSetSizes(x.m_globalSize,x.m_localSize);
  VecAXPY(m_v,alpha,x.getVector());
  return;
}
 
void
DistributedVector::scale(CepsReal alpha)
{
  VecScale(m_v,alpha);
}
 
void 
DistributedVector::pointWiseScale(const DistributedVector& v, CepsBool divide)
{
  CEPS_ABORT_IF(not sameLocalRange(v),
    "Trying to multiply vectors with different local ranges"
  );
  const CepsMathScalar *src;
 
  getLocalData();
  VecGetArrayRead(v.m_v,&src);
  if(divide)
  {
    for (CepsInt i=0; i<m_hi-m_lo; ++i)
      m_localData[i] /= src[i];
  }
  else 
  {
    for (CepsInt i=0; i<m_hi-m_lo; ++i)
          m_localData[i] *= src[i];
  }
  VecRestoreArrayRead(v.m_v,&src);
  releaseLocalData();
 
  return;
}
 
void
DistributedVector::mult(const DistributedMatrix& A, const DistributedVector& x)
{
  CEPS_ABORT_IF (ceps::isSamePtr(this, &x),
    "cannot compute x = A.x when using PETSc."
  );
 
  // check sizes and allocate if not initialized
  CEPS_ABORT_IF(not x.m_initialized,
    "x vector is not properly initialized."
  );
  checkAndSetSizes(x.m_globalSize,x.m_localSize);
 
  CepsInt matrixMSize, matrixNSize;
  A.getSize(&matrixMSize, &matrixNSize);
  CEPS_ABORT_IF(x.m_globalSize != matrixMSize, 
    "cannot compute A.x as A has size (" << matrixMSize << "," << matrixNSize
    << ") and vector has size " << x.m_globalSize
  );
 
  // this is not necessary, but keep it here for performance issues
  CepsGlobalIndex matrixLo, matrixHi;
  A.getLocalRange(&matrixLo, &matrixHi);
  CEPS_ABORT_IF((m_lo != matrixLo) or (m_hi != matrixHi), 
    "vector and matrix have different local ranges ("
    << m_lo << "," << m_hi << ") and (" << matrixLo
    << "," << matrixHi << ") respectively."
  );
 
  MatMult(const_cast<PetscMatrix>(A.getMatrix()),x.m_v,m_v);
  m_isAssembled = true;
  return;
}
 
void
DistributedVector::abs()
{
  VecAbs(m_v);
  return;
}
 
CepsMathScalar
DistributedVector::dot(DistributedVector &x)
{
  CepsMathScalar dot_val;
  VecDot (m_v,x.getVector(), &dot_val);
  return dot_val;
}
 
CepsMathScalar
DistributedVector::lInfNorm() const
{
  CepsMathScalar norm;
  VecNorm(m_v,NORM_INFINITY,&norm);
  return norm;
}
 
CepsMathScalar
DistributedVector::l1Norm() const
{
  CepsMathScalar norm;
  VecNorm (m_v,NORM_1,&norm);
  return norm;
}
 
CepsMathScalar
DistributedVector::l2Norm() const
{
  CepsMathScalar norm;
  VecNorm (m_v,NORM_2,&norm);
  return norm;
}
 
void
DistributedVector::view() const
{
  VecView(m_v,PETSC_VIEWER_STDOUT_WORLD);
  return;
}
 
// Operators
 
DistributedVector &
DistributedVector::operator+=(const DistributedVector &v)
{
  CEPS_ABORT_IF(not sameLocalRange(v),
    "Trying to add vectors with different local ranges"
  );
  const CepsMathScalar *src;
 
  getLocalData ();
  VecGetArrayRead (v.m_v,&src);
  for (CepsInt i=0; i< m_hi-m_lo; ++i)
    m_localData[i] += src[i];
  VecRestoreArrayRead(v.m_v,&src);
  releaseLocalData();
 
  return *this;
}
 
DistributedVector &
DistributedVector::operator+=(const CepsMathScalar &scalar)
{
  VecShift(m_v,scalar);
  return *this;
}
 
DistributedVector &
DistributedVector::operator-=(const DistributedVector &v)
{
  CEPS_ABORT_IF(not sameLocalRange(v),
    "Trying to add vectors with different local ranges"
  );
  const CepsMathScalar *src;
 
  getLocalData();
  VecGetArrayRead(v.m_v,&src);
  for (CepsInt i=0; i<m_hi-m_lo; ++i)
    m_localData[i] -= src[i];
 
  VecRestoreArrayRead(v.m_v,&src);
  releaseLocalData();
 
  return *this;
}
 
DistributedVector &
DistributedVector::operator-=(const CepsMathScalar &scalar)
{
  return this->operator+=(-scalar);
}
 
DistributedVector &
DistributedVector::operator*=(const CepsMathScalar &scalar)
{
  VecScale(m_v,scalar);
  return *this;
}
 
// ---------------------------------------------------------------------------------------------------------==
// Private, protected
// ---------------------------------------------------------------------------------------------------------==
 
void
DistributedVector::checkAndSetSizes(CepsInt gSize, CepsInt lSize)
{
  if (m_initialized)
  {
    CEPS_ABORT_IF (gSize != m_globalSize or lSize != m_localSize,
      "DistributedVector::checkAndSetSizes(): calling vector is not of proper sizes."
    );
  }
  else
  {
    createEmptyVector();
    setSize(gSize,lSize);
  }
  return;
}
 
void
DistributedVector::createEmptyVector()
{
  if (not m_vecCreated)
  {
    VecCreate(ceps::getCommunicator(),&(this->m_v));
    VecSetType(m_v,VECSTANDARD);
    m_vecCreated = true;
  }
  return;
}