PETScMatrix.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/DistributedMatrix.hpp"
#include "linearAlgebra/DistributedVector.hpp"
 
DistributedMatrix::DistributedMatrix() :
    m_A          (nullptr),
    m_lo         (0U),
    m_hi         (0U),
    m_properlySet(false)
{
  // Users MUST call DistributedMatrix::setNonZeroStructure(...)
  // and setSize
  // after calling this constructor to actually create the matrix
  MatCreate(ceps::getCommunicator(), &m_A);
  MatSetType(m_A, MATAIJ);
  MatSetUp(m_A);
}
 
DistributedMatrix::DistributedMatrix(CepsInt M, CepsInt N, CepsInt m, CepsInt n) :
    m_A          (nullptr),
    m_lo         (0U),
    m_hi         (0U),
    m_properlySet(false)
{
  // Users MUST call DistributedMatrix::setNonZeroStructure(...)
  // after calling this constructor to actually create the matrix
  MatCreate(ceps::getCommunicator(),&m_A);
  setSize(M, N, m, n);
  MatSetType(m_A, MATAIJ);
  MatSetUp(m_A);
}
 
DistributedMatrix::DistributedMatrix(const DistributedMatrix &templateMatrix,CepsBool copyValues) :
    m_A          (nullptr),
    m_lo         (0U),
    m_hi         (0U),
    m_properlySet(true)
{
  templateMatrix.duplicate(*this, copyValues);
  MatGetOwnershipRange(m_A, &m_lo, &m_hi);
  MatSetUp(m_A);
}
 
DistributedMatrix::DistributedMatrix (const CepsString &file) :
    m_A          (nullptr),
    m_lo         (0U),
    m_hi         (0U),
    m_properlySet(false)
{
  MatCreate(ceps::getCommunicator(), &m_A);
  load(file);
  MatSetUp(m_A);
}
 
DistributedMatrix::~DistributedMatrix()
{
  if (ceps::isValidPtr(m_A))
    MatDestroy(&m_A);
}
 
void
DistributedMatrix::ignoreOffProcEntries(CepsBool flag)
{
  PetscBool pflag = flag ? PETSC_TRUE : PETSC_FALSE;
  MatSetOption (m_A, MAT_IGNORE_OFF_PROC_ENTRIES, pflag);
  return;
}
 
void
DistributedMatrix::setSymmetric(CepsBool flag)
{
  PetscBool pflag = flag ? PETSC_TRUE : PETSC_FALSE;
 
  MatSetOption(m_A, MAT_SYMMETRIC, pflag);
  MatSetOption(m_A, MAT_SYMMETRY_ETERNAL, pflag);
  return;
}
 
void
DistributedMatrix::setSize(CepsInt M, CepsInt N, CepsInt m, CepsInt n)
{
  CEPS_ABORT_IF( not ((M >= 0 && N >= 0) or (m >= 0 && n >= 0)),
   "Invalid matrix sizes provided : got " << CEPS_DISPLAY_VAL (M, N, m, n)
  );
  MatSetSizes(m_A, m, n, M, N);
  return;
}
 
void
DistributedMatrix::setNonZeroStructure(const CepsInt d_nnz[], const CepsInt o_nnz[])
{
  if (ceps::isParallel())
  {
    MatMPIAIJSetPreallocation(m_A, 0, d_nnz, 0, o_nnz);
  }
  else
    MatSeqAIJSetPreallocation(m_A, 0, d_nnz);
  MatGetOwnershipRange(m_A, &m_lo, &m_hi);
  // now ok to use this matrix
  m_properlySet = true;
  return;
}
 
void
DistributedMatrix::setNonZeroStructure(CepsInt d_nz, CepsInt o_nz)
{
  if (ceps::isParallel())
    MatMPIAIJSetPreallocation(m_A, d_nz, nullptr, o_nz, nullptr);
  else
    MatSeqAIJSetPreallocation(m_A, d_nz, nullptr);
  MatGetOwnershipRange(m_A, &m_lo, &m_hi);
  // now ok to use this matrix
  m_properlySet = true;
  return;
}
 
void
DistributedMatrix::setValue(CepsMathScalar value, CepsGlobalIndex row, CepsGlobalIndex column)
{
  MatSetValue(m_A, row, column, value, INSERT_VALUES);
  return;
}
 
void
DistributedMatrix::addValue(CepsMathScalar value, CepsGlobalIndex row, CepsGlobalIndex column)
{
  MatSetValue(m_A, row, column, value, ADD_VALUES);
  return;
}
 
void
DistributedMatrix::getValues(
  CepsMathScalar        *values,
  CepsInt                nbRows,
  CepsInt                nbColumns,
  const CepsGlobalIndex *rowIndices,
  const CepsGlobalIndex *columnIndices
)
{
  MatGetValues(m_A, nbRows, rowIndices, nbColumns, columnIndices, values);
  return;
}
 
void
DistributedMatrix::setValues(
  const CepsMathScalar  *values,
  CepsInt                nbRows,
  CepsInt                nbColumns,
  const CepsGlobalIndex *rowIndices,
  const CepsGlobalIndex *columnIndices
)
{
  MatSetValues(m_A, nbRows, rowIndices, nbColumns, columnIndices, values, INSERT_VALUES);
  return;
}
 
void
DistributedMatrix::addValues(
  const CepsMathScalar  *values,
  CepsInt                nbRows,
  CepsInt                nbColumns,
  const CepsGlobalIndex *rowIndices,
  const CepsGlobalIndex *columnIndices
)
{
  MatSetValues(m_A, nbRows, rowIndices, nbColumns, columnIndices, values, ADD_VALUES);
  return;
}
 
void
DistributedMatrix::zeroRowsAndColumns(
  CepsInt          nbRows,
  CepsGlobalIndex *rows,
  CepsMathScalar   diagonalValue
)
{
  MatSetOption(m_A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
  MatZeroRowsColumns(m_A, nbRows, rows, diagonalValue, nullptr, nullptr);
  return;
}
 
void
DistributedMatrix::zeroRows(CepsInt nbRows, CepsGlobalIndex *rows, CepsMathScalar diagonalValue)
{
  MatSetOption(m_A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
  MatSetOption(m_A, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_FALSE);
  MatZeroRows (m_A, nbRows, rows, diagonalValue, PETSC_NULL, PETSC_NULL);
  return;
}
 
void
DistributedMatrix::zero()
{
  MatZeroEntries(m_A);
  return;
}
 
void
DistributedMatrix::identity()
{
  CepsInt mGlobal, nGlobal;
  MatGetSize(m_A, &mGlobal, &nGlobal);
  CEPS_ABORT_IF(mGlobal != nGlobal,
    "cannot set non-square matrix to identity.\nCurrent size is ("
    << mGlobal << "," << nGlobal << ")"
  );
 
  // proceed by setting non-zero structure: 1 value per row in
  // diagonal, 0 elsewhere
  this->setNonZeroStructure(1,0);
  CepsGlobalIndex i;
  for (i=m_lo; i<m_hi; i++)
    setValue(1.0,i,i);
  return;
}
 
void
DistributedMatrix::add(DistributedMatrix &X, CepsMathScalar a, CepsBool sameNonZeroStructure)
{
  if (sameNonZeroStructure)
    MatAXPY(m_A, a, X.m_A, SAME_NONZERO_PATTERN);
  else
    MatAXPY(m_A, a, X.m_A, DIFFERENT_NONZERO_PATTERN);
  return;
}
 
CepsMathScalar
DistributedMatrix::lInfNorm() const
{
  assertAssembled("computing L infty norm");
  CepsReal res = 0.0;
  MatNorm(m_A, NORM_INFINITY, &res);
  return res;
}
 
CepsMathScalar
DistributedMatrix::l1Norm() const
{
  assertAssembled("computing L 1 norm");
  CepsReal res = 0.0;
  MatNorm(m_A, NORM_1, &res);
  return res;
}
 
CepsMathScalar
DistributedMatrix::l2Norm() const
{
  assertAssembled("computing L 2 norm");
  CepsReal res = 0.0;
  MatNorm(m_A, NORM_FROBENIUS, &res);
  return res;
}
 
void
DistributedMatrix::copy(DistributedMatrix &A) const
{
  // Not to be used on a non-finalized matrix. Call
  // DistributedMatrix::finalize() before copying. <br>
  assertAssembled("calling copy");
  MatCopy (m_A, A.m_A, DIFFERENT_NONZERO_PATTERN);
  return;
}
 
void
DistributedMatrix::duplicate(DistributedMatrix &A, CepsBool flag) const
{
  // Not to be used on a non-finalized matrix. Call
  // DistributedMatrix::finalize() before copying. <br>
  if (flag)
    MatDuplicate(m_A, MAT_COPY_VALUES, &(A.m_A));
  else
    MatDuplicate(m_A, MAT_DO_NOT_COPY_VALUES, &(A.m_A));
  A.m_lo          = m_lo;
  A.m_hi          = m_hi;
  A.m_properlySet = true;
  return;
}
 
void
DistributedMatrix::getSize(CepsInt *M, CepsInt *N) const
{
  checkProperlySet ();
  MatGetSize (m_A, M, N);
  return;
}
 
void
DistributedMatrix::getLocalSize(CepsInt *m, CepsInt *n) const
{
  checkProperlySet ();
  MatGetLocalSize (m_A, m, n);
  return;
}
 
void
DistributedMatrix::getLocalRange(CepsGlobalIndex *lo, CepsGlobalIndex *hi) const
{
  checkProperlySet ();
  *lo = m_lo;
  *hi = m_hi;
}
 
void
DistributedMatrix::getLocalRow(
  const CepsMathScalar  **values,
  CepsGlobalIndex         row,
  CepsInt                *nbNonZero,
  const CepsGlobalIndex **cols
) const
{
  CEPS_ABORT_IF (row < m_lo or row >= m_hi,
    "Trying to access row " << row << "on process that owns range " << m_lo << "to " << m_hi
  );
  MatGetRow (m_A, row, nbNonZero, cols, values);
  return;
}
 
void
DistributedMatrix::releaseLocalRow(
  const CepsMathScalar  **values,
  CepsGlobalIndex         row,
  CepsInt                *nbNonZero,
  const CepsGlobalIndex **cols
)
{
  MatRestoreRow (m_A, row, nbNonZero, cols, values);
  return;
}
 
void
DistributedMatrix::getLocalRow (CepsVector<CepsMathScalar> &rowValues, CepsGlobalIndex row)
{
  /*
   * PETSc implementation:
   * Due to the way PETSc is designed, we cannot directly return
   * a pointer to the buffer holding the local values in the matrix;
   * that would require 2 function calls: one to get the pointer
   * and another to release it to ensure that the "application does not
   * bleed memory"(c.f PETSc doc).
   * What we do here is we <B>copy</B> the non-zero values from the local buffer
   * into a vector.
   * Obviously, the values in the matrix are not altered.
   */
 
  checkProperlySet ();
  assertAssembled("accessing memory");
  CEPS_ABORT_IF (row < m_lo or row >= m_hi,
    "Trying to access row " << row << "on process that owns range " << m_lo << "to " << m_hi
  );
 
 
  rowValues.clear();
 
  CepsInt nbNonZero;
 
  const CepsMathScalar *vals;
  MatGetRow(m_A, row, &nbNonZero, nullptr, &vals);
 
  rowValues.reserve(nbNonZero);
  // copy the values
  rowValues.insert(rowValues.begin(),vals,vals+nbNonZero);
 
  // restore array
  MatRestoreRow(m_A, row, &nbNonZero, nullptr, &vals);
}
 
void
DistributedMatrix::getEntireLocalRow(CepsVector<CepsMathScalar> &rowValues, CepsGlobalIndex row)
  const
{
 
  checkProperlySet();
  assertAssembled("accessing memory");
  CEPS_ABORT_IF (row < m_lo or row >= m_hi,
    "Trying to access row " << row << " on process that owns range " << m_lo << " to " << m_hi
  );
 
 
  rowValues.clear();
  // resize it to global size
  CepsInt M, N;
  this->getSize(&M, &N);
  rowValues.resize(N, 0.0);
 
  CepsInt                nbNonZero;
  const CepsGlobalIndex *cols;
  const CepsMathScalar  *vals;
  MatGetRow(m_A, row, &nbNonZero, &cols, &vals);
 
  // iterate on non zero values
  for (CepsGlobalIndex i = 0; i < nbNonZero; i++)
    rowValues[cols[i]] = vals[i];
 
  // restore array
  MatRestoreRow(m_A, row, &nbNonZero, &cols, &vals);
  return;
}
 
void
DistributedMatrix::getColumnAsDistributedVector(DistributedVector &v, CepsGlobalIndex column) const
{
  // Compatibility check
  CepsGlobalIndex loV, loA, hiV, hiA;
  v.getLocalRange(&loV, &hiV);
  getLocalRange(&loA, &hiA);
  CEPS_ABORT_IF (loV != loA or hiV != hiA,
    "Incompatible decompositions of vector and matrix for column "
    << "insertion." << std::endl << "Matrix has local range (" << loA << "," << hiA
    << ") whereas vector has local range (" << loV << "," << hiV << ")"
  );
 
  CepsInt          nbCols     = 1;
  CepsInt          nbRows     = hiA - loA;
  CepsGlobalIndex *rowIndices = new CepsGlobalIndex[nbRows];
  for (CepsGlobalIndex i = loA; i < hiA; ++i)
    rowIndices[i - loA] = i;
  CepsGlobalIndex colIndices[1] = {column};
 
  v.getLocalData();
  MatGetValues(m_A, nbRows, rowIndices, nbCols, colIndices, v.m_localData);
  v.releaseLocalData();
  ceps::destroyTabular(rowIndices);
  return;
}
 
void
DistributedMatrix::setColumnAsDistributedVector(DistributedVector &v, CepsGlobalIndex column)
{
  // Compatibility check
  CepsGlobalIndex loV, loA, hiV, hiA;
  v.getLocalRange (&loV, &hiV);
  getLocalRange (&loA, &hiA);
  CEPS_ABORT_IF (loV != loA or hiV != hiA,
    "Incompatible decompositions of vector and matrix for column "
    << "insertion." << std::endl << "Matrix has local range (" << loA << "," << hiA
    << ") whereas vector has local range (" << loV << "," << hiV << ")"
  );
 
  CepsInt          nbCols     = 1;
  CepsInt          nbRows     = hiA - loA;
  CepsGlobalIndex *rowIndices = new CepsGlobalIndex[nbRows];
  for (CepsGlobalIndex i = loA; i < hiA; ++i)
    rowIndices[i - loA] = i;
  CepsGlobalIndex colIndices[1] = {column};
 
  v.getLocalData();
  MatSetValues(m_A, nbRows, rowIndices, nbCols, colIndices, v.m_localData, INSERT_VALUES);
  finalize();
  v.releaseLocalData();
  ceps::destroyTabular(rowIndices);
  return;
}
 
void 
DistributedMatrix::getDiagonalAsDistributedVector(DistributedVector& v) const
{
  // Compatibility check
  CepsGlobalIndex loV, loA, hiV, hiA;
  v.getLocalRange(&loV, &hiV);
  getLocalRange(&loA, &hiA);
  CEPS_ABORT_IF (loV != loA or hiV != hiA,
    "Incompatible decompositions of vector and matrix for column "
    << "insertion." << std::endl << "Matrix has local range (" << loA << "," << hiA
    << ") whereas vector has local range (" << loV << "," << hiV << ")"
  );
 
  v.getLocalData();
  MatGetDiagonal(m_A, v.m_v);
  v.releaseLocalData();
}
 
void 
DistributedMatrix::getDiagonalFootPrintAsDistributedVector(DistributedVector &v, CepsReal epsilon) const
{
  getDiagonalAsDistributedVector(v);
  CepsInt lo = 0, hi = 0;
  v.getLocalRange(&lo, &hi);
  v.getLocalData();
  for (CepsInt i = lo; i < hi; i++)
    v[i] = static_cast<CepsReal>(std::abs(v[i]) > epsilon);
  v.releaseLocalData();
}
 
void
DistributedMatrix::save(const CepsString &fileName)
{
  PetscViewer fd;
  PetscViewerBinaryOpen(ceps::getCommunicator (), fileName.c_str (), FILE_MODE_WRITE, &fd);
  MatView(m_A, fd);
  PetscViewerDestroy(&fd);
  return;
}
 
void
DistributedMatrix::load(const CepsString &fileName)
{
  PetscViewer fd;
  PetscViewerBinaryOpen(ceps::getCommunicator(),fileName.c_str(),FILE_MODE_READ,&fd);
  MatLoad(m_A,fd);
  PetscViewerDestroy(&fd);
 
  m_properlySet = true;
  return;
}
 
void
DistributedMatrix::beginAssembly()
{
  MatAssemblyBegin (m_A, MAT_FINAL_ASSEMBLY);
  return;
}
 
void
DistributedMatrix::endAssembly()
{
  MatAssemblyEnd(m_A, MAT_FINAL_ASSEMBLY);
  return;
}
 
void
DistributedMatrix::finalize()
{
  beginAssembly();
  endAssembly  ();
  return;
}
 
void
DistributedMatrix::flush()
{
  MatAssemblyBegin(m_A, MAT_FLUSH_ASSEMBLY);
  MatAssemblyEnd  (m_A, MAT_FLUSH_ASSEMBLY);
  return;
}
 
CepsBool
DistributedMatrix::isSymmetric() const
{
  PetscBool set, flag = PETSC_FALSE;
  MatIsSymmetricKnown(m_A, &set, &flag);
 
  return (flag == PETSC_FALSE ? false : true);
}
 
CepsBool
DistributedMatrix::isAssembled() const
{
  PetscBool assembled;
 
  MatAssembled(m_A, &assembled);
 
  return (assembled == PETSC_FALSE ? false : true);
}
 
void
DistributedMatrix::assertAssembled(CepsString info) const
{
  CepsString message = "Not to be used on a non-finalized matrix. ";
  if (not info.empty())
    message += "Call DistributedMatrix::finalize() before " + info;
  CEPS_ABORT_IF(not isAssembled(),message);
  return;
}
 
void
DistributedMatrix::view() const
{
  checkProperlySet();
  MatView(m_A, PETSC_VIEWER_STDOUT_WORLD);
  return;
}
 
void
DistributedMatrix::checkProperlySet() const
{
  CEPS_ABORT_IF (not m_properlySet, 
    "canot use this matrix as it is not properly set.\nCheck if "
    "sizes and memory preallocation are correct"
  );
  return;
}
 
DistributedMatrix& 
DistributedMatrix::operator*= (const CepsMathScalar &scalar)
{
  MatScale (m_A, scalar);
  return *this;
}