filtered_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 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_filtered_matrix_h
#  define dealii_filtered_matrix_h



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

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

#  include <deal.II/lac/pointer_matrix.h>
#  include <deal.II/lac/vector_memory.h>

#  include <algorithm>
#  include <vector>

DEAL_II_NAMESPACE_OPEN

template <class VectorType>
class FilteredMatrixBlock;

template <typename VectorType>
class DEAL_II_DEPRECATED FilteredMatrix : public Subscriptor
{
public:
  class const_iterator;

  using size_type = types::global_dof_index;

  class Accessor
  {
    Accessor(const FilteredMatrix<VectorType> *matrix, const size_type index);

  public:
    size_type
    row() const;

    size_type
    column() const;

    double
    value() const;

  private:
    void
    advance();

    const FilteredMatrix<VectorType> *matrix;

    size_type index;
    /*
     * Make enclosing class a
     * friend.
     */
    friend class const_iterator;
  };

  class const_iterator
  {
  public:
    const_iterator(const FilteredMatrix<VectorType> *matrix,
                   const size_type                   index);

    const_iterator &
    operator++();

    const_iterator &
    operator++(int);

    const Accessor &operator*() const;

    const Accessor *operator->() const;

    bool
    operator==(const const_iterator &) const;
    bool
    operator!=(const const_iterator &) const;

    bool
    operator<(const const_iterator &) const;

    bool
    operator>(const const_iterator &) const;

  private:
    Accessor accessor;
  };

  using IndexValuePair = std::pair<size_type, double>;

  FilteredMatrix();

  FilteredMatrix(const FilteredMatrix &fm);

  template <typename MatrixType>
  FilteredMatrix(const MatrixType &matrix,
                 const bool        expect_constrained_source = false);

  FilteredMatrix &
  operator=(const FilteredMatrix &fm);

  template <typename MatrixType>
  void
  initialize(const MatrixType &m, const bool expect_constrained_source = false);

  void
  clear();

  void
  add_constraint(const size_type i, const double v);

  template <class ConstraintList>
  void
  add_constraints(const ConstraintList &new_constraints);

  void
  clear_constraints();

  DEAL_II_DEPRECATED
  void
  apply_constraints(VectorType &v, const bool matrix_is_symmetric) const;
  void
  apply_constraints(VectorType &v) const;

  void
  vmult(VectorType &dst, const VectorType &src) const;

  void
  Tvmult(VectorType &dst, const VectorType &src) const;

  void
  vmult_add(VectorType &dst, const VectorType &src) const;

  void
  Tvmult_add(VectorType &dst, const VectorType &src) const;

  const_iterator
  begin() const;
  const_iterator
  end() const;

  std::size_t
  memory_consumption() const;

private:
  bool expect_constrained_source;

  using const_index_value_iterator =
    typename std::vector<IndexValuePair>::const_iterator;

  struct PairComparison
  {
    bool
    operator()(const IndexValuePair &i1, const IndexValuePair &i2) const;
  };

  std::shared_ptr<PointerMatrixBase<VectorType>> matrix;

  std::vector<IndexValuePair> constraints;

  void
  pre_filter(VectorType &v) const;

  void
  post_filter(const VectorType &in, VectorType &out) const;

  friend class Accessor;
  friend class FilteredMatrixBlock<VectorType>;
};

/*---------------------- Inline functions -----------------------------------*/


//--------------------------------Iterators--------------------------------------//

template <typename VectorType>
inline FilteredMatrix<VectorType>::Accessor::Accessor(
  const FilteredMatrix<VectorType> *matrix,
  const size_type                   index)
  : matrix(matrix)
  , index(index)
{
  Assert(index <= matrix->constraints.size(),
         ExcIndexRange(index, 0, matrix->constraints.size()));
}



template <typename VectorType>
inline types::global_dof_index
FilteredMatrix<VectorType>::Accessor::row() const
{
  return matrix->constraints[index].first;
}



template <typename VectorType>
inline types::global_dof_index
FilteredMatrix<VectorType>::Accessor::column() const
{
  return matrix->constraints[index].first;
}



template <typename VectorType>
inline double
FilteredMatrix<VectorType>::Accessor::value() const
{
  return matrix->constraints[index].second;
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::Accessor::advance()
{
  Assert(index < matrix->constraints.size(), ExcIteratorPastEnd());
  ++index;
}



template <typename VectorType>
inline FilteredMatrix<VectorType>::const_iterator::const_iterator(
  const FilteredMatrix<VectorType> *matrix,
  const size_type                   index)
  : accessor(matrix, index)
{}



template <typename VectorType>
inline typename FilteredMatrix<VectorType>::const_iterator &
FilteredMatrix<VectorType>::const_iterator::operator++()
{
  accessor.advance();
  return *this;
}


template <typename number>
inline const typename FilteredMatrix<number>::Accessor &
  FilteredMatrix<number>::const_iterator::operator*() const
{
  return accessor;
}


template <typename number>
inline const typename FilteredMatrix<number>::Accessor *
  FilteredMatrix<number>::const_iterator::operator->() const
{
  return &accessor;
}


template <typename number>
inline bool
FilteredMatrix<number>::const_iterator::
operator==(const const_iterator &other) const
{
  return (accessor.index == other.accessor.index &&
          accessor.matrix == other.accessor.matrix);
}


template <typename number>
inline bool
FilteredMatrix<number>::const_iterator::
operator!=(const const_iterator &other) const
{
  return !(*this == other);
}



//------------------------- FilteredMatrix ----------------------------------//

template <typename number>
inline typename FilteredMatrix<number>::const_iterator
FilteredMatrix<number>::begin() const
{
  return const_iterator(this, 0);
}


template <typename number>
inline typename FilteredMatrix<number>::const_iterator
FilteredMatrix<number>::end() const
{
  return const_iterator(this, constraints.size());
}


template <typename VectorType>
inline bool
FilteredMatrix<VectorType>::PairComparison::
operator()(const IndexValuePair &i1, const IndexValuePair &i2) const
{
  return (i1.first < i2.first);
}



template <typename VectorType>
template <typename MatrixType>
inline void
FilteredMatrix<VectorType>::initialize(const MatrixType &m, bool ecs)
{
  matrix.reset(new_pointer_matrix_base(m, VectorType()));

  expect_constrained_source = ecs;
}



template <typename VectorType>
inline FilteredMatrix<VectorType>::FilteredMatrix()
  : expect_constrained_source(false)
{}



template <typename VectorType>
inline FilteredMatrix<VectorType>::FilteredMatrix(const FilteredMatrix &fm)
  : Subscriptor()
  , expect_constrained_source(fm.expect_constrained_source)
  , matrix(fm.matrix)
  , constraints(fm.constraints)
{}



template <typename VectorType>
template <typename MatrixType>
inline FilteredMatrix<VectorType>::FilteredMatrix(const MatrixType &m,
                                                  const bool        ecs)
  : expect_constrained_source(false)
{
  initialize(m, ecs);
}



template <typename VectorType>
inline FilteredMatrix<VectorType> &
FilteredMatrix<VectorType>::operator=(const FilteredMatrix &fm)
{
  matrix                    = fm.matrix;
  expect_constrained_source = fm.expect_constrained_source;
  constraints               = fm.constraints;
  return *this;
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::add_constraint(const size_type index,
                                           const double    value)
{
  // add new constraint to end
  constraints.push_back(IndexValuePair(index, value));
}



template <typename VectorType>
template <class ConstraintList>
inline void
FilteredMatrix<VectorType>::add_constraints(
  const ConstraintList &new_constraints)
{
  // add new constraints to end
  const size_type old_size = constraints.size();
  constraints.reserve(old_size + new_constraints.size());
  constraints.insert(constraints.end(),
                     new_constraints.begin(),
                     new_constraints.end());
  // then merge the two arrays to
  // form one sorted one
  std::inplace_merge(constraints.begin(),
                     constraints.begin() + old_size,
                     constraints.end(),
                     PairComparison());
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::clear_constraints()
{
  // swap vectors to release memory
  std::vector<IndexValuePair> empty;
  constraints.swap(empty);
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::clear()
{
  clear_constraints();
  matrix.reset();
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::apply_constraints(
  VectorType &v,
  const bool /* matrix_is_symmetric */) const
{
  apply_constraints(v);
}


template <typename VectorType>
inline void
FilteredMatrix<VectorType>::apply_constraints(VectorType &v) const
{
  GrowingVectorMemory<VectorType>            mem;
  typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
  tmp_vector->reinit(v);
  const_index_value_iterator       i = constraints.begin();
  const const_index_value_iterator e = constraints.end();
  for (; i != e; ++i)
    {
      AssertIsFinite(i->second);
      (*tmp_vector)(i->first) = -i->second;
    }

  // This vmult is without bc, to get
  // the rhs correction in a correct
  // way.
  matrix->vmult_add(v, *tmp_vector);
  // finally set constrained
  // entries themselves
  for (i = constraints.begin(); i != e; ++i)
    {
      AssertIsFinite(i->second);
      v(i->first) = i->second;
    }
}


template <typename VectorType>
inline void
FilteredMatrix<VectorType>::pre_filter(VectorType &v) const
{
  // iterate over all constraints and
  // zero out value
  const_index_value_iterator       i = constraints.begin();
  const const_index_value_iterator e = constraints.end();
  for (; i != e; ++i)
    v(i->first) = 0;
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::post_filter(const VectorType &in,
                                        VectorType &      out) const
{
  // iterate over all constraints and
  // set value correctly
  const_index_value_iterator       i = constraints.begin();
  const const_index_value_iterator e = constraints.end();
  for (; i != e; ++i)
    {
      AssertIsFinite(in(i->first));
      out(i->first) = in(i->first);
    }
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
{
  if (!expect_constrained_source)
    {
      GrowingVectorMemory<VectorType>            mem;
      typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
      // first copy over src vector and
      // pre-filter
      tmp_vector->reinit(src, true);
      *tmp_vector = src;
      pre_filter(*tmp_vector);
      // then let matrix do its work
      matrix->vmult(dst, *tmp_vector);
    }
  else
    {
      matrix->vmult(dst, src);
    }

  // finally do post-filtering
  post_filter(src, dst);
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
{
  if (!expect_constrained_source)
    {
      GrowingVectorMemory<VectorType>            mem;
      typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
      // first copy over src vector and
      // pre-filter
      tmp_vector->reinit(src, true);
      *tmp_vector = src;
      pre_filter(*tmp_vector);
      // then let matrix do its work
      matrix->Tvmult(dst, *tmp_vector);
    }
  else
    {
      matrix->Tvmult(dst, src);
    }

  // finally do post-filtering
  post_filter(src, dst);
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::vmult_add(VectorType &      dst,
                                      const VectorType &src) const
{
  if (!expect_constrained_source)
    {
      GrowingVectorMemory<VectorType>            mem;
      typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
      // first copy over src vector and
      // pre-filter
      tmp_vector->reinit(src, true);
      *tmp_vector = src;
      pre_filter(*tmp_vector);
      // then let matrix do its work
      matrix->vmult_add(dst, *tmp_vector);
    }
  else
    {
      matrix->vmult_add(dst, src);
    }

  // finally do post-filtering
  post_filter(src, dst);
}



template <typename VectorType>
inline void
FilteredMatrix<VectorType>::Tvmult_add(VectorType &      dst,
                                       const VectorType &src) const
{
  if (!expect_constrained_source)
    {
      GrowingVectorMemory<VectorType>            mem;
      typename VectorMemory<VectorType>::Pointer tmp_vector(mem);
      // first copy over src vector and
      // pre-filter
      tmp_vector->reinit(src, true);
      *tmp_vector = src;
      pre_filter(*tmp_vector);
      // then let matrix do its work
      matrix->Tvmult_add(dst, *tmp_vector);
    }
  else
    {
      matrix->Tvmult_add(dst, src);
    }

  // finally do post-filtering
  post_filter(src, dst);
}



template <typename VectorType>
inline std::size_t
FilteredMatrix<VectorType>::memory_consumption() const
{
  return (MemoryConsumption::memory_consumption(matrix) +
          MemoryConsumption::memory_consumption(constraints));
}



DEAL_II_NAMESPACE_CLOSE

#endif
/*----------------------------   filtered_matrix.h ---------------------------*/