precondition.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1999 - 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_precondition_h
#define dealii_precondition_h

// This file contains simple preconditioners.

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

#include <deal.II/base/parallel.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>

#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/vector_memory.h>

DEAL_II_NAMESPACE_OPEN

// forward declarations

template <typename number>
class Vector;
template <typename number>
class SparseMatrix;
namespace LinearAlgebra
{
  namespace distributed
  {
    template <typename number>
    class Vector;
  }
} // namespace LinearAlgebra



class PreconditionIdentity : public Subscriptor
{
public:
  using size_type = types::global_dof_index;

  struct AdditionalData
  {
    AdditionalData() = default;
  };

  PreconditionIdentity();

  template <typename MatrixType>
  void
  initialize(const MatrixType &    matrix,
             const AdditionalData &additional_data = AdditionalData());

  template <class VectorType>
  void
  vmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  Tvmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  vmult_add(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  Tvmult_add(VectorType &, const VectorType &) const;

  void
  clear()
  {}

  size_type
  m() const;

  size_type
  n() const;

private:
  size_type n_rows;

  size_type n_columns;
};



class PreconditionRichardson : public Subscriptor
{
public:
  using size_type = types::global_dof_index;

  class AdditionalData
  {
  public:
    AdditionalData(const double relaxation = 1.);

    double relaxation;
  };

  PreconditionRichardson();

  void
  initialize(const AdditionalData &parameters);

  template <typename MatrixType>
  void
  initialize(const MatrixType &matrix, const AdditionalData &parameters);

  template <class VectorType>
  void
  vmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  Tvmult(VectorType &, const VectorType &) const;
  template <class VectorType>
  void
  vmult_add(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  Tvmult_add(VectorType &, const VectorType &) const;

  void
  clear()
  {}

  size_type
  m() const;

  size_type
  n() const;

private:
  double relaxation;

  size_type n_rows;

  size_type n_columns;
};



template <typename MatrixType = SparseMatrix<double>,
          class VectorType    = Vector<double>>
class PreconditionUseMatrix : public Subscriptor
{
public:
  using function_ptr = void (MatrixType::*)(VectorType &,
                                            const VectorType &) const;

  PreconditionUseMatrix(const MatrixType &M, const function_ptr method);

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

private:
  const MatrixType &matrix;

  const function_ptr precondition;
};



template <typename MatrixType = SparseMatrix<double>>
class PreconditionRelaxation : public Subscriptor
{
public:
  using size_type = typename MatrixType::size_type;

  class AdditionalData
  {
  public:
    AdditionalData(const double relaxation = 1.);

    double relaxation;
  };

  void
  initialize(const MatrixType &    A,
             const AdditionalData &parameters = AdditionalData());

  void
  clear();

  size_type
  m() const;

  size_type
  n() const;

protected:
  SmartPointer<const MatrixType, PreconditionRelaxation<MatrixType>> A;

  double relaxation;
};



template <typename MatrixType = SparseMatrix<double>>
class PreconditionJacobi : public PreconditionRelaxation<MatrixType>
{
public:
  using AdditionalData =
    typename PreconditionRelaxation<MatrixType>::AdditionalData;

  template <class VectorType>
  void
  vmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  Tvmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  step(VectorType &x, const VectorType &rhs) const;

  template <class VectorType>
  void
  Tstep(VectorType &x, const VectorType &rhs) const;
};


template <typename MatrixType = SparseMatrix<double>>
class PreconditionSOR : public PreconditionRelaxation<MatrixType>
{
public:
  using AdditionalData =
    typename PreconditionRelaxation<MatrixType>::AdditionalData;

  template <class VectorType>
  void
  vmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  Tvmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  step(VectorType &x, const VectorType &rhs) const;

  template <class VectorType>
  void
  Tstep(VectorType &x, const VectorType &rhs) const;
};



template <typename MatrixType = SparseMatrix<double>>
class PreconditionSSOR : public PreconditionRelaxation<MatrixType>
{
public:
  using AdditionalData =
    typename PreconditionRelaxation<MatrixType>::AdditionalData;

  using size_type = typename MatrixType::size_type;

  using BaseClass = PreconditionRelaxation<MatrixType>;


  void
  initialize(const MatrixType &                        A,
             const typename BaseClass::AdditionalData &parameters =
               typename BaseClass::AdditionalData());

  template <class VectorType>
  void
  vmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  Tvmult(VectorType &, const VectorType &) const;


  template <class VectorType>
  void
  step(VectorType &x, const VectorType &rhs) const;

  template <class VectorType>
  void
  Tstep(VectorType &x, const VectorType &rhs) const;

private:
  std::vector<std::size_t> pos_right_of_diagonal;
};


template <typename MatrixType = SparseMatrix<double>>
class PreconditionPSOR : public PreconditionRelaxation<MatrixType>
{
public:
  using size_type = typename MatrixType::size_type;

  class AdditionalData
  {
  public:
    AdditionalData(
      const std::vector<size_type> &permutation,
      const std::vector<size_type> &inverse_permutation,
      const typename PreconditionRelaxation<MatrixType>::AdditionalData
        &parameters =
          typename PreconditionRelaxation<MatrixType>::AdditionalData());

    const std::vector<size_type> &permutation;
    const std::vector<size_type> &inverse_permutation;
    typename PreconditionRelaxation<MatrixType>::AdditionalData parameters;
  };

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

  void
  initialize(const MatrixType &A, const AdditionalData &additional_data);

  template <class VectorType>
  void
  vmult(VectorType &, const VectorType &) const;

  template <class VectorType>
  void
  Tvmult(VectorType &, const VectorType &) const;

private:
  const std::vector<size_type> *permutation;
  const std::vector<size_type> *inverse_permutation;
};



template <typename MatrixType         = SparseMatrix<double>,
          typename VectorType         = Vector<double>,
          typename PreconditionerType = DiagonalMatrix<VectorType>>
class PreconditionChebyshev : public Subscriptor
{
public:
  using size_type = types::global_dof_index;

  // avoid warning about use of deprecated variables

  struct AdditionalData
  {
    AdditionalData(const unsigned int degree              = 0,
                   const double       smoothing_range     = 0.,
                   const bool         nonzero_starting    = false,
                   const unsigned int eig_cg_n_iterations = 8,
                   const double       eig_cg_residual     = 1e-2,
                   const double       max_eigenvalue      = 1);

    unsigned int degree;

    double smoothing_range;

    bool nonzero_starting DEAL_II_DEPRECATED;

    unsigned int eig_cg_n_iterations;

    double eig_cg_residual;

    double max_eigenvalue;

    VectorType matrix_diagonal_inverse DEAL_II_DEPRECATED;

    std::shared_ptr<PreconditionerType> preconditioner;
  };


  PreconditionChebyshev();

  void
  initialize(const MatrixType &    matrix,
             const AdditionalData &additional_data = AdditionalData());

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

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

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

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

  void
  clear();

  size_type
  m() const;

  size_type
  n() const;

private:
  SmartPointer<
    const MatrixType,
    PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>>
    matrix_ptr;

  mutable VectorType update1;

  mutable VectorType update2;

  mutable VectorType update3;

  AdditionalData data;

  double theta;

  double delta;

  bool eigenvalues_are_initialized;

  mutable Threads::Mutex mutex;

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

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

  void
  estimate_eigenvalues(const VectorType &src) const;
};



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

#ifndef DOXYGEN

inline PreconditionIdentity::PreconditionIdentity()
  : n_rows(0)
  , n_columns(0)
{}

template <typename MatrixType>
inline void
PreconditionIdentity::initialize(const MatrixType &matrix,
                                 const PreconditionIdentity::AdditionalData &)
{
  n_rows    = matrix.m();
  n_columns = matrix.n();
}


template <class VectorType>
inline void
PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
{
  dst = src;
}



template <class VectorType>
inline void
PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
{
  dst = src;
}

template <class VectorType>
inline void
PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
{
  dst += src;
}



template <class VectorType>
inline void
PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
{
  dst += src;
}

inline PreconditionIdentity::size_type
PreconditionIdentity::m() const
{
  Assert(n_rows != 0, ExcNotInitialized());
  return n_rows;
}

inline PreconditionIdentity::size_type
PreconditionIdentity::n() const
{
  Assert(n_columns != 0, ExcNotInitialized());
  return n_columns;
}

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

inline PreconditionRichardson::AdditionalData::AdditionalData(
  const double relaxation)
  : relaxation(relaxation)
{}


inline PreconditionRichardson::PreconditionRichardson()
  : relaxation(0)
  , n_rows(0)
  , n_columns(0)
{
  AdditionalData add_data;
  relaxation = add_data.relaxation;
}



inline void
PreconditionRichardson::initialize(
  const PreconditionRichardson::AdditionalData &parameters)
{
  relaxation = parameters.relaxation;
}



template <typename MatrixType>
inline void
PreconditionRichardson::initialize(
  const MatrixType &                            matrix,
  const PreconditionRichardson::AdditionalData &parameters)
{
  relaxation = parameters.relaxation;
  n_rows     = matrix.m();
  n_columns  = matrix.n();
}



template <class VectorType>
inline void
PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
{
  static_assert(
    std::is_same<size_type, typename VectorType::size_type>::value,
    "PreconditionRichardson and VectorType must have the same size_type.");

  dst.equ(relaxation, src);
}



template <class VectorType>
inline void
PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
{
  static_assert(
    std::is_same<size_type, typename VectorType::size_type>::value,
    "PreconditionRichardson and VectorType must have the same size_type.");

  dst.equ(relaxation, src);
}

template <class VectorType>
inline void
PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
{
  static_assert(
    std::is_same<size_type, typename VectorType::size_type>::value,
    "PreconditionRichardson and VectorType must have the same size_type.");

  dst.add(relaxation, src);
}



template <class VectorType>
inline void
PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
{
  static_assert(
    std::is_same<size_type, typename VectorType::size_type>::value,
    "PreconditionRichardson and VectorType must have the same size_type.");

  dst.add(relaxation, src);
}

inline PreconditionRichardson::size_type
PreconditionRichardson::m() const
{
  Assert(n_rows != 0, ExcNotInitialized());
  return n_rows;
}

inline PreconditionRichardson::size_type
PreconditionRichardson::n() const
{
  Assert(n_columns != 0, ExcNotInitialized());
  return n_columns;
}

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

template <typename MatrixType>
inline void
PreconditionRelaxation<MatrixType>::initialize(const MatrixType &    rA,
                                               const AdditionalData &parameters)
{
  A          = &rA;
  relaxation = parameters.relaxation;
}


template <typename MatrixType>
inline void
PreconditionRelaxation<MatrixType>::clear()
{
  A = nullptr;
}

template <typename MatrixType>
inline typename PreconditionRelaxation<MatrixType>::size_type
PreconditionRelaxation<MatrixType>::m() const
{
  Assert(A != nullptr, ExcNotInitialized());
  return A->m();
}

template <typename MatrixType>
inline typename PreconditionRelaxation<MatrixType>::size_type
PreconditionRelaxation<MatrixType>::n() const
{
  Assert(A != nullptr, ExcNotInitialized());
  return A->n();
}

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

template <typename MatrixType>
template <class VectorType>
inline void
PreconditionJacobi<MatrixType>::vmult(VectorType &      dst,
                                      const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionJacobi and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->precondition_Jacobi(dst, src, this->relaxation);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionJacobi<MatrixType>::Tvmult(VectorType &      dst,
                                       const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionJacobi and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->precondition_Jacobi(dst, src, this->relaxation);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionJacobi<MatrixType>::step(VectorType &      dst,
                                     const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionJacobi and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->Jacobi_step(dst, src, this->relaxation);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionJacobi<MatrixType>::Tstep(VectorType &      dst,
                                      const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionJacobi<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionJacobi and VectorType must have the same size_type.");

  step(dst, src);
}



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

template <typename MatrixType>
template <class VectorType>
inline void
PreconditionSOR<MatrixType>::vmult(VectorType &dst, const VectorType &src) const
{
  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
                             typename VectorType::size_type>::value,
                "PreconditionSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->precondition_SOR(dst, src, this->relaxation);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionSOR<MatrixType>::Tvmult(VectorType &      dst,
                                    const VectorType &src) const
{
  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
                             typename VectorType::size_type>::value,
                "PreconditionSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->precondition_TSOR(dst, src, this->relaxation);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionSOR<MatrixType>::step(VectorType &dst, const VectorType &src) const
{
  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
                             typename VectorType::size_type>::value,
                "PreconditionSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->SOR_step(dst, src, this->relaxation);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionSOR<MatrixType>::Tstep(VectorType &dst, const VectorType &src) const
{
  static_assert(std::is_same<typename PreconditionSOR<MatrixType>::size_type,
                             typename VectorType::size_type>::value,
                "PreconditionSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->TSOR_step(dst, src, this->relaxation);
}



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

template <typename MatrixType>
inline void
PreconditionSSOR<MatrixType>::initialize(
  const MatrixType &                        rA,
  const typename BaseClass::AdditionalData &parameters)
{
  this->PreconditionRelaxation<MatrixType>::initialize(rA, parameters);

  // in case we have a SparseMatrix class, we can extract information about
  // the diagonal.
  const SparseMatrix<typename MatrixType::value_type> *mat =
    dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
      &*this->A);

  // calculate the positions first after the diagonal.
  if (mat != nullptr)
    {
      const size_type n = this->A->n();
      pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
      for (size_type row = 0; row < n; ++row)
        {
          // find the first element in this line which is on the right of the
          // diagonal.  we need to precondition with the elements on the left
          // only. note: the first entry in each line denotes the diagonal
          // element, which we need not check.
          typename SparseMatrix<typename MatrixType::value_type>::const_iterator
            it = mat->begin(row) + 1;
          for (; it < mat->end(row); ++it)
            if (it->column() > row)
              break;
          pos_right_of_diagonal[row] = it - mat->begin();
        }
    }
}


template <typename MatrixType>
template <class VectorType>
inline void
PreconditionSSOR<MatrixType>::vmult(VectorType &      dst,
                                    const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionSSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionSSOR<MatrixType>::Tvmult(VectorType &      dst,
                                     const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionSSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->precondition_SSOR(dst, src, this->relaxation, pos_right_of_diagonal);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionSSOR<MatrixType>::step(VectorType &dst, const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionSSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  this->A->SSOR_step(dst, src, this->relaxation);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionSSOR<MatrixType>::Tstep(VectorType &      dst,
                                    const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionSSOR<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionSSOR and VectorType must have the same size_type.");

  step(dst, src);
}



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

template <typename MatrixType>
inline void
PreconditionPSOR<MatrixType>::initialize(
  const MatrixType &                                                 rA,
  const std::vector<size_type> &                                     p,
  const std::vector<size_type> &                                     ip,
  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
{
  permutation         = &p;
  inverse_permutation = &ip;
  PreconditionRelaxation<MatrixType>::initialize(rA, parameters);
}


template <typename MatrixType>
inline void
PreconditionPSOR<MatrixType>::initialize(const MatrixType &    A,
                                         const AdditionalData &additional_data)
{
  initialize(A,
             additional_data.permutation,
             additional_data.inverse_permutation,
             additional_data.parameters);
}


template <typename MatrixType>
template <typename VectorType>
inline void
PreconditionPSOR<MatrixType>::vmult(VectorType &      dst,
                                    const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionPSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  dst = src;
  this->A->PSOR(dst, *permutation, *inverse_permutation, this->relaxation);
}



template <typename MatrixType>
template <class VectorType>
inline void
PreconditionPSOR<MatrixType>::Tvmult(VectorType &      dst,
                                     const VectorType &src) const
{
  static_assert(
    std::is_same<typename PreconditionPSOR<MatrixType>::size_type,
                 typename VectorType::size_type>::value,
    "PreconditionPSOR and VectorType must have the same size_type.");

  Assert(this->A != nullptr, ExcNotInitialized());
  dst = src;
  this->A->TPSOR(dst, *permutation, *inverse_permutation, this->relaxation);
}

template <typename MatrixType>
PreconditionPSOR<MatrixType>::AdditionalData::AdditionalData(
  const std::vector<size_type> &permutation,
  const std::vector<size_type> &inverse_permutation,
  const typename PreconditionRelaxation<MatrixType>::AdditionalData &parameters)
  : permutation(permutation)
  , inverse_permutation(inverse_permutation)
  , parameters(parameters)
{}


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


template <typename MatrixType, class VectorType>
PreconditionUseMatrix<MatrixType, VectorType>::PreconditionUseMatrix(
  const MatrixType & M,
  const function_ptr method)
  : matrix(M)
  , precondition(method)
{}



template <typename MatrixType, class VectorType>
void
PreconditionUseMatrix<MatrixType, VectorType>::vmult(
  VectorType &      dst,
  const VectorType &src) const
{
  (matrix.*precondition)(dst, src);
}

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

template <typename MatrixType>
inline PreconditionRelaxation<MatrixType>::AdditionalData::AdditionalData(
  const double relaxation)
  : relaxation(relaxation)
{}



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

namespace internal
{
  namespace PreconditionChebyshevImplementation
  {
    // for deal.II vectors, perform updates for Chebyshev preconditioner all
    // at once to reduce memory transfer. Here, we select between general
    // vectors and deal.II vectors where we expand the loop over the (local)
    // size of the vector

    // generic part for non-deal.II vectors
    template <typename VectorType, typename PreconditionerType>
    inline void
    vector_updates(const VectorType &        src,
                   const PreconditionerType &preconditioner,
                   const bool                start_zero,
                   const double              factor1,
                   const double              factor2,
                   VectorType &              update1,
                   VectorType &              update2,
                   VectorType &              update3,
                   VectorType &              dst)
    {
      if (start_zero)
        {
          update1.equ(factor2, src);
          preconditioner.vmult(dst, update1);
          update1.equ(-1., dst);
        }
      else
        {
          update2 -= src;
          preconditioner.vmult(update3, update2);
          update2 = update3;
          if (factor1 == 0.)
            update1.equ(factor2, update2);
          else
            update1.sadd(factor1, factor2, update2);
          dst -= update1;
        }
    }

    // worker routine for deal.II vectors. Because of vectorization, we need
    // to put the loop into an extra structure because the virtual function of
    // VectorUpdatesRange prevents the compiler from applying vectorization.
    template <typename Number>
    struct VectorUpdater
    {
      VectorUpdater(const Number *src,
                    const Number *matrix_diagonal_inverse,
                    const bool    start_zero,
                    const Number  factor1,
                    const Number  factor2,
                    Number *      update1,
                    Number *      update2,
                    Number *      dst)
        : src(src)
        , matrix_diagonal_inverse(matrix_diagonal_inverse)
        , do_startup(factor1 == Number())
        , start_zero(start_zero)
        , factor1(factor1)
        , factor2(factor2)
        , update1(update1)
        , update2(update2)
        , dst(dst)
      {}

      void
      apply_to_subrange(const std::size_t begin, const std::size_t end) const
      {
        // To circumvent a bug in gcc
        // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
        // copies of the variables factor1 and factor2 and do not check based on
        // factor1.
        const Number factor1 = this->factor1;
        const Number factor2 = this->factor2;
        if (do_startup)
          {
            if (start_zero)
              DEAL_II_OPENMP_SIMD_PRAGMA
            for (std::size_t i = begin; i < end; ++i)
              {
                dst[i]     = factor2 * src[i] * matrix_diagonal_inverse[i];
                update1[i] = -dst[i];
              }
            else DEAL_II_OPENMP_SIMD_PRAGMA for (std::size_t i = begin; i < end;
                                                 ++i)
            {
              update1[i] =
                ((update2[i] - src[i]) * factor2 * matrix_diagonal_inverse[i]);
              dst[i] -= update1[i];
            }
          }
        else
          DEAL_II_OPENMP_SIMD_PRAGMA
        for (std::size_t i = begin; i < end; ++i)
          {
            const Number update =
              factor1 * update1[i] +
              factor2 * ((update2[i] - src[i]) * matrix_diagonal_inverse[i]);
            update1[i] = update;
            dst[i] -= update;
          }
      }

      const Number *  src;
      const Number *  matrix_diagonal_inverse;
      const bool      do_startup;
      const bool      start_zero;
      const Number    factor1;
      const Number    factor2;
      mutable Number *update1;
      mutable Number *update2;
      mutable Number *dst;
    };

    template <typename Number>
    struct VectorUpdatesRange : public parallel::ParallelForInteger
    {
      VectorUpdatesRange(const VectorUpdater<Number> &updater,
                         const std::size_t            size)
        : updater(updater)
      {
        if (size < internal::VectorImplementation::minimum_parallel_grain_size)
          apply_to_subrange(0, size);
        else
          apply_parallel(
            0,
            size,
            internal::VectorImplementation::minimum_parallel_grain_size);
      }

      ~VectorUpdatesRange() override = default;

      virtual void
      apply_to_subrange(const std::size_t begin,
                        const std::size_t end) const override
      {
        updater.apply_to_subrange(begin, end);
      }

      const VectorUpdater<Number> &updater;
    };

    // selection for diagonal matrix around deal.II vector
    template <typename Number>
    inline void
    vector_updates(const ::Vector<Number> &                src,
                   const DiagonalMatrix<::Vector<Number>> &jacobi,
                   const bool                                      start_zero,
                   const double                                    factor1,
                   const double                                    factor2,
                   ::Vector<Number> &                      update1,
                   ::Vector<Number> &                      update2,
                   ::Vector<Number> &,
                   ::Vector<Number> &dst)
    {
      VectorUpdater<Number> upd(src.begin(),
                                jacobi.get_vector().begin(),
                                start_zero,
                                factor1,
                                factor2,
                                update1.begin(),
                                update2.begin(),
                                dst.begin());
      VectorUpdatesRange<Number>(upd, src.size());
    }

    // selection for diagonal matrix around parallel deal.II vector
    template <typename Number>
    inline void
    vector_updates(
      const LinearAlgebra::distributed::Vector<Number> &                src,
      const DiagonalMatrix<LinearAlgebra::distributed::Vector<Number>> &jacobi,
      const bool                                  start_zero,
      const double                                factor1,
      const double                                factor2,
      LinearAlgebra::distributed::Vector<Number> &update1,
      LinearAlgebra::distributed::Vector<Number> &update2,
      LinearAlgebra::distributed::Vector<Number> &,
      LinearAlgebra::distributed::Vector<Number> &dst)
    {
      VectorUpdater<Number> upd(src.begin(),
                                jacobi.get_vector().begin(),
                                start_zero,
                                factor1,
                                factor2,
                                update1.begin(),
                                update2.begin(),
                                dst.begin());
      VectorUpdatesRange<Number>(upd, src.local_size());
    }

    template <typename MatrixType,
              typename VectorType,
              typename PreconditionerType>
    inline void
    initialize_preconditioner(
      const MatrixType &                   matrix,
      std::shared_ptr<PreconditionerType> &preconditioner,
      VectorType &)
    {
      (void)matrix;
      (void)preconditioner;
      AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
    }

    template <typename MatrixType, typename VectorType>
    inline void
    initialize_preconditioner(
      const MatrixType &                           matrix,
      std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner,
      VectorType &                                 diagonal_inverse)
    {
      if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
        {
          if (preconditioner.get() == nullptr)
            preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();

          Assert(
            preconditioner->m() == 0,
            ExcMessage(
              "Preconditioner appears to be initialized but not sized correctly"));

          // Check if we can initialize from vector that then gets set to zero
          // as the matrix will own the memory
          preconditioner->reinit(diagonal_inverse);
          {
            VectorType empty_vector;
            diagonal_inverse.reinit(empty_vector);
          }

          // This part only works in serial
          if (preconditioner->m() != matrix.m())
            {
              preconditioner->get_vector().reinit(matrix.m());
              for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
                preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
            }
        }
    }

    template <typename VectorType>
    void
    set_initial_guess(VectorType &vector)
    {
      vector = 1. / std::sqrt(static_cast<double>(vector.size()));
      if (vector.locally_owned_elements().is_element(0))
        vector(0) = 0.;
    }

    template <typename Number>
    void
    set_initial_guess(::Vector<Number> &vector)
    {
      // Choose a high-frequency mode consisting of numbers between 0 and 1
      // that is cheap to compute (cheaper than random numbers) but avoids
      // obviously re-occurring numbers in multi-component systems by choosing
      // a period of 11
      for (unsigned int i = 0; i < vector.size(); ++i)
        vector(i) = i % 11;

      const Number mean_value = vector.mean_value();
      vector.add(-mean_value);
    }

    template <typename Number>
    void
    set_initial_guess(
      ::LinearAlgebra::distributed::Vector<Number> &vector)
    {
      // Choose a high-frequency mode consisting of numbers between 0 and 1
      // that is cheap to compute (cheaper than random numbers) but avoids
      // obviously re-occurring numbers in multi-component systems by choosing
      // a period of 11.
      // Make initial guess robust with respect to number of processors
      // by operating on the global index.
      types::global_dof_index first_local_range = 0;
      if (!vector.locally_owned_elements().is_empty())
        first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
      for (unsigned int i = 0; i < vector.local_size(); ++i)
        vector.local_element(i) = (i + first_local_range) % 11;

      const Number mean_value = vector.mean_value();
      vector.add(-mean_value);
    }

    struct EigenvalueTracker
    {
    public:
      void
      slot(const std::vector<double> &eigenvalues)
      {
        values = eigenvalues;
      }

      std::vector<double> values;
    };
  } // namespace PreconditionChebyshevImplementation
} // namespace internal



// avoid warning about deprecated variable nonzero_starting

template <typename MatrixType, class VectorType, typename PreconditionerType>
inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
  AdditionalData::AdditionalData(const unsigned int degree,
                                 const double       smoothing_range,
                                 const bool         nonzero_starting,
                                 const unsigned int eig_cg_n_iterations,
                                 const double       eig_cg_residual,
                                 const double       max_eigenvalue)
  : degree(degree)
  , smoothing_range(smoothing_range)
  , nonzero_starting(nonzero_starting)
  , eig_cg_n_iterations(eig_cg_n_iterations)
  , eig_cg_residual(eig_cg_residual)
  , max_eigenvalue(max_eigenvalue)
{}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
  PreconditionChebyshev()
  : theta(1.)
  , delta(1.)
  , eigenvalues_are_initialized(false)
{
  static_assert(
    std::is_same<size_type, typename VectorType::size_type>::value,
    "PreconditionChebyshev and VectorType must have the same size_type.");
}


// avoid warning about deprecated variable
// AdditionalData::matrix_diagonal_inverse

template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::initialize(
  const MatrixType &    matrix,
  const AdditionalData &additional_data)
{
  matrix_ptr = &matrix;
  data       = additional_data;
  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
    matrix, data.preconditioner, data.matrix_diagonal_inverse);
  eigenvalues_are_initialized = false;
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::clear()
{
  eigenvalues_are_initialized = false;
  theta = delta = 1.0;
  matrix_ptr    = nullptr;
  {
    VectorType empty_vector;
    data.matrix_diagonal_inverse.reinit(empty_vector);
    update1.reinit(empty_vector);
    update2.reinit(empty_vector);
    update3.reinit(empty_vector);
  }
  data.preconditioner.reset();
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
  estimate_eigenvalues(const VectorType &src) const
{
  Assert(eigenvalues_are_initialized == false, ExcInternalError());
  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());

  update1.reinit(src);
  update2.reinit(src, true);

  // calculate largest eigenvalue using a hand-tuned CG iteration on the
  // matrix weighted by its diagonal. we start with a vector that consists of
  // ones only, weighted by the length.
  double max_eigenvalue, min_eigenvalue;
  if (data.eig_cg_n_iterations > 0)
    {
      Assert(data.eig_cg_n_iterations > 2,
             ExcMessage(
               "Need to set at least two iterations to find eigenvalues."));

      // set a very strict tolerance to force at least two iterations
      ReductionControl control(
        data.eig_cg_n_iterations,
        std::sqrt(
          std::numeric_limits<typename VectorType::value_type>::epsilon()),
        1e-10,
        false,
        false);

      internal::PreconditionChebyshevImplementation::EigenvalueTracker
                           eigenvalue_tracker;
      SolverCG<VectorType> solver(control);
      solver.connect_eigenvalues_slot(std::bind(
        &internal::PreconditionChebyshevImplementation::EigenvalueTracker::slot,
        &eigenvalue_tracker,
        std::placeholders::_1));

      // set an initial guess which is close to the constant vector but where
      // one entry is different to trigger high frequencies
      internal::PreconditionChebyshevImplementation::set_initial_guess(update2);

      try
        {
          solver.solve(*matrix_ptr, update1, update2, *data.preconditioner);
        }
      catch (SolverControl::NoConvergence &)
        {}

      // read the eigenvalues from the attached eigenvalue tracker
      if (eigenvalue_tracker.values.empty())
        min_eigenvalue = max_eigenvalue = 1;
      else
        {
          min_eigenvalue = eigenvalue_tracker.values.front();

          // include a safety factor since the CG method will in general not
          // be converged
          max_eigenvalue = 1.2 * eigenvalue_tracker.values.back();
        }
    }
  else
    {
      max_eigenvalue = data.max_eigenvalue;
      min_eigenvalue = data.max_eigenvalue / data.smoothing_range;
    }

  const double alpha = (data.smoothing_range > 1. ?
                          max_eigenvalue / data.smoothing_range :
                          std::min(0.9 * max_eigenvalue, min_eigenvalue));

  // in case the user set the degree to invalid unsigned int, we have to
  // determine the number of necessary iterations from the Chebyshev error
  // estimate, given the target tolerance specified by smoothing_range. This
  // estimate is based on the error formula given in section 5.1 of
  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
  if (data.degree == numbers::invalid_unsigned_int)
    {
      const double actual_range = max_eigenvalue / alpha;
      const double sigma        = (1. - std::sqrt(1. / actual_range)) /
                           (1. + std::sqrt(1. / actual_range));
      const double eps = data.smoothing_range;
      const_cast<
        PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(
        this)
        ->data.degree = 1 + std::log(1. / eps + std::sqrt(1. / eps / eps - 1)) /
                              std::log(1. / sigma);
    }

  const_cast<
    PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
    ->delta = (max_eigenvalue - alpha) * 0.5;
  const_cast<
    PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
    ->theta = (max_eigenvalue + alpha) * 0.5;

  // We do not need the third auxiliary vector in case we have a
  // DiagonalMatrix as preconditioner and use deal.II's own vectors
  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
        false ||
      (std::is_same<VectorType,
                    ::Vector<typename VectorType::value_type>>::value ==
         false &&
       std::is_same<VectorType,
                    LinearAlgebra::distributed::Vector<
                      typename VectorType::value_type>>::value == false))
    update3.reinit(src, true);

  const_cast<
    PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
    ->eigenvalues_are_initialized = true;
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
  do_chebyshev_loop(VectorType &dst, const VectorType &src) const
{
  // if delta is zero, we do not need to iterate because the updates will be
  // zero
  if (std::abs(delta) < 1e-40)
    return;

  double rhok = delta / theta, sigma = theta / delta;
  for (unsigned int k = 0; k < data.degree; ++k)
    {
      matrix_ptr->vmult(update2, dst);
      const double rhokp   = 1. / (2. * sigma - rhok);
      const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
      rhok = rhokp;
      internal::PreconditionChebyshevImplementation::vector_updates(
        src,
        *data.preconditioner,
        false,
        factor1,
        factor2,
        update1,
        update2,
        update3,
        dst);
    }
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
  do_transpose_chebyshev_loop(VectorType &dst, const VectorType &src) const
{
  double rhok = delta / theta, sigma = theta / delta;
  for (unsigned int k = 0; k < data.degree; ++k)
    {
      matrix_ptr->Tvmult(update2, dst);
      const double rhokp   = 1. / (2. * sigma - rhok);
      const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
      rhok = rhokp;
      internal::PreconditionChebyshevImplementation::vector_updates(
        src,
        *data.preconditioner,
        false,
        factor1,
        factor2,
        update1,
        update2,
        update3,
        dst);
    }
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::vmult(
  VectorType &      dst,
  const VectorType &src) const
{
  Threads::Mutex::ScopedLock lock(mutex);
  if (eigenvalues_are_initialized == false)
    estimate_eigenvalues(src);

  internal::PreconditionChebyshevImplementation::vector_updates(
    src,
    *data.preconditioner,
    true,
    0.,
    1. / theta,
    update1,
    update2,
    update3,
    dst);

  do_chebyshev_loop(dst, src);
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tvmult(
  VectorType &      dst,
  const VectorType &src) const
{
  Threads::Mutex::ScopedLock lock(mutex);
  if (eigenvalues_are_initialized == false)
    estimate_eigenvalues(src);

  internal::PreconditionChebyshevImplementation::vector_updates(
    src,
    *data.preconditioner,
    true,
    0.,
    1. / theta,
    update1,
    update2,
    update3,
    dst);

  do_transpose_chebyshev_loop(dst, src);
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::step(
  VectorType &      dst,
  const VectorType &src) const
{
  Threads::Mutex::ScopedLock lock(mutex);
  if (eigenvalues_are_initialized == false)
    estimate_eigenvalues(src);

  matrix_ptr->vmult(update2, dst);
  internal::PreconditionChebyshevImplementation::vector_updates(
    src,
    *data.preconditioner,
    false,
    0.,
    1. / theta,
    update1,
    update2,
    update3,
    dst);

  do_chebyshev_loop(dst, src);
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline void
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tstep(
  VectorType &      dst,
  const VectorType &src) const
{
  Threads::Mutex::ScopedLock lock(mutex);
  if (eigenvalues_are_initialized == false)
    estimate_eigenvalues(src);

  matrix_ptr->Tvmult(update2, dst);
  internal::PreconditionChebyshevImplementation::vector_updates(
    src,
    *data.preconditioner,
    false,
    0.,
    1. / theta,
    update1,
    update2,
    update3,
    dst);

  do_transpose_chebyshev_loop(dst, src);
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline typename PreconditionChebyshev<MatrixType,
                                      VectorType,
                                      PreconditionerType>::size_type
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::m() const
{
  Assert(matrix_ptr != nullptr, ExcNotInitialized());
  return matrix_ptr->m();
}



template <typename MatrixType, typename VectorType, typename PreconditionerType>
inline typename PreconditionChebyshev<MatrixType,
                                      VectorType,
                                      PreconditionerType>::size_type
PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::n() const
{
  Assert(matrix_ptr != nullptr, ExcNotInitialized());
  return matrix_ptr->n();
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif