DistributedVector.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"
 
#ifdef CEPS_USE_PETSC
  // Use the wrapper around PETSc Vec
  #include "linearAlgebra/PETScVector.cpp"
#endif  // CEPS_USE_PETSC
 
CepsString
DistributedVector::getUnit() const
{
  return m_unit;
}
 
void
DistributedVector::setUnit(const CepsString& u)
{
  m_unit = u;
}
 
CepsBool
DistributedVector::isAssembled() const
{
  return m_isAssembled;
}
 
void
DistributedVector::checkNanOrInf(CepsString message)
{
  this->getLocalData();
  for (CepsGlobalIndex i=m_lo;i<m_hi;i++)
    ceps::checkNanOrInf(this->operator[](i),message+" coefficient #"+ceps::toString(i)+" ");
  this->releaseLocalData();
}
 
// -------------------------------------------------------------------------====
// Subvectors operations
// -------------------------------------------------------------------------====
 
void
DistributedVector::addSubVector(
  const CepsMathDynamic1D           &subVec,
  const CepsVector<CepsGlobalIndex> &rowIndices
)
{
  addValues(subVec.data(), rowIndices.size(), rowIndices.data());
  return;
}
 
void
DistributedVector::insertSubVector(
  const CepsMathDynamic1D           &subVec,
  const CepsVector<CepsGlobalIndex> &rowIndices
)
{
  setValues(subVec.data(), rowIndices.size(), rowIndices.data());
  return;
}
 
// -------------------------------------------------------------------------====
// OPERATORS
// -------------------------------------------------------------------------====
 
DistributedVector &
DistributedVector::operator=(const DistributedVector &v)
{
  v.duplicate(*this,true);
  return *this;
}
 
// Access: non const version
inline CepsReal &
DistributedVector::operator[] (CepsIndex globalIndex)
{
  CEPS_ABORT_IF(not m_mayReadData,
    "trying to access to memory without required called to getLocalData()."
  );
 
  CEPS_ABORT_IF((globalIndex < m_lo) or (globalIndex >= m_hi),
    "out-of-bounds index " << globalIndex << " while range of vector is ( low-high " << m_lo << ","
                           << m_hi << ")"
  );
 
  return m_localData[globalIndex-m_lo];
}
 
// Access: const version
inline CepsReal
DistributedVector::operator[] (CepsGlobalIndex globalIndex) const
{
  //CEPS_ABORT_IF (
  //  not m_mayReadData, "trying to access to memory without required called to getLocalData()."
  //);
  // get/realeaseLocalData is not compatible with const...
 
  CEPS_ABORT_IF((globalIndex < m_lo) or (globalIndex >= m_hi),
    "out-of-bounds index " << globalIndex << " while range of vector is ( low-high " << m_lo << ","
                           << m_hi << ")"
  );
 
  return m_localData[globalIndex-m_lo];
}
 
DistributedVector
operator+ (const DistributedVector &lhs, const DistributedVector &rhs)
{
  return DistributedVector(lhs) += rhs;
}
 
DistributedVector
operator+ (const DistributedVector &lhs, const CepsReal &scalar)
{
  return DistributedVector(lhs) += scalar;
}
 
DistributedVector
operator-(const DistributedVector &lhs, const DistributedVector &rhs)
{
  return DistributedVector(lhs) -= rhs;
}
 
DistributedVector
operator-(const DistributedVector &lhs, const CepsReal &scalar)
{
  return DistributedVector(lhs) -= scalar;
}
 
DistributedVector
operator*(const DistributedVector &lhs, const CepsReal &scalar)
{
  return DistributedVector(lhs) *= scalar;
}
 
DistributedVector
operator*(const CepsReal &scalar, const DistributedVector &rhs)
{
  return DistributedVector(rhs) *= scalar;
}
 
DistributedVector
operator*(const DistributedMatrix &A, const DistributedVector &x)
{
  DistributedVector res(x);
  res.mult(A,x);
  return res;
}
 
void
DistributedVector::fillFromVector(const CepsVector<CepsReal> &src)
{
  // Gather local sizes of vectors
  CepsUInt M     = ceps::getGridSize();
  int     *sizes = nullptr;
 
  if (ceps::isMaster())
    sizes = new int[M];
  MPI_Gather(
    &m_localSize,1,CEPS_MPI_INT,sizes,1,MPI_INT,CEPS_MASTER_PROC,ceps::getCommunicator()
  );
 
  // Create offset array
  int *displ = nullptr;
  if (ceps::isMaster())
  {
    displ    = new int[M];
    displ[0] = 0;
    for (CepsUInt i = 1; i < M; ++i)
      displ[i] = displ[i - 1] + sizes[i - 1];
  }
 
  // Send values directly into the vector
  getLocalData();
  MPI_Scatterv(
    src.data(), sizes, displ, CEPS_MPI_REAL, m_localData, m_localSize, CEPS_MPI_REAL,
    CEPS_MASTER_PROC, ceps::getCommunicator()
  );
  releaseLocalData();
 
  if (ceps::isMaster())
  {
    ceps::destroyTabular(sizes);
    ceps::destroyTabular(displ);
  }
  return;
}