precondition_block.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 1999 - 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.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;
 
      // Make enclosing class a friend.
      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)
  , bs(matrix->block_size())
  , a_block(row / bs)
  , b_iterator(&matrix->inverse(0), 0, 0)
  , b_end(&matrix->inverse(0), 0, 0)
{
  // 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();
 
  AssertIndexRange(a_block, 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
{
  AssertIndexRange(r, 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
{
  AssertIndexRange(r, this->A->m());
  return const_iterator(this, r + 1);
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif