matrix_block.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2007 - 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_matrix_block_h
#define dealii_matrix_block_h

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

#include <deal.II/algorithms/any_data.h>

#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/base/smartpointer.h>

#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>

#include <memory>


DEAL_II_NAMESPACE_OPEN

template <typename MatrixType>
class MatrixBlock;

namespace internal
{
  template <typename MatrixType>
  void
  reinit(MatrixBlock<MatrixType> &v, const BlockSparsityPattern &p);

  template <typename number>
  void
  reinit(MatrixBlock<::SparseMatrix<number>> &v,
         const BlockSparsityPattern &               p);
} // namespace internal

template <typename MatrixType>
class MatrixBlock : public Subscriptor
{
public:
  using size_type = types::global_dof_index;

  using value_type = typename MatrixType::value_type;

  MatrixBlock();

  MatrixBlock(const MatrixBlock<MatrixType> &M);

  MatrixBlock(size_type i, size_type j);

  void
  reinit(const BlockSparsityPattern &sparsity);

  operator MatrixType &();
  operator const MatrixType &() const;

  void
  add(const size_type                       i,
      const size_type                       j,
      const typename MatrixType::value_type value);

  template <typename number>
  void
  add(const std::vector<size_type> &indices,
      const FullMatrix<number> &    full_matrix,
      const bool                    elide_zero_values = true);

  template <typename number>
  void
  add(const std::vector<size_type> &row_indices,
      const std::vector<size_type> &col_indices,
      const FullMatrix<number> &    full_matrix,
      const bool                    elide_zero_values = true);

  template <typename number>
  void
  add(const size_type               row_index,
      const std::vector<size_type> &col_indices,
      const std::vector<number> &   values,
      const bool                    elide_zero_values = true);

  template <typename number>
  void
  add(const size_type  row,
      const size_type  n_cols,
      const size_type *col_indices,
      const number *   values,
      const bool       elide_zero_values      = true,
      const bool       col_indices_are_sorted = false);

  template <class VectorType>
  void
  vmult(VectorType &w, const VectorType &v) const;

  template <class VectorType>
  void
  vmult_add(VectorType &w, const VectorType &v) const;

  template <class VectorType>
  void
  Tvmult(VectorType &w, const VectorType &v) const;

  template <class VectorType>
  void
  Tvmult_add(VectorType &w, const VectorType &v) const;

  std::size_t
  memory_consumption() const;

  DeclException2(ExcBlockIndexMismatch,
                 size_type,
                 size_type,
                 << "Block index " << arg1 << " does not match " << arg2);

  size_type row;
  size_type column;

  MatrixType matrix;

private:
  BlockIndices row_indices;
  BlockIndices column_indices;

  template <class OTHER_MatrixType>
  friend void
  ::internal::reinit(MatrixBlock<OTHER_MatrixType> &,
                           const BlockSparsityPattern &);

  template <typename number>
  friend void
  internal::reinit(MatrixBlock<::SparseMatrix<number>> &v,
                   const BlockSparsityPattern &               p);
};


template <typename MatrixType>
class MatrixBlockVector : private AnyData
{
public:
  using size_type = types::global_dof_index;

  using value_type = MatrixBlock<MatrixType>;

  using ptr_type = std::shared_ptr<value_type>;

  void
  add(size_type row, size_type column, const std::string &name);

  void
  reinit(const BlockSparsityPattern &sparsity);

  void
  clear(bool really_clean = false);

  std::size_t
  memory_consumption() const;

  const value_type &
  block(size_type i) const;

  value_type &
  block(size_type i);

  MatrixType &
  matrix(size_type i);

  using AnyData::name;
  using AnyData::size;
  using AnyData::subscribe;
  using AnyData::unsubscribe;
};


template <typename MatrixType>
class MGMatrixBlockVector : public Subscriptor
{
public:
  using size_type = types::global_dof_index;

  using value_type = MGLevelObject<MatrixBlock<MatrixType>>;
  MGMatrixBlockVector(const bool edge_matrices      = false,
                      const bool edge_flux_matrices = false);

  unsigned int
  size() const;

  void
  add(size_type row, size_type column, const std::string &name);

  void
  reinit_matrix(const MGLevelObject<BlockSparsityPattern> &sparsity);
  void
  reinit_edge(const MGLevelObject<BlockSparsityPattern> &sparsity);
  void
  reinit_edge_flux(const MGLevelObject<BlockSparsityPattern> &sparsity);

  void
  clear(bool really_clean = false);

  const value_type &
  block(size_type i) const;

  value_type &
  block(size_type i);

  const value_type &
  block_in(size_type i) const;

  value_type &
  block_in(size_type i);

  const value_type &
  block_out(size_type i) const;

  value_type &
  block_out(size_type i);

  const value_type &
  block_up(size_type i) const;

  value_type &
  block_up(size_type i);

  const value_type &
  block_down(size_type i) const;

  value_type &
  block_down(size_type i);

  std::size_t
  memory_consumption() const;

private:
  void
  clear_object(AnyData &);

  const bool edge_matrices;

  const bool edge_flux_matrices;

  AnyData matrices;
  AnyData matrices_in;
  AnyData matrices_out;
  AnyData flux_matrices_down;
  AnyData flux_matrices_up;
};


//----------------------------------------------------------------------//

namespace internal
{
  template <typename MatrixType>
  void
  reinit(MatrixBlock<MatrixType> &v, const BlockSparsityPattern &p)
  {
    v.row_indices    = p.get_row_indices();
    v.column_indices = p.get_column_indices();
  }


  template <typename number>
  void
  reinit(MatrixBlock<::SparseMatrix<number>> &v,
         const BlockSparsityPattern &               p)
  {
    v.row_indices    = p.get_row_indices();
    v.column_indices = p.get_column_indices();
    v.matrix.reinit(p.block(v.row, v.column));
  }
} // namespace internal


template <typename MatrixType>
inline MatrixBlock<MatrixType>::MatrixBlock()
  : row(numbers::invalid_size_type)
  , column(numbers::invalid_size_type)
{}


template <typename MatrixType>
inline MatrixBlock<MatrixType>::MatrixBlock(const MatrixBlock<MatrixType> &M)
  : Subscriptor()
  , row(M.row)
  , column(M.column)
  , matrix(M.matrix)
  , row_indices(M.row_indices)
  , column_indices(M.column_indices)
{}


template <typename MatrixType>
inline MatrixBlock<MatrixType>::MatrixBlock(size_type i, size_type j)
  : row(i)
  , column(j)
{}


template <typename MatrixType>
inline void
MatrixBlock<MatrixType>::reinit(const BlockSparsityPattern &sparsity)
{
  internal::reinit(*this, sparsity);
}


template <typename MatrixType>
inline MatrixBlock<MatrixType>::operator MatrixType &()
{
  return matrix;
}


template <typename MatrixType>
inline MatrixBlock<MatrixType>::operator const MatrixType &() const
{
  return matrix;
}


template <typename MatrixType>
inline void
MatrixBlock<MatrixType>::add(const size_type                       gi,
                             const size_type                       gj,
                             const typename MatrixType::value_type value)
{
  Assert(row_indices.size() != 0, ExcNotInitialized());
  Assert(column_indices.size() != 0, ExcNotInitialized());

  const std::pair<unsigned int, size_type> bi = row_indices.global_to_local(gi);
  const std::pair<unsigned int, size_type> bj =
    column_indices.global_to_local(gj);

  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));

  matrix.add(bi.second, bj.second, value);
}


template <typename MatrixType>
template <typename number>
inline void
MatrixBlock<MatrixType>::add(const std::vector<size_type> &r_indices,
                             const std::vector<size_type> &c_indices,
                             const FullMatrix<number> &    values,
                             const bool                    elide_zero_values)
{
  Assert(row_indices.size() != 0, ExcNotInitialized());
  Assert(column_indices.size() != 0, ExcNotInitialized());

  AssertDimension(r_indices.size(), values.m());
  AssertDimension(c_indices.size(), values.n());

  for (size_type i = 0; i < row_indices.size(); ++i)
    add(r_indices[i],
        c_indices.size(),
        c_indices.data(),
        &values(i, 0),
        elide_zero_values);
}


template <typename MatrixType>
template <typename number>
inline void
MatrixBlock<MatrixType>::add(const size_type  b_row,
                             const size_type  n_cols,
                             const size_type *col_indices,
                             const number *   values,
                             const bool,
                             const bool)
{
  Assert(row_indices.size() != 0, ExcNotInitialized());
  Assert(column_indices.size() != 0, ExcNotInitialized());

  const std::pair<unsigned int, size_type> bi =
    row_indices.global_to_local(b_row);

  // In debug mode, we check whether
  // all indices are in the correct
  // block.

  // Actually, for the time being, we
  // leave it at this. While it may
  // not be the most efficient way,
  // it is at least thread safe.
  //#ifdef DEBUG
  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));

  for (size_type j = 0; j < n_cols; ++j)
    {
      const std::pair<unsigned int, size_type> bj =
        column_indices.global_to_local(col_indices[j]);
      Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));

      matrix.add(bi.second, bj.second, values[j]);
    }
  //#endif
}


template <typename MatrixType>
template <typename number>
inline void
MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
                             const FullMatrix<number> &    values,
                             const bool                    elide_zero_values)
{
  Assert(row_indices.size() != 0, ExcNotInitialized());
  Assert(column_indices.size() != 0, ExcNotInitialized());

  AssertDimension(indices.size(), values.m());
  Assert(values.n() == values.m(), ExcNotQuadratic());

  for (size_type i = 0; i < indices.size(); ++i)
    add(indices[i],
        indices.size(),
        indices.data(),
        &values(i, 0),
        elide_zero_values);
}



template <typename MatrixType>
template <typename number>
inline void
MatrixBlock<MatrixType>::add(const size_type               row,
                             const std::vector<size_type> &col_indices,
                             const std::vector<number> &   values,
                             const bool                    elide_zero_values)
{
  Assert(row_indices.size() != 0, ExcNotInitialized());
  Assert(column_indices.size() != 0, ExcNotInitialized());

  AssertDimension(col_indices.size(), values.size());
  add(row,
      col_indices.size(),
      col_indices.data(),
      values.data(),
      elide_zero_values);
}


template <typename MatrixType>
template <class VectorType>
inline void
MatrixBlock<MatrixType>::vmult(VectorType &w, const VectorType &v) const
{
  matrix.vmult(w, v);
}


template <typename MatrixType>
template <class VectorType>
inline void
MatrixBlock<MatrixType>::vmult_add(VectorType &w, const VectorType &v) const
{
  matrix.vmult_add(w, v);
}


template <typename MatrixType>
template <class VectorType>
inline void
MatrixBlock<MatrixType>::Tvmult(VectorType &w, const VectorType &v) const
{
  matrix.Tvmult(w, v);
}


template <typename MatrixType>
template <class VectorType>
inline void
MatrixBlock<MatrixType>::Tvmult_add(VectorType &w, const VectorType &v) const
{
  matrix.Tvmult_add(w, v);
}


template <typename MatrixType>
inline std::size_t
MatrixBlock<MatrixType>::memory_consumption() const
{
  return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
          sizeof(matrix));
}

//----------------------------------------------------------------------//

template <typename MatrixType>
inline void
MatrixBlockVector<MatrixType>::add(size_type          row,
                                   size_type          column,
                                   const std::string &name)
{
  ptr_type p(new value_type(row, column));
  AnyData::add(p, name);
}


template <typename MatrixType>
inline void
MatrixBlockVector<MatrixType>::reinit(const BlockSparsityPattern &sparsity)
{
  for (size_type i = 0; i < this->size(); ++i)
    {
      block(i).reinit(sparsity);
    }
}


template <typename MatrixType>
inline void
MatrixBlockVector<MatrixType>::clear(bool really_clean)
{
  if (really_clean)
    {
      Assert(false, ExcNotImplemented());
    }
  else
    {
      for (size_type i = 0; i < this->size(); ++i)
        matrix(i).clear();
    }
}



template <typename MatrixType>
inline const MatrixBlock<MatrixType> &
MatrixBlockVector<MatrixType>::block(size_type i) const
{
  return *this->read<ptr_type>(i);
}


template <typename MatrixType>
inline MatrixBlock<MatrixType> &
MatrixBlockVector<MatrixType>::block(size_type i)
{
  return *this->entry<ptr_type>(i);
}


template <typename MatrixType>
inline MatrixType &
MatrixBlockVector<MatrixType>::matrix(size_type i)
{
  return this->entry<ptr_type>(i)->matrix;
}



//----------------------------------------------------------------------//

template <typename MatrixType>
inline MGMatrixBlockVector<MatrixType>::MGMatrixBlockVector(const bool e,
                                                            const bool f)
  : edge_matrices(e)
  , edge_flux_matrices(f)
{}


template <typename MatrixType>
inline unsigned int
MGMatrixBlockVector<MatrixType>::size() const
{
  return matrices.size();
}


template <typename MatrixType>
inline void
MGMatrixBlockVector<MatrixType>::add(size_type          row,
                                     size_type          column,
                                     const std::string &name)
{
  MGLevelObject<MatrixBlock<MatrixType>> p(0, 1);
  p[0].row    = row;
  p[0].column = column;

  matrices.add(p, name);
  if (edge_matrices)
    {
      matrices_in.add(p, name);
      matrices_out.add(p, name);
    }
  if (edge_flux_matrices)
    {
      flux_matrices_up.add(p, name);
      flux_matrices_down.add(p, name);
    }
}


template <typename MatrixType>
inline const MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block(size_type i) const
{
  return *matrices.read<const MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block(size_type i)
{
  return *matrices.entry<MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline const MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block_in(size_type i) const
{
  return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block_in(size_type i)
{
  return *matrices_in.entry<MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline const MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block_out(size_type i) const
{
  return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block_out(size_type i)
{
  return *matrices_out.entry<MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline const MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block_up(size_type i) const
{
  return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block_up(size_type i)
{
  return *flux_matrices_up.entry<MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline const MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block_down(size_type i) const
{
  return *flux_matrices_down.read<const MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline MGLevelObject<MatrixBlock<MatrixType>> &
MGMatrixBlockVector<MatrixType>::block_down(size_type i)
{
  return *flux_matrices_down.entry<MGLevelObject<MatrixType> *>(i);
}


template <typename MatrixType>
inline void
MGMatrixBlockVector<MatrixType>::reinit_matrix(
  const MGLevelObject<BlockSparsityPattern> &sparsity)
{
  for (size_type i = 0; i < this->size(); ++i)
    {
      MGLevelObject<MatrixBlock<MatrixType>> &o   = block(i);
      const size_type                         row = o[o.min_level()].row;
      const size_type                         col = o[o.min_level()].column;

      o.resize(sparsity.min_level(), sparsity.max_level());
      for (size_type level = o.min_level(); level <= o.max_level(); ++level)
        {
          o[level].row    = row;
          o[level].column = col;
          internal::reinit(o[level], sparsity[level]);
        }
    }
}


template <typename MatrixType>
inline void
MGMatrixBlockVector<MatrixType>::reinit_edge(
  const MGLevelObject<BlockSparsityPattern> &sparsity)
{
  for (size_type i = 0; i < this->size(); ++i)
    {
      MGLevelObject<MatrixBlock<MatrixType>> &o   = block(i);
      const size_type                         row = o[o.min_level()].row;
      const size_type                         col = o[o.min_level()].column;

      block_in(i).resize(sparsity.min_level(), sparsity.max_level());
      block_out(i).resize(sparsity.min_level(), sparsity.max_level());
      for (size_type level = o.min_level(); level <= o.max_level(); ++level)
        {
          block_in(i)[level].row    = row;
          block_in(i)[level].column = col;
          internal::reinit(block_in(i)[level], sparsity[level]);
          block_out(i)[level].row    = row;
          block_out(i)[level].column = col;
          internal::reinit(block_out(i)[level], sparsity[level]);
        }
    }
}


template <typename MatrixType>
inline void
MGMatrixBlockVector<MatrixType>::reinit_edge_flux(
  const MGLevelObject<BlockSparsityPattern> &sparsity)
{
  for (size_type i = 0; i < this->size(); ++i)
    {
      MGLevelObject<MatrixBlock<MatrixType>> &o   = block(i);
      const size_type                         row = o[o.min_level()].row;
      const size_type                         col = o[o.min_level()].column;

      block_up(i).resize(sparsity.min_level(), sparsity.max_level());
      block_down(i).resize(sparsity.min_level(), sparsity.max_level());
      for (size_type level = o.min_level(); level <= o.max_level(); ++level)
        {
          block_up(i)[level].row    = row;
          block_up(i)[level].column = col;
          internal::reinit(block_up(i)[level], sparsity[level]);
          block_down(i)[level].row    = row;
          block_down(i)[level].column = col;
          internal::reinit(block_down(i)[level], sparsity[level]);
        }
    }
}


template <typename MatrixType>
inline void
MGMatrixBlockVector<MatrixType>::clear_object(AnyData &mo)
{
  for (size_type i = 0; i < mo.size(); ++i)
    {
      MGLevelObject<MatrixBlock<MatrixType>> &o =
        mo.entry<MGLevelObject<MatrixType> *>(i);
      for (size_type level = o.min_level(); level <= o.max_level(); ++level)
        o[level].matrix.clear();
    }
}


template <typename MatrixType>
inline void
MGMatrixBlockVector<MatrixType>::clear(bool really_clean)
{
  if (really_clean)
    {
      Assert(false, ExcNotImplemented());
    }
  else
    {
      clear_object(matrices);
      clear_object(matrices_in);
      clear_object(matrices_out);
      clear_object(flux_matrices_up);
      clear_object(flux_matrices_down);
    }
}



DEAL_II_NAMESPACE_CLOSE

#endif