petsc_matrix_base.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.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/subscriptor.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 <petscmat.h>
 
#  include <cmath>
#  include <memory>
#  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
    {
    private:
      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 Subscriptor
  {
  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 std::vector<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:
    mutable std::vector<PetscInt> column_indices;
 
    mutable std::vector<PetscScalar> column_values;
 
 
    // To allow calling protected prepare_add() and prepare_set().
    template <class>
    friend class ::BlockMatrixBase;
  };
 
 
 
#  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)
  {
    prepare_action(VectorOperation::insert);
 
    const PetscInt  petsc_i = row;
    PetscInt const *col_index_ptr;
 
    PetscScalar const *col_value_ptr;
    int                n_columns;
 
    // If we don't elide zeros, the pointers are already available...
    if (elide_zero_values == false)
      {
        col_index_ptr = reinterpret_cast<const PetscInt *>(col_indices);
        col_value_ptr = values;
        n_columns     = n_cols;
      }
    else
      {
        // Otherwise, extract nonzero values in each row and get the
        // respective index.
        if (column_indices.size() < n_cols)
          {
            column_indices.resize(n_cols);
            column_values.resize(n_cols);
          }
 
        n_columns = 0;
        for (size_type j = 0; j < n_cols; ++j)
          {
            const PetscScalar value = values[j];
            AssertIsFinite(value);
            if (value != PetscScalar())
              {
                column_indices[n_columns] = col_indices[j];
                column_values[n_columns]  = value;
                n_columns++;
              }
          }
        AssertIndexRange(n_columns, n_cols + 1);
 
        col_index_ptr = column_indices.data();
        col_value_ptr = column_values.data();
      }
 
    const PetscErrorCode ierr = MatSetValues(matrix,
                                             1,
                                             &petsc_i,
                                             n_columns,
                                             col_index_ptr,
                                             col_value_ptr,
                                             INSERT_VALUES);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
  }
 
 
 
  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*/)
  {
    (void)elide_zero_values;
 
    prepare_action(VectorOperation::add);
 
    const PetscInt  petsc_i = row;
    PetscInt const *col_index_ptr;
 
    PetscScalar const *col_value_ptr;
    int                n_columns;
 
    // If we don't elide zeros, the pointers are already available...
    if (elide_zero_values == false)
      {
        col_index_ptr = reinterpret_cast<const PetscInt *>(col_indices);
        col_value_ptr = values;
        n_columns     = n_cols;
      }
    else
      {
        // Otherwise, extract nonzero values in each row and get the
        // respective index.
        if (column_indices.size() < n_cols)
          {
            column_indices.resize(n_cols);
            column_values.resize(n_cols);
          }
 
        n_columns = 0;
        for (size_type j = 0; j < n_cols; ++j)
          {
            const PetscScalar value = values[j];
            AssertIsFinite(value);
            if (value != PetscScalar())
              {
                column_indices[n_columns] = col_indices[j];
                column_values[n_columns]  = value;
                n_columns++;
              }
          }
        AssertIndexRange(n_columns, n_cols + 1);
 
        col_index_ptr = column_indices.data();
        col_value_ptr = column_values.data();
      }
 
    const PetscErrorCode ierr = MatSetValues(
      matrix, 1, &petsc_i, n_columns, col_index_ptr, col_value_ptr, 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 begin, end;
 
    const PetscErrorCode ierr =
      MatGetOwnershipRange(static_cast<const Mat &>(matrix), &begin, &end);
    AssertThrow(ierr == 0, ExcPETScError(ierr));
 
    return ((index >= static_cast<size_type>(begin)) &&
            (index < static_cast<size_type>(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