block_sparse_matrix_ez.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2002 - 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_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