petsc_matrix_base.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2004 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_petsc_matrix_base_h
#define dealii_petsc_matrix_base_h
 
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_PETSC
 
#  include <deal.II/base/enable_observer_pointer.h>
 
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/full_matrix.h>
#  include <deal.II/lac/petsc_compatibility.h>
#  include <deal.II/lac/petsc_vector_base.h>
#  include <deal.II/lac/vector_operation.h>
 
#  include <boost/container/small_vector.hpp>
 
#  include <petscmat.h>
 
#  include <cmath>
#  include <memory>
#  include <optional>
#  include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#  ifndef DOXYGEN
template <typename Matrix>
class BlockMatrixBase;
#  endif
 
 
namespace PETScWrappers
{
  // forward declarations
  class MatrixBase;
 
  namespace MatrixIterators
  {
    class const_iterator
    {
#  ifdef __CUDA_ARCH__
      // NVCC, at least until 12.6, fails to compile the
      // implementations of the nested Accessor class if
      // it is declared as private. Work around this by
      // making it public.
    public:
#  else
    private:
#  endif
      class Accessor
      {
      public:
        using size_type = types::global_dof_index;
 
        Accessor(const MatrixBase *matrix,
                 const size_type   row,
                 const size_type   index);
 
        size_type
        row() const;
 
        size_type
        index() const;
 
        size_type
        column() const;
 
        PetscScalar
        value() const;
 
        DeclException0(ExcBeyondEndOfMatrix);
        DeclException3(ExcAccessToNonlocalRow,
                       int,
                       int,
                       int,
                       << "You tried to access row " << arg1
                       << " of a distributed matrix, but only rows " << arg2
                       << " through " << arg3
                       << " are stored locally and can be accessed.");
 
      private:
        mutable MatrixBase *matrix;
 
        size_type a_row;
 
        size_type a_index;
 
        std::shared_ptr<const std::vector<size_type>> colnum_cache;
 
        std::shared_ptr<const std::vector<PetscScalar>> value_cache;
 
        void
        visit_present_row();
 
        // Make enclosing class a friend.
        friend class const_iterator;
      };
 
    public:
      using size_type = types::global_dof_index;
 
      const_iterator(const MatrixBase *matrix,
                     const size_type   row,
                     const size_type   index);
 
      const_iterator &
      operator++();
 
      const_iterator
      operator++(int);
 
      const Accessor &
      operator*() const;
 
      const Accessor *
      operator->() const;
 
      bool
      operator==(const const_iterator &) const;
      bool
      operator!=(const const_iterator &) const;
 
      bool
      operator<(const const_iterator &) const;
 
      DeclException2(ExcInvalidIndexWithinRow,
                     int,
                     int,
                     << "Attempt to access element " << arg2 << " of row "
                     << arg1 << " which doesn't have that many elements.");
 
    private:
      Accessor accessor;
    };
 
  } // namespace MatrixIterators
 
 
  class MatrixBase : public EnableObserverPointer
  {
  public:
    using const_iterator = MatrixIterators::const_iterator;
 
    using size_type = types::global_dof_index;
 
    using value_type = PetscScalar;
 
    MatrixBase();
 
    explicit MatrixBase(const Mat &);
 
    MatrixBase(const MatrixBase &) = delete;
 
    MatrixBase &
    operator=(const MatrixBase &) = delete;
 
    virtual ~MatrixBase() override;
 
    void
    reinit(Mat A);
 
    MatrixBase &
    operator=(const value_type d);
 
    void
    clear();
 
    void
    set(const size_type i, const size_type j, const PetscScalar value);
 
    void
    set(const std::vector<size_type>  &indices,
        const FullMatrix<PetscScalar> &full_matrix,
        const bool                     elide_zero_values = false);
 
    void
    set(const std::vector<size_type>  &row_indices,
        const std::vector<size_type>  &col_indices,
        const FullMatrix<PetscScalar> &full_matrix,
        const bool                     elide_zero_values = false);
 
    void
    set(const size_type                 row,
        const std::vector<size_type>   &col_indices,
        const std::vector<PetscScalar> &values,
        const bool                      elide_zero_values = false);
 
    void
    set(const size_type    row,
        const size_type    n_cols,
        const size_type   *col_indices,
        const PetscScalar *values,
        const bool         elide_zero_values = false);
 
    void
    add(const size_type i, const size_type j, const PetscScalar value);
 
    void
    add(const std::vector<size_type>  &indices,
        const FullMatrix<PetscScalar> &full_matrix,
        const bool                     elide_zero_values = true);
 
    void
    add(const std::vector<size_type>  &row_indices,
        const std::vector<size_type>  &col_indices,
        const FullMatrix<PetscScalar> &full_matrix,
        const bool                     elide_zero_values = true);
 
    void
    add(const size_type                 row,
        const std::vector<size_type>   &col_indices,
        const std::vector<PetscScalar> &values,
        const bool                      elide_zero_values = true);
 
    void
    add(const size_type    row,
        const size_type    n_cols,
        const size_type   *col_indices,
        const PetscScalar *values,
        const bool         elide_zero_values      = true,
        const bool         col_indices_are_sorted = false);
 
    void
    clear_row(const size_type row, const PetscScalar new_diag_value = 0);
 
    void
    clear_rows(const ArrayView<const size_type> &rows,
               const PetscScalar                 new_diag_value = 0);
 
    void
    clear_rows_columns(const std::vector<size_type> &row_and_column_indices,
                       const PetscScalar             new_diag_value = 0);
 
    void
    compress(const VectorOperation::values operation);
 
    PetscScalar
    operator()(const size_type i, const size_type j) const;
 
    PetscScalar
    el(const size_type i, const size_type j) const;
 
    PetscScalar
    diag_element(const size_type i) const;
 
    size_type
    m() const;
 
    size_type
    n() const;
 
    size_type
    local_size() const;
 
    std::pair<size_type, size_type>
    local_range() const;
 
    bool
    in_local_range(const size_type index) const;
 
    size_type
    local_domain_size() const;
 
    std::pair<size_type, size_type>
    local_domain() const;
 
    MPI_Comm
    get_mpi_communicator() const;
 
    std::uint64_t
    n_nonzero_elements() const;
 
    size_type
    row_length(const size_type row) const;
 
    PetscReal
    l1_norm() const;
 
    PetscReal
    linfty_norm() const;
 
    PetscReal
    frobenius_norm() const;
 
 
    PetscScalar
    matrix_norm_square(const VectorBase &v) const;
 
 
    PetscScalar
    matrix_scalar_product(const VectorBase &u, const VectorBase &v) const;
 
    PetscScalar
    trace() const;
 
    MatrixBase &
    operator*=(const PetscScalar factor);
 
    MatrixBase &
    operator/=(const PetscScalar factor);
 
 
    MatrixBase &
    add(const PetscScalar factor, const MatrixBase &other);
 
    void
    vmult(VectorBase &dst, const VectorBase &src) const;
 
    void
    Tvmult(VectorBase &dst, const VectorBase &src) const;
 
    void
    vmult_add(VectorBase &dst, const VectorBase &src) const;
 
    void
    Tvmult_add(VectorBase &dst, const VectorBase &src) const;
 
    PetscScalar
    residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const;
 
    const_iterator
    begin() const;
 
    const_iterator
    end() const;
 
    const_iterator
    begin(const size_type r) const;
 
    const_iterator
    end(const size_type r) const;
 
    operator Mat() const;
 
    Mat &
    petsc_matrix();
 
    void
    transpose();
 
    PetscBool
    is_symmetric(const double tolerance = 1.e-12);
 
    PetscBool
    is_hermitian(const double tolerance = 1.e-12);
 
    void
    write_ascii(const PetscViewerFormat format = PETSC_VIEWER_DEFAULT);
 
    void
    print(std::ostream &out, const bool alternative_output = false) const;
 
    std::size_t
    memory_consumption() const;
 
    DeclExceptionMsg(ExcSourceEqualsDestination,
                     "You are attempting an operation on two vectors that "
                     "are the same object, but the operation requires that the "
                     "two objects are in fact different.");
 
    DeclException2(ExcWrongMode,
                   int,
                   int,
                   << "You tried to do a "
                   << (arg1 == 1 ? "'set'" : (arg1 == 2 ? "'add'" : "???"))
                   << " operation but the matrix is currently in "
                   << (arg2 == 1 ? "'set'" : (arg2 == 2 ? "'add'" : "???"))
                   << " mode. You first have to call 'compress()'.");
 
  protected:
    Mat matrix;
 
    VectorOperation::values last_action;
 
    void
    prepare_action(const VectorOperation::values new_action);
 
    void
    assert_is_compressed();
 
    void
    prepare_add();
    void
    prepare_set();
 
    void
    mmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const;
 
    void
    Tmmult(MatrixBase &C, const MatrixBase &B, const VectorBase &V) const;
 
  private:
    // To allow calling protected prepare_add() and prepare_set().
    template <class>
    friend class ::BlockMatrixBase;
 
 
    void
    add_or_set(const VectorOperation::values &operation,
               const size_type                row,
               const size_type                n_cols,
               const size_type               *col_indices,
               const PetscScalar             *values,
               const bool                     elide_zero_values);
  };
 
 
 
#  ifndef DOXYGEN
  // ---------------------- inline and template functions ---------------------
 
 
  namespace MatrixIterators
  {
    inline const_iterator::Accessor::Accessor(const MatrixBase *matrix,
                                              const size_type   row,
                                              const size_type   index)
      : matrix(const_cast<MatrixBase *>(matrix))
      , a_row(row)
      , a_index(index)
    {
      visit_present_row();
    }
 
 
 
    inline const_iterator::Accessor::size_type
    const_iterator::Accessor::row() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return a_row;
    }
 
 
    inline const_iterator::Accessor::size_type
    const_iterator::Accessor::column() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return (*colnum_cache)[a_index];
    }
 
 
    inline const_iterator::Accessor::size_type
    const_iterator::Accessor::index() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return a_index;
    }
 
 
    inline PetscScalar
    const_iterator::Accessor::value() const
    {
      Assert(a_row < matrix->m(), ExcBeyondEndOfMatrix());
      return (*value_cache)[a_index];
    }
 
 
    inline const_iterator::const_iterator(const MatrixBase *matrix,
                                          const size_type   row,
                                          const size_type   index)
      : accessor(matrix, row, index)
    {}
 
 
 
    inline const_iterator &
    const_iterator::operator++()
    {
      Assert(accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
 
      ++accessor.a_index;
 
      // if at end of line: do one step, then cycle until we find a
      // row with a nonzero number of entries
      if (accessor.a_index >= accessor.colnum_cache->size())
        {
          accessor.a_index = 0;
          ++accessor.a_row;
 
          while ((accessor.a_row < accessor.matrix->m()) &&
                 (accessor.a_row < accessor.matrix->local_range().second) &&
                 (accessor.matrix->row_length(accessor.a_row) == 0))
            ++accessor.a_row;
 
          accessor.visit_present_row();
        }
      return *this;
    }
 
 
    inline const_iterator
    const_iterator::operator++(int)
    {
      const const_iterator old_state = *this;
      ++(*this);
      return old_state;
    }
 
 
    inline const const_iterator::Accessor &
    const_iterator::operator*() const
    {
      return accessor;
    }
 
 
    inline const const_iterator::Accessor *
    const_iterator::operator->() const
    {
      return &accessor;
    }
 
 
    inline bool
    const_iterator::operator==(const const_iterator &other) const
    {
      return (accessor.a_row == other.accessor.a_row &&
              accessor.a_index == other.accessor.a_index);
    }
 
 
    inline bool
    const_iterator::operator!=(const const_iterator &other) const
    {
      return !(*this == other);
    }
 
 
    inline bool
    const_iterator::operator<(const const_iterator &other) const
    {
      return (accessor.row() < other.accessor.row() ||
              (accessor.row() == other.accessor.row() &&
               accessor.index() < other.accessor.index()));
    }
 
  } // namespace MatrixIterators
 
 
 
  // Inline the set() and add()
  // functions, since they will be
  // called frequently, and the
  // compiler can optimize away
  // some unnecessary loops when
  // the sizes are given at
  // compile time.
  inline void
  MatrixBase::set(const size_type i, const size_type j, const PetscScalar value)
  {
    AssertIsFinite(value);
 
    set(i, 1, &j, &value, false);
  }
 
 
 
  inline void
  MatrixBase::set(const std::vector<size_type>  &indices,
                  const FullMatrix<PetscScalar> &values,
                  const bool                     elide_zero_values)
  {
    Assert(indices.size() == values.m(),
           ExcDimensionMismatch(indices.size(), values.m()));
    Assert(values.m() == values.n(), ExcNotQuadratic());
 
    for (size_type i = 0; i < indices.size(); ++i)
      set(indices[i],
          indices.size(),
          indices.data(),
          &values(i, 0),
          elide_zero_values);
  }
 
 
 
  inline void
  MatrixBase::set(const std::vector<size_type>  &row_indices,
                  const std::vector<size_type>  &col_indices,
                  const FullMatrix<PetscScalar> &values,
                  const bool                     elide_zero_values)
  {
    Assert(row_indices.size() == values.m(),
           ExcDimensionMismatch(row_indices.size(), values.m()));
    Assert(col_indices.size() == values.n(),
           ExcDimensionMismatch(col_indices.size(), values.n()));
 
    for (size_type i = 0; i < row_indices.size(); ++i)
      set(row_indices[i],
          col_indices.size(),
          col_indices.data(),
          &values(i, 0),
          elide_zero_values);
  }
 
 
 
  inline void
  MatrixBase::set(const size_type                 row,
                  const std::vector<size_type>   &col_indices,
                  const std::vector<PetscScalar> &values,
                  const bool                      elide_zero_values)
  {
    Assert(col_indices.size() == values.size(),
           ExcDimensionMismatch(col_indices.size(), values.size()));
 
    set(row,
        col_indices.size(),
        col_indices.data(),
        values.data(),
        elide_zero_values);
  }
 
 
 
  inline void
  MatrixBase::set(const size_type    row,
                  const size_type    n_cols,
                  const size_type   *col_indices,
                  const PetscScalar *values,
                  const bool         elide_zero_values)
  {
    add_or_set(VectorOperation::insert,
               row,
               n_cols,
               col_indices,
               values,
               elide_zero_values);
  }
 
 
 
  inline void
  MatrixBase::add(const size_type i, const size_type j, const PetscScalar value)
  {
    AssertIsFinite(value);
 
    if (value == PetscScalar())
      {
        // we have to check after using Insert/Add in any case to be
        // consistent with the MPI communication model, but we can save
        // some work if the addend is zero. However, these actions are done
        // in case we pass on to the other function.
        prepare_action(VectorOperation::add);
 
        return;
      }
    else
      add(i, 1, &j, &value, false);
  }
 
 
 
  inline void
  MatrixBase::add(const std::vector<size_type>  &indices,
                  const FullMatrix<PetscScalar> &values,
                  const bool                     elide_zero_values)
  {
    Assert(indices.size() == values.m(),
           ExcDimensionMismatch(indices.size(), values.m()));
    Assert(values.m() == values.n(), ExcNotQuadratic());
 
    for (size_type i = 0; i < indices.size(); ++i)
      add(indices[i],
          indices.size(),
          indices.data(),
          &values(i, 0),
          elide_zero_values);
  }
 
 
 
  inline void
  MatrixBase::add(const std::vector<size_type>  &row_indices,
                  const std::vector<size_type>  &col_indices,
                  const FullMatrix<PetscScalar> &values,
                  const bool                     elide_zero_values)
  {
    Assert(row_indices.size() == values.m(),
           ExcDimensionMismatch(row_indices.size(), values.m()));
    Assert(col_indices.size() == values.n(),
           ExcDimensionMismatch(col_indices.size(), values.n()));
 
    for (size_type i = 0; i < row_indices.size(); ++i)
      add(row_indices[i],
          col_indices.size(),
          col_indices.data(),
          &values(i, 0),
          elide_zero_values);
  }
 
 
 
  inline void
  MatrixBase::add(const size_type                 row,
                  const std::vector<size_type>   &col_indices,
                  const std::vector<PetscScalar> &values,
                  const bool                      elide_zero_values)
  {
    Assert(col_indices.size() == values.size(),
           ExcDimensionMismatch(col_indices.size(), values.size()));
 
    add(row,
        col_indices.size(),
        col_indices.data(),
        values.data(),
        elide_zero_values);
  }
 
 
 
  inline void
  MatrixBase::add(const size_type    row,
                  const size_type    n_cols,
                  const size_type   *col_indices,
                  const PetscScalar *values,
                  const bool         elide_zero_values,
                  const bool /*col_indices_are_sorted*/)
  {
    add_or_set(VectorOperation::add,
               row,
               n_cols,
               col_indices,
               values,
               elide_zero_values);
  }
 
 
 
  inline void
  MatrixBase::add_or_set(const VectorOperation::values &operation,
                         const size_type                row,
                         const size_type                n_cols,
                         const size_type               *col_indices,
                         const PetscScalar             *values,
                         const bool                     elide_zero_values)
  {
    prepare_action(operation);
 
    const auto petsc_row = static_cast<PetscInt>(row);
    AssertIntegerConversion(petsc_row, row);
 
    // Use 100 entries so that we can store a row of a matrix constructed with
    // FESystem<3>(FE_Q<3>(2), 3, FE_Q<3>(1), 1) (i.e., a common Stokes FE,
    // which has 89 DoFs) with room to spare without allocating memory in the
    // free store.
    //
    // Setting up small_vectors isn't free so only set up column_values if we
    // actually use it.
    boost::container::small_vector<PetscInt, 100> column_indices;
    std::optional<boost::container::small_vector<PetscScalar, 100>>
      column_values;
 
    const PetscScalar *values_ptr = nullptr;
    if (elide_zero_values == false)
      {
        column_indices.resize(n_cols);
 
        for (size_type j = 0; j < n_cols; ++j)
          {
            AssertIsFinite(values[j]);
            column_indices[j] = static_cast<PetscInt>(col_indices[j]);
            AssertIntegerConversion(column_indices[j], col_indices[j]);
          }
 
        values_ptr = values;
      }
    else
      {
        // Otherwise, extract nonzero values in each row and get the
        // respective index.
        column_values.emplace();
        for (size_type j = 0; j < n_cols; ++j)
          {
            const PetscScalar value = values[j];
            AssertIsFinite(value);
            if (value != PetscScalar())
              {
                column_indices.push_back(static_cast<PetscInt>(col_indices[j]));
                AssertIntegerConversion(column_indices.back(), col_indices[j]);
                column_values->push_back(value);
              }
          }
        values_ptr = column_values->data();
      }
 
    const auto petsc_n_columns = static_cast<PetscInt>(column_indices.size());
    AssertIntegerConversion(petsc_n_columns, column_indices.size());
 
    Assert(operation == VectorOperation::insert ||
             operation == VectorOperation::add,
           ExcInternalError());
    const PetscErrorCode ierr =
      MatSetValues(matrix,
                   1,
                   &petsc_row,
                   petsc_n_columns,
                   column_indices.data(),
                   values_ptr,
                   operation == VectorOperation::insert ? INSERT_VALUES :
                                                          ADD_VALUES);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }
 
 
 
  inline PetscScalar
  MatrixBase::operator()(const size_type i, const size_type j) const
  {
    return el(i, j);
  }
 
 
 
  inline MatrixBase::const_iterator
  MatrixBase::begin() const
  {
    Assert(
      (in_local_range(0) && in_local_range(m() - 1)),
      ExcMessage(
        "begin() and end() can only be called on a processor owning the entire matrix. If this is a distributed matrix, use begin(row) and end(row) instead."));
 
    // find the first non-empty row in order to make sure that the returned
    // iterator points to something useful
    size_type first_nonempty_row = 0;
    while ((first_nonempty_row < m()) && (row_length(first_nonempty_row) == 0))
      ++first_nonempty_row;
 
    return const_iterator(this, first_nonempty_row, 0);
  }
 
 
  inline MatrixBase::const_iterator
  MatrixBase::end() const
  {
    Assert(
      (in_local_range(0) && in_local_range(m() - 1)),
      ExcMessage(
        "begin() and end() can only be called on a processor owning the entire matrix. If this is a distributed matrix, use begin(row) and end(row) instead."));
 
    return const_iterator(this, m(), 0);
  }
 
 
  inline MatrixBase::const_iterator
  MatrixBase::begin(const size_type r) const
  {
    Assert(in_local_range(r),
           ExcIndexRange(r, local_range().first, local_range().second));
 
    if (row_length(r) > 0)
      return const_iterator(this, r, 0);
    else
      return end(r);
  }
 
 
  inline MatrixBase::const_iterator
  MatrixBase::end(const size_type r) const
  {
    Assert(in_local_range(r),
           ExcIndexRange(r, local_range().first, local_range().second));
 
    // place the iterator on the first entry past this line, or at the
    // end of the matrix
    //
    // in the parallel case, we need to put it on the first entry of
    // the first row after the locally owned range. this of course
    // doesn't exist, but we can nevertheless create such an
    // iterator. we need to check whether 'i' is past the locally
    // owned range of rows first, before we ask for the length of the
    // row since the latter query leads to an exception in case we ask
    // for a row that is not locally owned
    for (size_type i = r + 1; i < m(); ++i)
      if (i == local_range().second || (row_length(i) > 0))
        return const_iterator(this, i, 0);
 
    // if there is no such line, then take the
    // end iterator of the matrix
    // we don't allow calling end() directly for distributed matrices so we need
    // to copy the code without the assertion.
    return {this, m(), 0};
  }
 
 
 
  inline bool
  MatrixBase::in_local_range(const size_type index) const
  {
    PetscInt petsc_begin, petsc_end;
 
    const PetscErrorCode ierr =
      MatGetOwnershipRange(static_cast<const Mat &>(matrix),
                           &petsc_begin,
                           &petsc_end);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
 
    const auto begin = static_cast<size_type>(petsc_begin);
    AssertIntegerConversion(begin, petsc_begin);
    const auto end = static_cast<size_type>(petsc_end);
    AssertIntegerConversion(end, petsc_end);
 
    return ((index >= begin) && (index < end));
  }
 
 
 
  inline void
  MatrixBase::prepare_action(const VectorOperation::values new_action)
  {
    if (last_action == VectorOperation::unknown)
      last_action = new_action;
 
    Assert(last_action == new_action, ExcWrongMode(last_action, new_action));
  }
 
 
 
  inline void
  MatrixBase::assert_is_compressed()
  {
    // compress() sets the last action to none, which allows us to check if
    // there are pending add/insert operations:
    AssertThrow(last_action == VectorOperation::unknown,
                ExcMessage("Error: missing compress() call."));
  }
 
 
 
  inline void
  MatrixBase::prepare_add()
  {
    prepare_action(VectorOperation::add);
  }
 
 
 
  inline void
  MatrixBase::prepare_set()
  {
    prepare_action(VectorOperation::insert);
  }
 
  inline MPI_Comm
  MatrixBase::get_mpi_communicator() const
  {
    return PetscObjectComm(reinterpret_cast<PetscObject>(matrix));
  }
 
#  endif // DOXYGEN
} // namespace PETScWrappers
 
 
DEAL_II_NAMESPACE_CLOSE
 
 
#endif // DEAL_II_WITH_PETSC
 
#endif