petsc_sparse_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2018 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_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 <vector>

DEAL_II_NAMESPACE_OPEN
// forward declaration
template <typename MatrixType>
class BlockMatrixBase;


namespace PETScWrappers
{
  class SparseMatrix : public MatrixBase
  {
  public:
    struct Traits
    {
      static const bool zero_addition_can_be_elided = true;
    };

    SparseMatrix();

    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);

    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);

    virtual const MPI_Comm &
    get_mpi_communicator() const override;

    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:
    SparseMatrix(const SparseMatrix &) = delete;
    SparseMatrix &
    operator=(const SparseMatrix &) = delete;

    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);

    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();

      ~SparseMatrix() override;

      DEAL_II_DEPRECATED
      SparseMatrix(const MPI_Comm &communicator,
                   const size_type m,
                   const size_type n,
                   const size_type local_rows,
                   const size_type local_columns,
                   const size_type n_nonzero_per_row,
                   const bool      is_symmetric              = false,
                   const size_type n_offdiag_nonzero_per_row = 0);

      DEAL_II_DEPRECATED
      SparseMatrix(const MPI_Comm &              communicator,
                   const size_type               m,
                   const size_type               n,
                   const size_type               local_rows,
                   const size_type               local_columns,
                   const std::vector<size_type> &row_lengths,
                   const bool                    is_symmetric = false,
                   const std::vector<size_type> &offdiag_row_lengths =
                     std::vector<size_type>());

      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);

      DEAL_II_DEPRECATED
      void
      reinit(const MPI_Comm &communicator,
             const size_type m,
             const size_type n,
             const size_type local_rows,
             const size_type local_columns,
             const size_type n_nonzero_per_row,
             const bool      is_symmetric              = false,
             const size_type n_offdiag_nonzero_per_row = 0);

      DEAL_II_DEPRECATED
      void
      reinit(const MPI_Comm &              communicator,
             const size_type               m,
             const size_type               n,
             const size_type               local_rows,
             const size_type               local_columns,
             const std::vector<size_type> &row_lengths,
             const bool                    is_symmetric = false,
             const std::vector<size_type> &offdiag_row_lengths =
               std::vector<size_type>());

      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_rows,
             const IndexSet &           local_columns,
             const SparsityPatternType &sparsity_pattern,
             const MPI_Comm &           communicator);

      void
      reinit(const SparseMatrix &other);

      virtual const MPI_Comm &
      get_mpi_communicator() const override;

      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:
      MPI_Comm communicator;

      DEAL_II_DEPRECATED
      void
      do_reinit(const size_type m,
                const size_type n,
                const size_type local_rows,
                const size_type local_columns,
                const size_type n_nonzero_per_row,
                const bool      is_symmetric              = false,
                const size_type n_offdiag_nonzero_per_row = 0);

      DEAL_II_DEPRECATED
      void
      do_reinit(const size_type               m,
                const size_type               n,
                const size_type               local_rows,
                const size_type               local_columns,
                const std::vector<size_type> &row_lengths,
                const bool                    is_symmetric = false,
                const std::vector<size_type> &offdiag_row_lengths =
                  std::vector<size_type>());

      template <typename SparsityPatternType>
      void
      do_reinit(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 IndexSet &           local_rows,
                const IndexSet &           local_columns,
                const SparsityPatternType &sparsity_pattern);

      friend class BlockMatrixBase<SparseMatrix>;
    };



    // -------- template and inline functions ----------

    inline const MPI_Comm &
    SparseMatrix::get_mpi_communicator() const
    {
      return communicator;
    }
  } // namespace MPI
} // namespace PETScWrappers

DEAL_II_NAMESPACE_CLOSE

#  endif // DEAL_II_WITH_PETSC

#endif
/*--------------------------- petsc_sparse_matrix.h -------------------------*/