filtered_matrix.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2016 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__filtered_matrix_h
#define dealii__filtered_matrix_h



#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/lac/pointer_matrix.h>
#include <deal.II/lac/vector_memory.h>
#include <vector>
#include <algorithm>

DEAL_II_NAMESPACE_OPEN

template <class VectorType> class FilteredMatrixBlock;

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

  typedef types::global_dof_index size_type;

  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;
  };

  typedef std::pair<size_type, double> IndexValuePair;

  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 ();

  void apply_constraints (VectorType &v,
                          const bool matrix_is_symmetric) const DEAL_II_DEPRECATED;
  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;

  typedef typename std::vector<IndexValuePair>::const_iterator const_index_value_iterator;

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

  std_cxx11::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;
      VectorType *tmp_vector = mem.alloc();
      // 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);
      mem.free(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;
      VectorType *tmp_vector = mem.alloc();
      // 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);
      mem.free(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;
      VectorType *tmp_vector = mem.alloc();
      // 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);
      mem.free(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;
      VectorType *tmp_vector = mem.alloc();
      // 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);
      mem.free(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     ---------------------------*/