block_sparse_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2000 - 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_block_sparse_matrix_h
#define dealii_block_sparse_matrix_h


#include <deal.II/base/config.h>

#include <deal.II/lac/block_matrix_base.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_vector.h>
#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/sparse_matrix.h>

#include <cmath>

DEAL_II_NAMESPACE_OPEN


template <typename number>
class BlockSparseMatrix : public BlockMatrixBase<SparseMatrix<number>>
{
public:
  using BaseClass = BlockMatrixBase<SparseMatrix<number>>;

  using BlockType = typename BaseClass::BlockType;

  using value_type      = typename BaseClass::value_type;
  using pointer         = typename BaseClass::pointer;
  using const_pointer   = typename BaseClass::const_pointer;
  using reference       = typename BaseClass::reference;
  using const_reference = typename BaseClass::const_reference;
  using size_type       = typename BaseClass::size_type;
  using iterator        = typename BaseClass::iterator;
  using const_iterator  = typename BaseClass::const_iterator;

  BlockSparseMatrix() = default;

  BlockSparseMatrix(const BlockSparsityPattern &sparsity);

  virtual ~BlockSparseMatrix() override;



  BlockSparseMatrix &
  operator=(const BlockSparseMatrix &);

  BlockSparseMatrix &
  operator=(const double d);

  void
  clear();

  virtual void
  reinit(const BlockSparsityPattern &sparsity);

  bool
  empty() const;

  size_type
  get_row_length(const size_type row) const;

  size_type
  n_nonzero_elements() const;

  size_type
  n_actually_nonzero_elements(const double threshold = 0.0) const;

  const BlockSparsityPattern &
  get_sparsity_pattern() const;

  std::size_t
  memory_consumption() const;

  template <typename block_number>
  void
  vmult(BlockVector<block_number> &      dst,
        const BlockVector<block_number> &src) const;

  template <typename block_number, typename nonblock_number>
  void
  vmult(BlockVector<block_number> &    dst,
        const Vector<nonblock_number> &src) const;

  template <typename block_number, typename nonblock_number>
  void
  vmult(Vector<nonblock_number> &        dst,
        const BlockVector<block_number> &src) const;

  template <typename nonblock_number>
  void
  vmult(Vector<nonblock_number> &dst, const Vector<nonblock_number> &src) const;

  template <typename block_number>
  void
  Tvmult(BlockVector<block_number> &      dst,
         const BlockVector<block_number> &src) const;

  template <typename block_number, typename nonblock_number>
  void
  Tvmult(BlockVector<block_number> &    dst,
         const Vector<nonblock_number> &src) const;

  template <typename block_number, typename nonblock_number>
  void
  Tvmult(Vector<nonblock_number> &        dst,
         const BlockVector<block_number> &src) const;

  template <typename nonblock_number>
  void
  Tvmult(Vector<nonblock_number> &      dst,
         const Vector<nonblock_number> &src) const;

  template <class BlockVectorType>
  void
  precondition_Jacobi(BlockVectorType &      dst,
                      const BlockVectorType &src,
                      const number           omega = 1.) const;

  template <typename number2>
  void
  precondition_Jacobi(Vector<number2> &      dst,
                      const Vector<number2> &src,
                      const number           omega = 1.) const;

  void
  print_formatted(std::ostream &     out,
                  const unsigned int precision   = 3,
                  const bool         scientific  = true,
                  const unsigned int width       = 0,
                  const char *       zero_string = " ",
                  const double       denominator = 1.) const;

  DeclException0(ExcBlockDimensionMismatch);

private:
  SmartPointer<const BlockSparsityPattern, BlockSparseMatrix<number>>
    sparsity_pattern;
};



/* ------------------------- Template functions ---------------------- */



template <typename number>
inline BlockSparseMatrix<number> &
BlockSparseMatrix<number>::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;
}



template <typename number>
template <typename block_number>
inline void
BlockSparseMatrix<number>::vmult(BlockVector<block_number> &      dst,
                                 const BlockVector<block_number> &src) const
{
  BaseClass::vmult_block_block(dst, src);
}



template <typename number>
template <typename block_number, typename nonblock_number>
inline void
BlockSparseMatrix<number>::vmult(BlockVector<block_number> &    dst,
                                 const Vector<nonblock_number> &src) const
{
  BaseClass::vmult_block_nonblock(dst, src);
}



template <typename number>
template <typename block_number, typename nonblock_number>
inline void
BlockSparseMatrix<number>::vmult(Vector<nonblock_number> &        dst,
                                 const BlockVector<block_number> &src) const
{
  BaseClass::vmult_nonblock_block(dst, src);
}



template <typename number>
template <typename nonblock_number>
inline void
BlockSparseMatrix<number>::vmult(Vector<nonblock_number> &      dst,
                                 const Vector<nonblock_number> &src) const
{
  BaseClass::vmult_nonblock_nonblock(dst, src);
}



template <typename number>
template <typename block_number>
inline void
BlockSparseMatrix<number>::Tvmult(BlockVector<block_number> &      dst,
                                  const BlockVector<block_number> &src) const
{
  BaseClass::Tvmult_block_block(dst, src);
}



template <typename number>
template <typename block_number, typename nonblock_number>
inline void
BlockSparseMatrix<number>::Tvmult(BlockVector<block_number> &    dst,
                                  const Vector<nonblock_number> &src) const
{
  BaseClass::Tvmult_block_nonblock(dst, src);
}



template <typename number>
template <typename block_number, typename nonblock_number>
inline void
BlockSparseMatrix<number>::Tvmult(Vector<nonblock_number> &        dst,
                                  const BlockVector<block_number> &src) const
{
  BaseClass::Tvmult_nonblock_block(dst, src);
}



template <typename number>
template <typename nonblock_number>
inline void
BlockSparseMatrix<number>::Tvmult(Vector<nonblock_number> &      dst,
                                  const Vector<nonblock_number> &src) const
{
  BaseClass::Tvmult_nonblock_nonblock(dst, src);
}



template <typename number>
template <class BlockVectorType>
inline void
BlockSparseMatrix<number>::precondition_Jacobi(BlockVectorType &      dst,
                                               const BlockVectorType &src,
                                               const number omega) const
{
  Assert(this->n_block_rows() == this->n_block_cols(), ExcNotQuadratic());
  Assert(dst.n_blocks() == this->n_block_rows(),
         ExcDimensionMismatch(dst.n_blocks(), this->n_block_rows()));
  Assert(src.n_blocks() == this->n_block_cols(),
         ExcDimensionMismatch(src.n_blocks(), this->n_block_cols()));

  // do a diagonal preconditioning. uses only
  // the diagonal blocks of the matrix
  for (size_type i = 0; i < this->n_block_rows(); ++i)
    this->block(i, i).precondition_Jacobi(dst.block(i), src.block(i), omega);
}



template <typename number>
template <typename number2>
inline void
BlockSparseMatrix<number>::precondition_Jacobi(Vector<number2> &      dst,
                                               const Vector<number2> &src,
                                               const number omega) const
{
  // check number of blocks. the sizes of the
  // single block is checked in the function
  // we call
  Assert(this->n_block_cols() == 1,
         ExcMessage("This function only works if the matrix has "
                    "a single block"));
  Assert(this->n_block_rows() == 1,
         ExcMessage("This function only works if the matrix has "
                    "a single block"));

  // do a diagonal preconditioning. uses only
  // the diagonal blocks of the matrix
  this->block(0, 0).precondition_Jacobi(dst, src, omega);
}


DEAL_II_NAMESPACE_CLOSE

#endif // dealii_block_sparse_matrix_h