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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/table.h>
#  include <deal.II/base/template_constraints.h>
#  include <deal.II/lac/block_matrix_base.h>
#  include <deal.II/lac/trilinos_sparse_matrix.h>
#  include <deal.II/lac/trilinos_parallel_block_vector.h>
#  include <deal.II/lac/full_matrix.h>
#  include <deal.II/lac/exceptions.h>

#  include <cmath>

DEAL_II_NAMESPACE_OPEN

// forward declarations
class BlockSparsityPattern;
template <typename number> class BlockSparseMatrix;


namespace TrilinosWrappers
{

  class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix>
  {
  public:
    typedef BlockMatrixBase<SparseMatrix> BaseClass;

    typedef BaseClass::BlockType  BlockType;

    typedef BaseClass::value_type      value_type;
    typedef BaseClass::pointer         pointer;
    typedef BaseClass::const_pointer   const_pointer;
    typedef BaseClass::reference       reference;
    typedef BaseClass::const_reference const_reference;
    typedef BaseClass::size_type       size_type;
    typedef BaseClass::iterator        iterator;
    typedef BaseClass::const_iterator  const_iterator;

    BlockSparseMatrix () = default;

    ~BlockSparseMatrix ();

    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:
        typedef PayloadBlockType BlockType;

        template <typename... Args>
        TrilinosBlockPayload (const Args &...)
        {
          static_assert(typeid(PayloadBlockType)==typeid(internal::LinearOperatorImplementation::TrilinosPayload),
                        "TrilinosBlockPayload can only accept a payload of type TrilinosPayload.");
        }
      };

    } /*namespace BlockLinearOperator*/
  } /* namespace internal */


}/* namespace TrilinosWrappers */


DEAL_II_NAMESPACE_CLOSE

#endif    // DEAL_II_WITH_TRILINOS

#endif    // dealii_trilinos_block_sparse_matrix_h