petsc_sparse_matrix.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_sparse_matrix_h
#define dealii_petsc_sparse_matrix_h
 
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_PETSC
 
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/petsc_matrix_base.h>
#  include <deal.II/lac/petsc_vector.h>
 
#  include <petscis.h>
#  include <petscistypes.h>
 
#  include <vector>
 
DEAL_II_NAMESPACE_OPEN
// forward declaration
#  ifndef DOXYGEN
template <typename MatrixType>
class BlockMatrixBase;
#  endif
 
namespace PETScWrappers
{
  class SparseMatrix : public MatrixBase
  {
  public:
    struct Traits
    {
      static const bool zero_addition_can_be_elided = true;
    };
 
    SparseMatrix();
 
    explicit SparseMatrix(const Mat &);
 
    SparseMatrix(const size_type m,
                 const size_type n,
                 const size_type n_nonzero_per_row,
                 const bool      is_symmetric = false);
 
    SparseMatrix(const size_type               m,
                 const size_type               n,
                 const std::vector<size_type> &row_lengths,
                 const bool                    is_symmetric = false);
 
    template <typename SparsityPatternType>
    explicit SparseMatrix(const SparsityPatternType &sparsity_pattern,
                          const bool preset_nonzero_locations = true);
 
    SparseMatrix &
    operator=(const double d);
 
    SparseMatrix(const SparseMatrix &) = delete;
 
    SparseMatrix &
    operator=(const SparseMatrix &) = delete;
 
    void
    reinit(const size_type m,
           const size_type n,
           const size_type n_nonzero_per_row,
           const bool      is_symmetric = false);
 
    void
    reinit(const size_type               m,
           const size_type               n,
           const std::vector<size_type> &row_lengths,
           const bool                    is_symmetric = false);
 
    template <typename SparsityPatternType>
    void
    reinit(const SparsityPatternType &sparsity_pattern,
           const bool                 preset_nonzero_locations = true);
 
    size_t
    m() const;
 
    size_t
    n() const;
 
    void
    mmult(SparseMatrix       &C,
          const SparseMatrix &B,
          const MPI::Vector  &V = MPI::Vector()) const;
 
    void
    Tmmult(SparseMatrix       &C,
           const SparseMatrix &B,
           const MPI::Vector  &V = MPI::Vector()) const;
 
  private:
    void
    do_reinit(const size_type m,
              const size_type n,
              const size_type n_nonzero_per_row,
              const bool      is_symmetric = false);
 
    void
    do_reinit(const size_type               m,
              const size_type               n,
              const std::vector<size_type> &row_lengths,
              const bool                    is_symmetric = false);
 
    template <typename SparsityPatternType>
    void
    do_reinit(const SparsityPatternType &sparsity_pattern,
              const bool                 preset_nonzero_locations);
 
    // To allow calling protected prepare_add() and prepare_set().
    friend class BlockMatrixBase<SparseMatrix>;
  };
 
  namespace MPI
  {
    class SparseMatrix : public MatrixBase
    {
    public:
      using size_type = types::global_dof_index;
 
      struct Traits
      {
        static const bool zero_addition_can_be_elided = false;
      };
 
      SparseMatrix();
 
      explicit SparseMatrix(const Mat &);
 
      ~SparseMatrix() override;
 
      template <typename SparsityPatternType>
      SparseMatrix(const MPI_Comm                communicator,
                   const SparsityPatternType    &sparsity_pattern,
                   const std::vector<size_type> &local_rows_per_process,
                   const std::vector<size_type> &local_columns_per_process,
                   const unsigned int            this_process,
                   const bool preset_nonzero_locations = true);
 
      SparseMatrix &
      operator=(const value_type d);
 
 
      void
      copy_from(const SparseMatrix &other);
 
      template <typename SparsityPatternType>
      void
      reinit(const MPI_Comm                communicator,
             const SparsityPatternType    &sparsity_pattern,
             const std::vector<size_type> &local_rows_per_process,
             const std::vector<size_type> &local_columns_per_process,
             const unsigned int            this_process,
             const bool                    preset_nonzero_locations = true);
 
      template <typename SparsityPatternType>
      void
      reinit(const IndexSet            &local_partitioning,
             const SparsityPatternType &sparsity_pattern,
             const MPI_Comm             communicator);
 
      template <typename SparsityPatternType>
      void
      reinit(const IndexSet            &local_rows,
             const IndexSet            &local_columns,
             const SparsityPatternType &sparsity_pattern,
             const MPI_Comm             communicator);
 
      void
      reinit(const SparseMatrix &other);
 
      template <typename SparsityPatternType>
      void
      reinit(const IndexSet            &local_rows,
             const IndexSet            &local_active_rows,
             const IndexSet            &local_columns,
             const IndexSet            &local_active_columns,
             const SparsityPatternType &sparsity_pattern,
             const MPI_Comm             communicator);
 
      DeclException2(ExcLocalRowsTooLarge,
                     int,
                     int,
                     << "The number of local rows " << arg1
                     << " must be larger than the total number of rows "
                     << arg2);
      PetscScalar
      matrix_norm_square(const Vector &v) const;
 
      PetscScalar
      matrix_scalar_product(const Vector &u, const Vector &v) const;
 
      IndexSet
      locally_owned_domain_indices() const;
 
      IndexSet
      locally_owned_range_indices() const;
 
      void
      mmult(SparseMatrix       &C,
            const SparseMatrix &B,
            const MPI::Vector  &V = MPI::Vector()) const;
 
      void
      Tmmult(SparseMatrix       &C,
             const SparseMatrix &B,
             const MPI::Vector  &V = MPI::Vector()) const;
 
    private:
      template <typename SparsityPatternType>
      void
      do_reinit(const MPI_Comm                comm,
                const SparsityPatternType    &sparsity_pattern,
                const std::vector<size_type> &local_rows_per_process,
                const std::vector<size_type> &local_columns_per_process,
                const unsigned int            this_process,
                const bool                    preset_nonzero_locations);
 
      template <typename SparsityPatternType>
      void
      do_reinit(const MPI_Comm             comm,
                const IndexSet            &local_rows,
                const IndexSet            &local_columns,
                const SparsityPatternType &sparsity_pattern);
 
      template <typename SparsityPatternType>
      void
      do_reinit(const MPI_Comm             comm,
                const IndexSet            &local_rows,
                const IndexSet            &local_active_rows,
                const IndexSet            &local_columns,
                const IndexSet            &local_active_columns,
                const SparsityPatternType &sparsity_pattern);
 
      // To allow calling protected prepare_add() and prepare_set().
      friend class BlockMatrixBase<SparseMatrix>;
    };
 
  } // namespace MPI
} // namespace PETScWrappers
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_PETSC
 
#endif