block_sparse_matrix_ez.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2002 - 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_ez_h
#define dealii_block_sparse_matrix_ez_h


//TODO: Derive BlockSparseMatrixEZ from BlockMatrixBase, like all the other block matrices as well; this would allow to instantiate a few functions with this template argument as well (in particular ConstraintMatrix::distribute_local_to_global)

#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/table.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/sparse_matrix_ez.h>

DEAL_II_NAMESPACE_OPEN

template <typename Number> class BlockVector;

template <typename Number>
class BlockSparseMatrixEZ : public Subscriptor
{
public:
  typedef types::global_dof_index size_type;

  BlockSparseMatrixEZ () = default;

  BlockSparseMatrixEZ (const unsigned int block_rows,
                       const unsigned int block_cols);

  BlockSparseMatrixEZ (const BlockSparseMatrixEZ<Number> &);

  BlockSparseMatrixEZ &operator = (const BlockSparseMatrixEZ<Number> &);

  BlockSparseMatrixEZ &operator = (const double d);


  void clear ();

  void reinit (const unsigned int n_block_rows,
               const unsigned int n_block_cols);
  void collect_sizes ();

  SparseMatrixEZ<Number> &
  block (const unsigned int row,
         const unsigned int column);


  const SparseMatrixEZ<Number> &
  block (const unsigned int row,
         const unsigned int column) const;

  unsigned int n_block_rows () const;

  unsigned int n_block_cols () const;

  bool empty () const;

  size_type m () const;

  size_type n () const;

  void set (const size_type i,
            const size_type j,
            const Number value);

  void add (const size_type i, const size_type j,
            const Number value);


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

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

  template <typename somenumber>
  void vmult_add (BlockVector<somenumber>       &dst,
                  const BlockVector<somenumber> &src) const;

  template <typename somenumber>
  void Tvmult_add (BlockVector<somenumber>       &dst,
                   const BlockVector<somenumber> &src) const;


  template <class StreamType>
  void print_statistics (StreamType &s, bool full = false);

private:
  BlockIndices    row_indices;

  BlockIndices    column_indices;

  Table<2, SparseMatrixEZ<Number> > blocks;
};

/*----------------------------------------------------------------------*/


template <typename Number>
inline
unsigned int
BlockSparseMatrixEZ<Number>::n_block_rows () const
{
  return row_indices.size();
}



template <typename Number>
inline
unsigned int
BlockSparseMatrixEZ<Number>::n_block_cols () const
{
  return column_indices.size();
}



template <typename Number>
inline
SparseMatrixEZ<Number> &
BlockSparseMatrixEZ<Number>::block (const unsigned int row,
                                    const unsigned int column)
{
  Assert (row<n_block_rows(), ExcIndexRange (row, 0, n_block_rows()));
  Assert (column<n_block_cols(), ExcIndexRange (column, 0, n_block_cols()));

  return blocks[row][column];
}



template <typename Number>
inline
const SparseMatrixEZ<Number> &
BlockSparseMatrixEZ<Number>::block (const unsigned int row,
                                    const unsigned int column) const
{
  Assert (row<n_block_rows(), ExcIndexRange (row, 0, n_block_rows()));
  Assert (column<n_block_cols(), ExcIndexRange (column, 0, n_block_cols()));

  return blocks[row][column];
}



template <typename Number>
inline
typename BlockSparseMatrixEZ<Number>::size_type
BlockSparseMatrixEZ<Number>::m () const
{
  return row_indices.total_size();
}



template <typename Number>
inline
typename BlockSparseMatrixEZ<Number>::size_type
BlockSparseMatrixEZ<Number>::n () const
{
  return column_indices.total_size();
}



template <typename Number>
inline
void
BlockSparseMatrixEZ<Number>::set (const size_type i,
                                  const size_type j,
                                  const Number value)
{

  AssertIsFinite(value);

  const std::pair<size_type,size_type>
  row_index = row_indices.global_to_local (i),
  col_index = column_indices.global_to_local (j);
  block(row_index.first,col_index.first).set (row_index.second,
                                              col_index.second,
                                              value);
}



template <typename Number>
inline
void
BlockSparseMatrixEZ<Number>::add (const size_type i,
                                  const size_type j,
                                  const Number value)
{

  AssertIsFinite(value);

  const std::pair<unsigned int,size_type>
  row_index = row_indices.global_to_local (i),
  col_index = column_indices.global_to_local (j);
  block(row_index.first,col_index.first).add (row_index.second,
                                              col_index.second,
                                              value);
}


template <typename Number>
template <typename somenumber>
void
BlockSparseMatrixEZ<Number>::vmult (BlockVector<somenumber>       &dst,
                                    const BlockVector<somenumber> &src) const
{
  Assert (dst.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert (src.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(src.n_blocks(), n_block_cols()));

  dst = 0.;

  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      block(row,col).vmult_add (dst.block(row),
                                src.block(col));
}



template <typename Number>
template <typename somenumber>
void
BlockSparseMatrixEZ<Number>::
vmult_add (BlockVector<somenumber>       &dst,
           const BlockVector<somenumber> &src) const
{
  Assert (dst.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert (src.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(src.n_blocks(), n_block_cols()));

  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      block(row,col).vmult_add (dst.block(row),
                                src.block(col));
}




template <typename Number>
template <typename somenumber>
void
BlockSparseMatrixEZ<Number>::
Tvmult (BlockVector<somenumber>       &dst,
        const BlockVector<somenumber> &src) const
{
  Assert (dst.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
  Assert (src.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(src.n_blocks(), n_block_rows()));

  dst = 0.;

  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      block(row,col).Tvmult_add (dst.block(col),
                                 src.block(row));
}



template <typename Number>
template <typename somenumber>
void
BlockSparseMatrixEZ<Number>::
Tvmult_add (BlockVector<somenumber>       &dst,
            const BlockVector<somenumber> &src) const
{
  Assert (dst.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
  Assert (src.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(src.n_blocks(), n_block_rows()));

  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      block(row,col).Tvmult_add (dst.block(col),
                                 src.block(row));
}


template <typename number>
template <class StreamType>
inline
void
BlockSparseMatrixEZ<number>::print_statistics (StreamType &out, bool full)
{
  size_type used_total = 0;
  size_type allocated_total = 0;
  size_type reserved_total = 0;
  std::vector<size_type> used_by_line_total;

  size_type used;
  size_type allocated;
  size_type reserved;
  std::vector<size_type> used_by_line;

  for (size_type i=0; i<n_block_rows(); ++i)
    for (size_type j=0; j<n_block_cols(); ++j)
      {
        used_by_line.clear();
        out << "block:\t" << i << '\t' << j << std::endl;
        block(i,j).compute_statistics (used, allocated, reserved,
                                       used_by_line, full);

        out << "used:" << used << std::endl
            << "allocated:" << allocated << std::endl
            << "reserved:" << reserved << std::endl;

        used_total += used;
        allocated_total += allocated;
        reserved_total += reserved;

        if (full)
          {
            used_by_line_total.resize(used_by_line.size());
            for (size_type i=0; i< used_by_line.size(); ++i)
              if (used_by_line[i] != 0)
                {
                  out << "row-entries\t" << i
                      << "\trows\t" << used_by_line[i]
                      << std::endl;
                  used_by_line_total[i] += used_by_line[i];
                }
          }
      }
  out << "Total" << std::endl
      << "used:" << used_total << std::endl
      << "allocated:" << allocated_total << std::endl
      << "reserved:" << reserved_total << std::endl;
  for (size_type i=0; i< used_by_line_total.size(); ++i)
    if (used_by_line_total[i] != 0)
      {
        out << "row-entries\t" << i
            << "\trows\t" << used_by_line_total[i]
            << std::endl;
      }
}


DEAL_II_NAMESPACE_CLOSE

#endif //dealii_block_sparse_matrix_ez_h