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

#ifndef dealii_matrix_block_h
#define dealii_matrix_block_h

#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/mg_level_object.h>
#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/algorithms/any_data.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);
}

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

  typedef typename MatrixType::value_type 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:
  typedef types::global_dof_index size_type;

  typedef MatrixBlock<MatrixType> value_type;

  typedef std::shared_ptr<value_type> ptr_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::subscribe;
  using AnyData::unsubscribe;
  using AnyData::size;
  using AnyData::name;
};


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

  typedef MGLevelObject<MatrixBlock<MatrixType> > value_type;
  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));
  }
}


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