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

#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
// AffineConstraints::distribute_local_to_global)

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

#include <deal.II/base/exceptions.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/table.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:
  using size_type = types::global_dof_index;

  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