precondition_block.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 2017 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_precondition_block_h
#define dealii_precondition_block_h


#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/lac/precondition_block_base.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN

template <typename MatrixType, typename inverse_type = typename MatrixType::value_type>
class PreconditionBlock
  : public virtual Subscriptor,
    protected PreconditionBlockBase<inverse_type>
{
private:
  typedef typename MatrixType::value_type number;

  typedef inverse_type value_type;

public:
  typedef types::global_dof_index size_type;

  class AdditionalData
  {
  public:
    AdditionalData (const size_type block_size,
                    const double relaxation = 1.,
                    const bool invert_diagonal = true,
                    const bool same_diagonal = false);

    double relaxation;

    size_type block_size;

    bool invert_diagonal;

    bool same_diagonal;
    typename PreconditionBlockBase<inverse_type>::Inversion inversion;

    double threshold;
  };


  PreconditionBlock(bool store_diagonals = false);

  ~PreconditionBlock() = default;

  void initialize (const MatrixType &A,
                   const AdditionalData parameters);
protected:
  void initialize (const MatrixType &A,
                   const std::vector<size_type> &permutation,
                   const std::vector<size_type> &inverse_permutation,
                   const AdditionalData parameters);

  void set_permutation(const std::vector<size_type> &permutation,
                       const std::vector<size_type> &inverse_permutation);

  void invert_permuted_diagblocks();
public:
  void clear();

  bool empty () const;

  value_type el(size_type i,
                size_type j) const;

  void invert_diagblocks();

  template <typename number2>
  void forward_step (Vector<number2>       &dst,
                     const Vector<number2> &prev,
                     const Vector<number2> &src,
                     const bool             transpose_diagonal) const;

  template <typename number2>
  void backward_step (Vector<number2>       &dst,
                      const Vector<number2> &prev,
                      const Vector<number2> &src,
                      const bool             transpose_diagonal) const;


  size_type block_size () const;

  std::size_t memory_consumption () const;

  DeclException2 (ExcWrongBlockSize,
                  int, int,
                  << "The blocksize " << arg1
                  << " and the size of the matrix " << arg2
                  << " do not match.");

  DeclException0 (ExcInverseMatricesAlreadyExist);


protected:
  size_type blocksize;

  SmartPointer<const MatrixType,PreconditionBlock<MatrixType,inverse_type> > A;
  double relaxation;

  std::vector<size_type> permutation;

  std::vector<size_type> inverse_permutation;
};



template <typename MatrixType, typename inverse_type = typename MatrixType::value_type>
class PreconditionBlockJacobi : public virtual Subscriptor,
  private PreconditionBlock<MatrixType, inverse_type>
{
private:
  typedef typename MatrixType::value_type number;

public:
  typedef types::global_dof_index size_type;

  class const_iterator
  {
  private:
    class Accessor
    {
    public:
      Accessor (const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix,
                const size_type row);

      size_type row() const;

      size_type column() const;

      inverse_type value() const;

    protected:
      const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix;

      size_type bs;

      size_type a_block;

      typename FullMatrix<inverse_type>::const_iterator b_iterator;

      typename FullMatrix<inverse_type>::const_iterator b_end;

      friend class const_iterator;
    };

  public:
    const_iterator(const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix,
                   const size_type row);

    const_iterator &operator++ ();

    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;

  private:
    Accessor accessor;
  };

  using typename PreconditionBlock<MatrixType,inverse_type>::AdditionalData;
  using PreconditionBlock<MatrixType, inverse_type>::initialize;
  using PreconditionBlock<MatrixType, inverse_type>::clear;
  using PreconditionBlock<MatrixType, inverse_type>::empty;
  using PreconditionBlock<MatrixType, inverse_type>::el;
  using PreconditionBlock<MatrixType, inverse_type>::invert_diagblocks;
  using PreconditionBlock<MatrixType, inverse_type>::block_size;
  using PreconditionBlockBase<inverse_type>::size;
  using PreconditionBlockBase<inverse_type>::inverse;
  using PreconditionBlockBase<inverse_type>::inverse_householder;
  using PreconditionBlockBase<inverse_type>::inverse_svd;
  using PreconditionBlockBase<inverse_type>::log_statistics;
  using PreconditionBlock<MatrixType, inverse_type>::set_permutation;

  template <typename number2>
  void vmult (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void Tvmult (Vector<number2> &, const Vector<number2> &) const;
  template <typename number2>
  void vmult_add (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void Tvmult_add (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void step (Vector<number2> &dst, const Vector<number2> &rhs) const;

  template <typename number2>
  void Tstep (Vector<number2> &dst, const Vector<number2> &rhs) const;

  const_iterator begin () const;

  const_iterator end () const;

  const_iterator begin (const size_type r) const;

  const_iterator end (const size_type r) const;


private:
  template <typename number2>
  void do_vmult (Vector<number2> &,
                 const Vector<number2> &,
                 bool adding) const;

  friend class Accessor;
  friend class const_iterator;
};



template <typename MatrixType, typename inverse_type = typename MatrixType::value_type>
class PreconditionBlockSOR : public virtual Subscriptor,
  protected PreconditionBlock<MatrixType, inverse_type>
{
public:
  typedef types::global_dof_index size_type;

  PreconditionBlockSOR();

  typedef typename MatrixType::value_type number;

  using typename PreconditionBlock<MatrixType,inverse_type>::AdditionalData;
  using PreconditionBlock<MatrixType, inverse_type>::initialize;
  using PreconditionBlock<MatrixType, inverse_type>::clear;
  using PreconditionBlockBase<inverse_type>::size;
  using PreconditionBlockBase<inverse_type>::inverse;
  using PreconditionBlockBase<inverse_type>::inverse_householder;
  using PreconditionBlockBase<inverse_type>::inverse_svd;
  using PreconditionBlock<MatrixType, inverse_type>::invert_diagblocks;
  using PreconditionBlock<MatrixType, inverse_type>::set_permutation;
  using PreconditionBlockBase<inverse_type>::log_statistics;

  template <typename number2>
  void vmult (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void vmult_add (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void Tvmult (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void Tvmult_add (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void step (Vector<number2> &dst, const Vector<number2> &rhs) const;

  template <typename number2>
  void Tstep (Vector<number2> &dst, const Vector<number2> &rhs) const;

protected:
  PreconditionBlockSOR(bool store);

  template <typename number2>
  void forward (Vector<number2> &,
                const Vector<number2> &,
                const bool transpose_diagonal,
                const bool adding) const;

  template <typename number2>
  void backward (Vector<number2> &,
                 const Vector<number2> &,
                 const bool transpose_diagonal,
                 const bool adding) const;
};


template <typename MatrixType, typename inverse_type = typename MatrixType::value_type>
class PreconditionBlockSSOR : public virtual Subscriptor,
  private PreconditionBlockSOR<MatrixType, inverse_type>
{
public:
  typedef types::global_dof_index size_type;

  typedef typename MatrixType::value_type number;

  PreconditionBlockSSOR ();

  // Keep AdditionalData accessible
  using typename PreconditionBlockSOR<MatrixType,inverse_type>::AdditionalData;

  // The following are the
  // functions of the base classes
  // which we want to keep
  // accessible.
  using PreconditionBlockSOR<MatrixType,inverse_type>::initialize;
  using PreconditionBlockSOR<MatrixType,inverse_type>::clear;
  using PreconditionBlockBase<inverse_type>::size;
  using PreconditionBlockBase<inverse_type>::inverse;
  using PreconditionBlockBase<inverse_type>::inverse_householder;
  using PreconditionBlockBase<inverse_type>::inverse_svd;
  using PreconditionBlockBase<inverse_type>::log_statistics;
  using PreconditionBlockSOR<MatrixType,inverse_type>::set_permutation;
  using PreconditionBlockSOR<MatrixType, inverse_type>::empty;
  using PreconditionBlockSOR<MatrixType, inverse_type>::el;
  using PreconditionBlockSOR<MatrixType,inverse_type>::invert_diagblocks;

  template <typename number2>
  void vmult (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void Tvmult (Vector<number2> &, const Vector<number2> &) const;

  template <typename number2>
  void step (Vector<number2> &dst, const Vector<number2> &rhs) const;

  template <typename number2>
  void Tstep (Vector<number2> &dst, const Vector<number2> &rhs) const;
};

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

#ifndef DOXYGEN

template <typename MatrixType, typename inverse_type>
inline bool
PreconditionBlock<MatrixType, inverse_type>::empty () const
{
  if (A == nullptr)
    return true;
  return A->empty();
}


template <typename MatrixType, typename inverse_type>
inline inverse_type
PreconditionBlock<MatrixType, inverse_type>::el (
  size_type i,
  size_type j) const
{
  const size_type bs = blocksize;
  const unsigned int nb = i/bs;

  const FullMatrix<inverse_type> &B = this->inverse(nb);

  const size_type ib = i % bs;
  const size_type jb = j % bs;

  if (jb + nb*bs != j)
    {
      return 0.;
    }

  return B(ib, jb);
}

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

template <typename MatrixType, typename inverse_type>
inline
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::Accessor::
Accessor (const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix,
          const size_type row)
  :
  matrix(matrix),
  b_iterator(&matrix->inverse(0), 0, 0),
  b_end(&matrix->inverse(0), 0, 0)
{
  bs = matrix->block_size();
  a_block = row / bs;

  // This is the end accessor, which
  // does not have a valid block.
  if (a_block == matrix->size())
    return;

  const size_type r = row % bs;

  b_iterator = matrix->inverse(a_block).begin(r);
  b_end = matrix->inverse(a_block).end();

  Assert (a_block < matrix->size(),
          ExcIndexRange(a_block, 0, matrix->size()));
}


template <typename MatrixType, typename inverse_type>
inline
typename PreconditionBlockJacobi<MatrixType, inverse_type>::size_type
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::Accessor::row() const
{
  Assert (a_block < matrix->size(),
          ExcIteratorPastEnd());

  return bs * a_block + b_iterator->row();
}


template <typename MatrixType, typename inverse_type>
inline
typename PreconditionBlockJacobi<MatrixType, inverse_type>::size_type
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::Accessor::column() const
{
  Assert (a_block < matrix->size(),
          ExcIteratorPastEnd());

  return bs * a_block + b_iterator->column();
}


template <typename MatrixType, typename inverse_type>
inline
inverse_type
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::Accessor::value() const
{
  Assert (a_block < matrix->size(),
          ExcIteratorPastEnd());

  return b_iterator->value();
}


template <typename MatrixType, typename inverse_type>
inline
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
const_iterator(const PreconditionBlockJacobi<MatrixType, inverse_type> *matrix,
               const size_type row)
  :
  accessor(matrix, row)
{}


template <typename MatrixType, typename inverse_type>
inline
typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator &
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::operator++ ()
{
  Assert (*this != accessor.matrix->end(), ExcIteratorPastEnd());

  ++accessor.b_iterator;
  if (accessor.b_iterator == accessor.b_end)
    {
      ++accessor.a_block;

      if (accessor.a_block < accessor.matrix->size())
        {
          accessor.b_iterator = accessor.matrix->inverse(accessor.a_block).begin();
          accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
        }
    }
  return *this;
}


template <typename MatrixType, typename inverse_type>
inline
const typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::Accessor &
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::operator* () const
{
  return accessor;
}


template <typename MatrixType, typename inverse_type>
inline
const typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::Accessor *
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::operator-> () const
{
  return &accessor;
}


template <typename MatrixType, typename inverse_type>
inline
bool
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
operator == (const const_iterator &other) const
{
  if (accessor.a_block == accessor.matrix->size() &&
      accessor.a_block == other.accessor.a_block)
    return true;

  if (accessor.a_block != other.accessor.a_block)
    return false;

  return (accessor.row() == other.accessor.row() &&
          accessor.column() == other.accessor.column());
}


template <typename MatrixType, typename inverse_type>
inline
bool
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
operator != (const const_iterator &other) const
{
  return ! (*this == other);
}


template <typename MatrixType, typename inverse_type>
inline
bool
PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator::
operator < (const const_iterator &other) const
{
  return (accessor.row() < other.accessor.row() ||
          (accessor.row() == other.accessor.row() &&
           accessor.column() < other.accessor.column()));
}


template <typename MatrixType, typename inverse_type>
inline
typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator
PreconditionBlockJacobi<MatrixType, inverse_type>::begin () const
{
  return const_iterator(this, 0);
}


template <typename MatrixType, typename inverse_type>
inline
typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator
PreconditionBlockJacobi<MatrixType, inverse_type>::end () const
{
  return const_iterator(this, this->size() * this->block_size());
}


template <typename MatrixType, typename inverse_type>
inline
typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator
PreconditionBlockJacobi<MatrixType, inverse_type>::begin (
  const size_type r) const
{
  Assert (r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
  return const_iterator(this, r);
}



template <typename MatrixType, typename inverse_type>
inline
typename PreconditionBlockJacobi<MatrixType, inverse_type>::const_iterator
PreconditionBlockJacobi<MatrixType, inverse_type>::end (
  const size_type r) const
{
  Assert (r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
  return const_iterator(this, r+1);
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif