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

#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/smartpointer.h>
#include <deal.II/base/subscriptor.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:
  using number = typename MatrixType::value_type;

  using value_type = inverse_type;

public:
  using size_type = types::global_dof_index;

  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() override = 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:
  using number = typename MatrixType::value_type;

public:
  using size_type = types::global_dof_index;

  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:
  using size_type = types::global_dof_index;

  PreconditionBlockSOR();

  using number = typename MatrixType::value_type;

  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:
  using size_type = types::global_dof_index;

  using number = typename MatrixType::value_type;

  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