trilinos_block_sparse_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2008 - 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_trilinos_block_sparse_matrix_h
#define dealii_trilinos_block_sparse_matrix_h
 
 
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_TRILINOS
 
#  include <deal.II/base/template_constraints.h>
 
#  include <deal.II/lac/block_matrix_base.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/full_matrix.h>
#  include <deal.II/lac/trilinos_parallel_block_vector.h>
#  include <deal.II/lac/trilinos_sparse_matrix.h>
 
#  include <cmath>
 
DEAL_II_NAMESPACE_OPEN
 
// forward declarations
class BlockSparsityPattern;
template <typename number>
class BlockSparseMatrix;
 
 
namespace TrilinosWrappers
{
  class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix>
  {
  public:
    using BaseClass = BlockMatrixBase<SparseMatrix>;
 
    using BlockType = BaseClass::BlockType;
 
    using value_type      = BaseClass::value_type;
    using pointer         = BaseClass::pointer;
    using const_pointer   = BaseClass::const_pointer;
    using reference       = BaseClass::reference;
    using const_reference = BaseClass::const_reference;
    using size_type       = BaseClass::size_type;
    using iterator        = BaseClass::iterator;
    using const_iterator  = BaseClass::const_iterator;
 
    BlockSparseMatrix() = default;
 
    ~BlockSparseMatrix() override;
 
    BlockSparseMatrix &
    operator=(const BlockSparseMatrix &) = default;
 
    BlockSparseMatrix &
    operator=(const double d);
 
    void
    reinit(const size_type n_block_rows, const size_type n_block_columns);
 
    template <typename BlockSparsityPatternType>
    void
    reinit(const std::vector<Epetra_Map> & input_maps,
           const BlockSparsityPatternType &block_sparsity_pattern,
           const bool                      exchange_data = false);
 
    template <typename BlockSparsityPatternType>
    void
    reinit(const std::vector<IndexSet> &   input_maps,
           const BlockSparsityPatternType &block_sparsity_pattern,
           const MPI_Comm &                communicator  = MPI_COMM_WORLD,
           const bool                      exchange_data = false);
 
    template <typename BlockSparsityPatternType>
    void
    reinit(const BlockSparsityPatternType &block_sparsity_pattern);
 
    DEAL_II_DEPRECATED
    void
    reinit(const std::vector<Epetra_Map> &            input_maps,
           const ::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
           const double                               drop_tolerance = 1e-13);
 
    void
    reinit(const ::BlockSparseMatrix<double> &deal_ii_sparse_matrix,
           const double                               drop_tolerance = 1e-13);
 
    bool
    is_compressed() const;
 
    void
    collect_sizes();
 
    size_type
    n_nonzero_elements() const;
 
    MPI_Comm
    get_mpi_communicator() const;
 
    DEAL_II_DEPRECATED
    std::vector<Epetra_Map>
    domain_partitioner() const;
 
    DEAL_II_DEPRECATED
    std::vector<Epetra_Map>
    range_partitioner() const;
 
    std::vector<IndexSet>
    locally_owned_domain_indices() const;
 
    std::vector<IndexSet>
    locally_owned_range_indices() const;
 
    template <typename VectorType1, typename VectorType2>
    void
    vmult(VectorType1 &dst, const VectorType2 &src) const;
 
    template <typename VectorType1, typename VectorType2>
    void
    Tvmult(VectorType1 &dst, const VectorType2 &src) const;
 
    TrilinosScalar
    residual(MPI::BlockVector &      dst,
             const MPI::BlockVector &x,
             const MPI::BlockVector &b) const;
 
    TrilinosScalar
    residual(MPI::BlockVector &      dst,
             const MPI::Vector &     x,
             const MPI::BlockVector &b) const;
 
    TrilinosScalar
    residual(MPI::Vector &           dst,
             const MPI::BlockVector &x,
             const MPI::Vector &     b) const;
 
    TrilinosScalar
    residual(MPI::Vector &      dst,
             const MPI::Vector &x,
             const MPI::Vector &b) const;
 
    using BlockMatrixBase<SparseMatrix>::clear;
 
    DeclException4(ExcIncompatibleRowNumbers,
                   int,
                   int,
                   int,
                   int,
                   << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
                   << ',' << arg4 << "] have differing row numbers.");
 
    DeclException4(ExcIncompatibleColNumbers,
                   int,
                   int,
                   int,
                   int,
                   << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
                   << ',' << arg4 << "] have differing column numbers.");
 
  private:
    template <typename VectorType1, typename VectorType2>
    void
    vmult(VectorType1 &      dst,
          const VectorType2 &src,
          const bool         transpose,
          const std::integral_constant<bool, true>,
          const std::integral_constant<bool, true>) const;
 
    template <typename VectorType1, typename VectorType2>
    void
    vmult(VectorType1 &      dst,
          const VectorType2 &src,
          const bool         transpose,
          const std::integral_constant<bool, false>,
          const std::integral_constant<bool, true>) const;
 
    template <typename VectorType1, typename VectorType2>
    void
    vmult(VectorType1 &      dst,
          const VectorType2 &src,
          const bool         transpose,
          const std::integral_constant<bool, true>,
          const std::integral_constant<bool, false>) const;
 
    template <typename VectorType1, typename VectorType2>
    void
    vmult(VectorType1 &      dst,
          const VectorType2 &src,
          const bool         transpose,
          const std::integral_constant<bool, false>,
          const std::integral_constant<bool, false>) const;
  };
 
 
 
  // ------------- inline and template functions -----------------
 
 
 
  inline BlockSparseMatrix &
  BlockSparseMatrix::operator=(const double d)
  {
    Assert(d == 0, ExcScalarAssignmentOnlyForZeroValue());
 
    for (size_type r = 0; r < this->n_block_rows(); ++r)
      for (size_type c = 0; c < this->n_block_cols(); ++c)
        this->block(r, c) = d;
 
    return *this;
  }
 
 
 
  inline bool
  BlockSparseMatrix::is_compressed() const
  {
    bool compressed = true;
    for (size_type row = 0; row < n_block_rows(); ++row)
      for (size_type col = 0; col < n_block_cols(); ++col)
        if (block(row, col).is_compressed() == false)
          {
            compressed = false;
            break;
          }
 
    return compressed;
  }
 
 
 
  template <typename VectorType1, typename VectorType2>
  inline void
  BlockSparseMatrix::vmult(VectorType1 &dst, const VectorType2 &src) const
  {
    vmult(dst,
          src,
          false,
          std::integral_constant<bool, IsBlockVector<VectorType1>::value>(),
          std::integral_constant<bool, IsBlockVector<VectorType2>::value>());
  }
 
 
 
  template <typename VectorType1, typename VectorType2>
  inline void
  BlockSparseMatrix::Tvmult(VectorType1 &dst, const VectorType2 &src) const
  {
    vmult(dst,
          src,
          true,
          std::integral_constant<bool, IsBlockVector<VectorType1>::value>(),
          std::integral_constant<bool, IsBlockVector<VectorType2>::value>());
  }
 
 
 
  template <typename VectorType1, typename VectorType2>
  inline void
  BlockSparseMatrix::vmult(VectorType1 &      dst,
                           const VectorType2 &src,
                           const bool         transpose,
                           std::integral_constant<bool, true>,
                           std::integral_constant<bool, true>) const
  {
    if (transpose == true)
      BaseClass::Tvmult_block_block(dst, src);
    else
      BaseClass::vmult_block_block(dst, src);
  }
 
 
 
  template <typename VectorType1, typename VectorType2>
  inline void
  BlockSparseMatrix::vmult(VectorType1 &      dst,
                           const VectorType2 &src,
                           const bool         transpose,
                           std::integral_constant<bool, false>,
                           std::integral_constant<bool, true>) const
  {
    if (transpose == true)
      BaseClass::Tvmult_nonblock_block(dst, src);
    else
      BaseClass::vmult_nonblock_block(dst, src);
  }
 
 
 
  template <typename VectorType1, typename VectorType2>
  inline void
  BlockSparseMatrix::vmult(VectorType1 &      dst,
                           const VectorType2 &src,
                           const bool         transpose,
                           std::integral_constant<bool, true>,
                           std::integral_constant<bool, false>) const
  {
    if (transpose == true)
      BaseClass::Tvmult_block_nonblock(dst, src);
    else
      BaseClass::vmult_block_nonblock(dst, src);
  }
 
 
 
  template <typename VectorType1, typename VectorType2>
  inline void
  BlockSparseMatrix::vmult(VectorType1 &      dst,
                           const VectorType2 &src,
                           const bool         transpose,
                           std::integral_constant<bool, false>,
                           std::integral_constant<bool, false>) const
  {
    if (transpose == true)
      BaseClass::Tvmult_nonblock_nonblock(dst, src);
    else
      BaseClass::vmult_nonblock_nonblock(dst, src);
  }
 
 
 
  inline std::vector<IndexSet>
  BlockSparseMatrix::locally_owned_domain_indices() const
  {
    Assert(this->n_block_cols() != 0, ExcNotInitialized());
    Assert(this->n_block_rows() != 0, ExcNotInitialized());
 
    std::vector<IndexSet> domain_indices;
    for (size_type c = 0; c < this->n_block_cols(); ++c)
      domain_indices.push_back(
        this->sub_objects[0][c]->locally_owned_domain_indices());
 
    return domain_indices;
  }
 
 
 
  inline std::vector<IndexSet>
  BlockSparseMatrix::locally_owned_range_indices() const
  {
    Assert(this->n_block_cols() != 0, ExcNotInitialized());
    Assert(this->n_block_rows() != 0, ExcNotInitialized());
 
    std::vector<IndexSet> range_indices;
    for (size_type r = 0; r < this->n_block_rows(); ++r)
      range_indices.push_back(
        this->sub_objects[r][0]->locally_owned_range_indices());
 
    return range_indices;
  }
 
 
 
  namespace internal
  {
    namespace BlockLinearOperatorImplementation
    {
      template <typename PayloadBlockType>
      class TrilinosBlockPayload
      {
      public:
        using BlockType = PayloadBlockType;
 
        template <typename... Args>
        TrilinosBlockPayload(const Args &...)
        {
          static_assert(
            std::is_same<
              PayloadBlockType,
              internal::LinearOperatorImplementation::TrilinosPayload>::value,
            "TrilinosBlockPayload can only accept a payload of type TrilinosPayload.");
        }
      };
 
    } // namespace BlockLinearOperatorImplementation
  }   /* namespace internal */
 
 
} /* namespace TrilinosWrappers */
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif // DEAL_II_WITH_TRILINOS
 
#endif // dealii_trilinos_block_sparse_matrix_h