block_sparse_matrix.h

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


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

#include <cmath>

DEAL_II_NAMESPACE_OPEN


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

  typedef typename BaseClass::BlockType  BlockType;

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

  BlockSparseMatrix () = default;

  BlockSparseMatrix (const BlockSparsityPattern &sparsity);

  virtual ~BlockSparseMatrix ();



  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