linear_operator.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2015 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_linear_operator_h
#define dealii_linear_operator_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
 
#include <deal.II/lac/vector_memory.h>
 
#include <array>
#include <functional>
#include <type_traits>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations:
#ifndef DOXYGEN
namespace internal
{
  namespace LinearOperatorImplementation
  {
    class EmptyPayload;
  }
} // namespace internal
 
template <typename Number>
class Vector;
 
class PreconditionIdentity;
 
template <typename Range  = Vector<double>,
          typename Domain = Range,
          typename Payload =
            internal::LinearOperatorImplementation::EmptyPayload>
class LinearOperator;
#endif
 
template <
  typename Range   = Vector<double>,
  typename Domain  = Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload,
  typename OperatorExemplar,
  typename Matrix>
LinearOperator<Range, Domain, Payload>
linear_operator(const OperatorExemplar &, const Matrix &);
 
template <
  typename Range   = Vector<double>,
  typename Domain  = Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload,
  typename Matrix>
LinearOperator<Range, Domain, Payload>
linear_operator(const Matrix &);
 
template <
  typename Range   = Vector<double>,
  typename Domain  = Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator<Range, Domain, Payload>
null_operator(const LinearOperator<Range, Domain, Payload> &);
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Range, Domain, Payload>
identity_operator(const LinearOperator<Range, Domain, Payload> &);
 
 
template <typename Range, typename Domain, typename Payload>
class LinearOperator : public Payload
{
public:
  LinearOperator(const Payload &payload = Payload())
    : Payload(payload)
    , is_null_operator(false)
  {
    vmult = [](Range &, const Domain &) {
      Assert(false,
             ExcMessage("Uninitialized LinearOperator<Range, "
                        "Domain>::vmult called"));
    };
 
    vmult_add = [](Range &, const Domain &) {
      Assert(false,
             ExcMessage("Uninitialized LinearOperator<Range, "
                        "Domain>::vmult_add called"));
    };
 
    Tvmult = [](Domain &, const Range &) {
      Assert(false,
             ExcMessage("Uninitialized LinearOperator<Range, "
                        "Domain>::Tvmult called"));
    };
 
    Tvmult_add = [](Domain &, const Range &) {
      Assert(false,
             ExcMessage("Uninitialized LinearOperator<Range, "
                        "Domain>::Tvmult_add called"));
    };
 
    reinit_range_vector = [](Range &, bool) {
      Assert(false,
             ExcMessage("Uninitialized LinearOperator<Range, "
                        "Domain>::reinit_range_vector method called"));
    };
 
    reinit_domain_vector = [](Domain &, bool) {
      Assert(false,
             ExcMessage("Uninitialized LinearOperator<Range, "
                        "Domain>::reinit_domain_vector method called"));
    };
  }
 
  LinearOperator(const LinearOperator<Range, Domain, Payload> &) = default;
 
  template <typename Op,
            typename = std::enable_if_t<
              !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
  LinearOperator(const Op &op)
  {
    *this = linear_operator<Range, Domain, Payload, Op>(op);
  }
 
  LinearOperator<Range, Domain, Payload> &
  operator=(const LinearOperator<Range, Domain, Payload> &) = default;
 
  template <typename Op,
            typename = std::enable_if_t<
              !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
  LinearOperator<Range, Domain, Payload> &
  operator=(const Op &op)
  {
    *this = linear_operator<Range, Domain, Payload, Op>(op);
    return *this;
  }
 
  std::function<void(Range &v, const Domain &u)> vmult;
 
  std::function<void(Range &v, const Domain &u)> vmult_add;
 
  std::function<void(Domain &v, const Range &u)> Tvmult;
 
  std::function<void(Domain &v, const Range &u)> Tvmult_add;
 
  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_range_vector;
 
  std::function<void(Domain &v, bool omit_zeroing_entries)>
    reinit_domain_vector;
 
  LinearOperator<Range, Domain, Payload> &
  operator+=(const LinearOperator<Range, Domain, Payload> &second_op)
  {
    *this = *this + second_op;
    return *this;
  }
 
  LinearOperator<Range, Domain, Payload> &
  operator-=(const LinearOperator<Range, Domain, Payload> &second_op)
  {
    *this = *this - second_op;
    return *this;
  }
 
  LinearOperator<Range, Domain, Payload> &
  operator*=(const LinearOperator<Domain, Domain, Payload> &second_op)
  {
    *this = *this * second_op;
    return *this;
  }
 
  LinearOperator<Range, Domain, Payload>
  operator*=(typename Domain::value_type number)
  {
    *this = *this * number;
    return *this;
  }
 
  bool is_null_operator;
 
};
 
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Range, Domain, Payload>
operator+(const LinearOperator<Range, Domain, Payload> &first_op,
          const LinearOperator<Range, Domain, Payload> &second_op)
{
  if (first_op.is_null_operator)
    {
      return second_op;
    }
  else if (second_op.is_null_operator)
    {
      return first_op;
    }
  else
    {
      LinearOperator<Range, Domain, Payload> return_op{
        static_cast<const Payload &>(first_op) +
        static_cast<const Payload &>(second_op)};
 
      return_op.reinit_range_vector  = first_op.reinit_range_vector;
      return_op.reinit_domain_vector = first_op.reinit_domain_vector;
 
      // ensure to have valid computation objects by catching first_op and
      // second_op by value
 
      return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
        first_op.vmult(v, u);
        second_op.vmult_add(v, u);
      };
 
      return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
        first_op.vmult_add(v, u);
        second_op.vmult_add(v, u);
      };
 
      return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
        second_op.Tvmult(v, u);
        first_op.Tvmult_add(v, u);
      };
 
      return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
        second_op.Tvmult_add(v, u);
        first_op.Tvmult_add(v, u);
      };
 
      return return_op;
    }
}
 
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Range, Domain, Payload>
operator-(const LinearOperator<Range, Domain, Payload> &first_op,
          const LinearOperator<Range, Domain, Payload> &second_op)
{
  if (first_op.is_null_operator)
    {
      return -1. * second_op;
    }
  else if (second_op.is_null_operator)
    {
      return first_op;
    }
  else
    {
      // implement with addition and scalar multiplication
      return first_op + (-1. * second_op);
    }
}
 
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Range, Domain, Payload>
operator*(typename Range::value_type                    number,
          const LinearOperator<Range, Domain, Payload> &op)
{
  static_assert(
    std::is_convertible_v<typename Range::value_type,
                          typename Domain::value_type>,
    "Range and Domain must have implicitly convertible 'value_type's");
 
  if (op.is_null_operator)
    {
      return op;
    }
  else if (number == 0.)
    {
      return null_operator(op);
    }
  else
    {
      LinearOperator<Range, Domain, Payload> return_op = op;
 
      // ensure to have valid computation objects by catching number and op by
      // value
 
      return_op.vmult = [number, op](Range &v, const Domain &u) {
        op.vmult(v, u);
        v *= number;
      };
 
      return_op.vmult_add = [number, op](Range &v, const Domain &u) {
        v /= number;
        op.vmult_add(v, u);
        v *= number;
      };
 
      return_op.Tvmult = [number, op](Domain &v, const Range &u) {
        op.Tvmult(v, u);
        v *= number;
      };
 
      return_op.Tvmult_add = [number, op](Domain &v, const Range &u) {
        v /= number;
        op.Tvmult_add(v, u);
        v *= number;
      };
 
      return return_op;
    }
}
 
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Range, Domain, Payload>
operator*(const LinearOperator<Range, Domain, Payload> &op,
          typename Domain::value_type                   number)
{
  static_assert(
    std::is_convertible_v<typename Range::value_type,
                          typename Domain::value_type>,
    "Range and Domain must have implicitly convertible 'value_type's");
 
  return number * op;
}
 
template <typename Range,
          typename Intermediate,
          typename Domain,
          typename Payload>
LinearOperator<Range, Domain, Payload>
operator*(const LinearOperator<Range, Intermediate, Payload>  &first_op,
          const LinearOperator<Intermediate, Domain, Payload> &second_op)
{
  if (first_op.is_null_operator || second_op.is_null_operator)
    {
      LinearOperator<Range, Domain, Payload> return_op;
      return_op.reinit_domain_vector = second_op.reinit_domain_vector;
      return_op.reinit_range_vector  = first_op.reinit_range_vector;
      return null_operator(return_op);
    }
  else
    {
      LinearOperator<Range, Domain, Payload> return_op{
        static_cast<const Payload &>(first_op) *
        static_cast<const Payload &>(second_op)};
 
      return_op.reinit_domain_vector = second_op.reinit_domain_vector;
      return_op.reinit_range_vector  = first_op.reinit_range_vector;
 
      // ensure to have valid computation objects by catching first_op and
      // second_op by value
 
      return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
        GrowingVectorMemory<Intermediate> vector_memory;
 
        typename VectorMemory<Intermediate>::Pointer i(vector_memory);
        second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
        second_op.vmult(*i, u);
        first_op.vmult(v, *i);
      };
 
      return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
        GrowingVectorMemory<Intermediate> vector_memory;
 
        typename VectorMemory<Intermediate>::Pointer i(vector_memory);
        second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
        second_op.vmult(*i, u);
        first_op.vmult_add(v, *i);
      };
 
      return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
        GrowingVectorMemory<Intermediate> vector_memory;
 
        typename VectorMemory<Intermediate>::Pointer i(vector_memory);
        first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
        first_op.Tvmult(*i, u);
        second_op.Tvmult(v, *i);
      };
 
      return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
        GrowingVectorMemory<Intermediate> vector_memory;
 
        typename VectorMemory<Intermediate>::Pointer i(vector_memory);
        first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
        first_op.Tvmult(*i, u);
        second_op.Tvmult_add(v, *i);
      };
 
      return return_op;
    }
}
 
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Domain, Range, Payload>
transpose_operator(const LinearOperator<Range, Domain, Payload> &op)
{
  LinearOperator<Domain, Range, Payload> return_op{op.transpose_payload()};
 
  return_op.reinit_range_vector  = op.reinit_domain_vector;
  return_op.reinit_domain_vector = op.reinit_range_vector;
 
  return_op.vmult      = op.Tvmult;
  return_op.vmult_add  = op.Tvmult_add;
  return_op.Tvmult     = op.vmult;
  return_op.Tvmult_add = op.vmult_add;
 
  return return_op;
}
 
 
template <typename Payload,
          typename Solver,
          typename Preconditioner,
          typename Range  = typename Solver::vector_type,
          typename Domain = Range>
LinearOperator<Domain, Range, Payload>
inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
                 Solver                                       &solver,
                 const Preconditioner                         &preconditioner)
{
  LinearOperator<Domain, Range, Payload> return_op{
    op.inverse_payload(solver, preconditioner)};
 
  return_op.reinit_range_vector  = op.reinit_domain_vector;
  return_op.reinit_domain_vector = op.reinit_range_vector;
 
  return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
    op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
    solver.solve(op, v, u, preconditioner);
  };
 
  return_op.vmult_add = [op, &solver, &preconditioner](Range        &v,
                                                       const Domain &u) {
    GrowingVectorMemory<Range> vector_memory;
 
    typename VectorMemory<Range>::Pointer v2(vector_memory);
    op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
    solver.solve(op, *v2, u, preconditioner);
    v += *v2;
  };
 
  return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
    op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
    solver.solve(transpose_operator(op), v, u, preconditioner);
  };
 
  return_op.Tvmult_add = [op, &solver, &preconditioner](Range        &v,
                                                        const Domain &u) {
    GrowingVectorMemory<Range> vector_memory;
 
    typename VectorMemory<Range>::Pointer v2(vector_memory);
    op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
    solver.solve(transpose_operator(op), *v2, u, preconditioner);
    v += *v2;
  };
 
  return return_op;
}
 
 
template <typename Payload,
          typename Solver,
          typename Range  = typename Solver::vector_type,
          typename Domain = Range>
LinearOperator<Domain, Range, Payload>
inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
                 Solver                                       &solver,
                 const LinearOperator<Range, Domain, Payload> &preconditioner)
{
  LinearOperator<Domain, Range, Payload> return_op{
    op.inverse_payload(solver, preconditioner)};
 
  return_op.reinit_range_vector  = op.reinit_domain_vector;
  return_op.reinit_domain_vector = op.reinit_range_vector;
 
  return_op.vmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
    op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
    solver.solve(op, v, u, preconditioner);
  };
 
  return_op.vmult_add = [op, &solver, preconditioner](Range        &v,
                                                      const Domain &u) {
    GrowingVectorMemory<Range> vector_memory;
 
    typename VectorMemory<Range>::Pointer v2(vector_memory);
    op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
    solver.solve(op, *v2, u, preconditioner);
    v += *v2;
  };
 
  return_op.Tvmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
    op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
    solver.solve(transpose_operator(op), v, u, preconditioner);
  };
 
  return_op.Tvmult_add = [op, &solver, preconditioner](Range        &v,
                                                       const Domain &u) {
    GrowingVectorMemory<Range> vector_memory;
 
    typename VectorMemory<Range>::Pointer v2(vector_memory);
    op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
    solver.solve(transpose_operator(op), *v2, u, preconditioner);
    v += *v2;
  };
 
  return return_op;
}
 
 
template <typename Payload,
          typename Solver,
          typename Range  = typename Solver::vector_type,
          typename Domain = Range>
LinearOperator<Domain, Range, Payload>
inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
                 Solver                                       &solver)
{
  return inverse_operator(op, solver, identity_operator(op));
}
 
 
template <typename Payload,
          typename Solver,
          typename Range  = typename Solver::vector_type,
          typename Domain = Range>
LinearOperator<Domain, Range, Payload>
inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
                 Solver                                       &solver,
                 const PreconditionIdentity &)
{
  return inverse_operator(op, solver);
}
 
template <
  typename Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator<Range, Range, Payload>
identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
{
  LinearOperator<Range, Range, Payload> return_op{Payload()};
 
  return_op.reinit_range_vector  = reinit_vector;
  return_op.reinit_domain_vector = reinit_vector;
 
  return_op.vmult = [](Range &v, const Range &u) { v = u; };
 
  return_op.vmult_add = [](Range &v, const Range &u) { v += u; };
 
  return_op.Tvmult = [](Range &v, const Range &u) { v = u; };
 
  return_op.Tvmult_add = [](Range &v, const Range &u) { v += u; };
 
  return return_op;
}
 
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Range, Domain, Payload>
identity_operator(const LinearOperator<Range, Domain, Payload> &op)
{
  auto return_op = identity_operator<Range, Payload>(op.reinit_range_vector);
  static_cast<Payload &>(return_op) = op.identity_payload();
 
  return return_op;
}
 
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Range, Domain, Payload>
null_operator(const LinearOperator<Range, Domain, Payload> &op)
{
  LinearOperator<Range, Domain, Payload> return_op{op.null_payload()};
 
  return_op.is_null_operator = true;
 
  return_op.reinit_range_vector  = op.reinit_range_vector;
  return_op.reinit_domain_vector = op.reinit_domain_vector;
 
  return_op.vmult = [](Range &v, const Domain &) { v = 0.; };
 
  return_op.vmult_add = [](Range &, const Domain &) {};
 
  return_op.Tvmult = [](Domain &v, const Range &) { v = 0.; };
 
  return_op.Tvmult_add = [](Domain &, const Range &) {};
 
  return return_op;
}
 
template <typename Range,
          typename Domain,
          typename Payload,
          typename VectorType>
LinearOperator<Range, Domain, Payload>
diagonal_operator(const LinearOperator<Range, Domain, Payload> &exemplar,
                  const VectorType                             &diagonal)
{
  LinearOperator<Range, Domain, Payload> return_op;
  return_op.reinit_range_vector  = exemplar.reinit_range_vector;
  return_op.reinit_domain_vector = exemplar.reinit_domain_vector;
 
  return_op.vmult = [&diagonal](Range &dst, const Domain &src) {
    Assert(src.size() == diagonal.size(),
           ExcDimensionMismatch(src.size(), diagonal.size()));
    Assert(dst.size() == diagonal.size(),
           ExcDimensionMismatch(dst.size(), diagonal.size()));
    dst = src;
    dst.scale(diagonal);
  };
 
  return_op.vmult_add = [&diagonal](Range &dst, const Domain &src) {
    Assert(src.size() == diagonal.size(),
           ExcDimensionMismatch(src.size(), diagonal.size()));
    Assert(dst.size() == diagonal.size(),
           ExcDimensionMismatch(dst.size(), diagonal.size()));
    Range tmp;
    tmp = src;
    tmp.scale(diagonal);
    dst += tmp;
  };
 
  return_op.Tvmult = [&diagonal](Domain &dst, const Range &src) {
    Assert(src.size() == diagonal.size(),
           ExcDimensionMismatch(src.size(), diagonal.size()));
    Assert(dst.size() == diagonal.size(),
           ExcDimensionMismatch(dst.size(), diagonal.size()));
    dst = src;
    dst.scale(diagonal);
  };
 
  return_op.Tvmult_add = [&diagonal](Domain &dst, const Range &src) {
    Assert(src.size() == diagonal.size(),
           ExcDimensionMismatch(src.size(), diagonal.size()));
    Assert(dst.size() == diagonal.size(),
           ExcDimensionMismatch(dst.size(), diagonal.size()));
    Domain tmp;
    tmp = src;
    tmp.scale(diagonal);
    dst += tmp;
  };
  return return_op;
}
 
 
 
template <
  typename Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload>
LinearOperator<Range, Range, Payload>
mean_value_filter(const std::function<void(Range &, bool)> &reinit_vector)
{
  LinearOperator<Range, Range, Payload> return_op{Payload()};
 
  return_op.reinit_range_vector  = reinit_vector;
  return_op.reinit_domain_vector = reinit_vector;
 
  return_op.vmult = [](Range &v, const Range &u) {
    const auto mean = u.mean_value();
 
    v = u;
    v.add(-mean);
  };
 
  return_op.vmult_add = [](Range &v, const Range &u) {
    const auto mean = u.mean_value();
 
    v += u;
    v.add(-mean);
  };
 
  return_op.Tvmult     = return_op.vmult_add;
  return_op.Tvmult_add = return_op.vmult_add;
 
  return return_op;
}
 
 
template <typename Range, typename Domain, typename Payload>
LinearOperator<Range, Domain, Payload>
mean_value_filter(const LinearOperator<Range, Domain, Payload> &op)
{
  auto return_op = mean_value_filter<Range, Payload>(op.reinit_range_vector);
  static_cast<Payload &>(return_op) = op.identity_payload();
 
  return return_op;
}
 
 
namespace internal
{
  namespace LinearOperatorImplementation
  {
    template <typename Vector>
    class ReinitHelper
    {
    public:
      template <typename Matrix>
      static void
      reinit_range_vector(const Matrix &matrix,
                          Vector       &v,
                          bool          omit_zeroing_entries)
      {
        v.reinit(matrix.m(), omit_zeroing_entries);
      }
 
      template <typename Matrix>
      static void
      reinit_domain_vector(const Matrix &matrix,
                           Vector       &v,
                           bool          omit_zeroing_entries)
      {
        v.reinit(matrix.n(), omit_zeroing_entries);
      }
    };
 
 
    class EmptyPayload
    {
    public:
      template <typename... Args>
      EmptyPayload(const Args &...)
      {}
 
 
      EmptyPayload
      identity_payload() const
      {
        return *this;
      }
 
 
      EmptyPayload
      null_payload() const
      {
        return *this;
      }
 
 
      EmptyPayload
      transpose_payload() const
      {
        return *this;
      }
 
 
      template <typename Solver, typename Preconditioner>
      EmptyPayload
      inverse_payload(Solver &, const Preconditioner &) const
      {
        return *this;
      }
    };
 
    inline EmptyPayload
    operator+(const EmptyPayload &, const EmptyPayload &)
    {
      return {};
    }
 
    inline EmptyPayload
    operator*(const EmptyPayload &, const EmptyPayload &)
    {
      return {};
    }
 
 
 
    // A trait class that determines whether type T provides public
    // (templated or non-templated) vmult_add member functions
    template <typename Range, typename Domain, typename T>
    class has_vmult_add_and_Tvmult_add
    {
      template <typename C>
      static std::false_type
      test(...);
 
      template <typename C>
      static auto
      test(Range *r, Domain *d)
        -> decltype(std::declval<C>().vmult_add(*r, *d),
                    std::declval<C>().Tvmult_add(*d, *r),
                    std::true_type());
 
    public:
      // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
      // otherwise it is std::false_type
 
      using type = decltype(test<T>(nullptr, nullptr));
    };
 
 
    // A helper function to apply a given vmult, or Tvmult to a vector with
    // intermediate storage
    template <typename Function, typename Range, typename Domain>
    void
    apply_with_intermediate_storage(Function      function,
                                    Range        &v,
                                    const Domain &u,
                                    bool          add)
    {
      GrowingVectorMemory<Range> vector_memory;
 
      typename VectorMemory<Range>::Pointer i(vector_memory);
      i->reinit(v, /*bool omit_zeroing_entries =*/true);
 
      function(*i, u);
 
      if (add)
        v += *i;
      else
        v = *i;
    }
 
 
    // A helper class to add a reduced matrix interface to a LinearOperator
    // (typically provided by Preconditioner classes)
    template <typename Range, typename Domain, typename Payload>
    class MatrixInterfaceWithoutVmultAdd
    {
    public:
      template <typename Matrix>
      void
      operator()(LinearOperator<Range, Domain, Payload> &op,
                 const Matrix                           &matrix)
      {
        op.vmult = [&matrix](Range &v, const Domain &u) {
          if (PointerComparison::equal(&v, &u))
            {
              // If v and u are the same memory location use intermediate
              // storage
              apply_with_intermediate_storage(
                [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
                v,
                u,
                /*bool add =*/false);
            }
          else
            {
              matrix.vmult(v, u);
            }
        };
 
        op.vmult_add = [&matrix](Range &v, const Domain &u) {
          // use intermediate storage to implement vmult_add with vmult
          apply_with_intermediate_storage(
            [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
            v,
            u,
            /*bool add =*/true);
        };
 
        op.Tvmult = [&matrix](Domain &v, const Range &u) {
          if (PointerComparison::equal(&v, &u))
            {
              // If v and u are the same memory location use intermediate
              // storage
              apply_with_intermediate_storage(
                [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
                v,
                u,
                /*bool add =*/false);
            }
          else
            {
              matrix.Tvmult(v, u);
            }
        };
 
        op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
          // use intermediate storage to implement Tvmult_add with Tvmult
          apply_with_intermediate_storage(
            [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
            v,
            u,
            /*bool add =*/true);
        };
      }
    };
 
 
    // A helper class to add the full matrix interface to a LinearOperator
    template <typename Range, typename Domain, typename Payload>
    class MatrixInterfaceWithVmultAdd
    {
    public:
      template <typename Matrix>
      void
      operator()(LinearOperator<Range, Domain, Payload> &op,
                 const Matrix                           &matrix)
      {
        // As above ...
 
        MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>().operator()(
          op, matrix);
 
        // ... but add native vmult_add and Tvmult_add variants:
 
        op.vmult_add = [&matrix](Range &v, const Domain &u) {
          if (PointerComparison::equal(&v, &u))
            {
              apply_with_intermediate_storage(
                [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
                v,
                u,
                /*bool add =*/true);
            }
          else
            {
              matrix.vmult_add(v, u);
            }
        };
 
        op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
          if (PointerComparison::equal(&v, &u))
            {
              apply_with_intermediate_storage(
                [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
                v,
                u,
                /*bool add =*/true);
            }
          else
            {
              matrix.Tvmult_add(v, u);
            }
        };
      }
    };
  } // namespace LinearOperatorImplementation
} // namespace internal
 
 
template <typename Range, typename Domain, typename Payload, typename Matrix>
LinearOperator<Range, Domain, Payload>
linear_operator(const Matrix &matrix)
{
  // implement with the more generic variant below...
  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
                                                                 matrix);
}
 
 
template <typename Range,
          typename Domain,
          typename Payload,
          typename OperatorExemplar,
          typename Matrix>
LinearOperator<Range, Domain, Payload>
linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
{
  using namespace internal::LinearOperatorImplementation;
  // Initialize the payload based on the input exemplar matrix
  LinearOperator<Range, Domain, Payload> return_op{
    Payload(operator_exemplar, matrix)};
 
  // Always store a reference to matrix and operator_exemplar in the lambda
  // functions. This ensures that a modification of the matrix after the
  // creation of a LinearOperator wrapper is respected - further a matrix
  // or an operator_exemplar cannot usually be copied...
 
  return_op.reinit_range_vector =
    [&operator_exemplar](Range &v, bool omit_zeroing_entries) {
      internal::LinearOperatorImplementation::ReinitHelper<
        Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
    };
 
  return_op.reinit_domain_vector = [&operator_exemplar](
                                     Domain &v, bool omit_zeroing_entries) {
    internal::LinearOperatorImplementation::ReinitHelper<
      Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
  };
 
  std::conditional_t<
    has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
    MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
    MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>()
    .
    operator()(return_op, matrix);
 
  return return_op;
}
 
 
 
template <typename Range, typename Domain, typename Payload, typename Matrix>
LinearOperator<Range, Domain, Payload>
linear_operator(const LinearOperator<Range, Domain, Payload> &operator_exemplar,
                const Matrix                                 &matrix)
{
  using namespace internal::LinearOperatorImplementation;
  // Initialize the payload based on the LinearOperator exemplar
  auto return_op = operator_exemplar;
 
  std::conditional_t<
    has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
    MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
    MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>()
    .
    operator()(return_op, matrix);
 
  return return_op;
}
 
 
#ifndef DOXYGEN
 
//
// Ensure that we never capture a reference to a temporary by accident.
// to avoid "stack use after free".
//
 
template <
  typename Range   = Vector<double>,
  typename Domain  = Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload,
  typename OperatorExemplar,
  typename Matrix,
  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
LinearOperator<Range, Domain, Payload>
linear_operator(const OperatorExemplar &, Matrix &&) = delete;
 
template <
  typename Range   = Vector<double>,
  typename Domain  = Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload,
  typename OperatorExemplar,
  typename Matrix,
  typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
  typename = std::enable_if_t<
    !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
LinearOperator<Range, Domain, Payload>
linear_operator(OperatorExemplar &&, const Matrix &) = delete;
 
template <
  typename Range   = Vector<double>,
  typename Domain  = Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload,
  typename OperatorExemplar,
  typename Matrix,
  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>,
  typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
  typename = std::enable_if_t<
    !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
LinearOperator<Range, Domain, Payload>
linear_operator(OperatorExemplar &&, Matrix &&) = delete;
 
template <
  typename Range   = Vector<double>,
  typename Domain  = Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload,
  typename Matrix,
  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
LinearOperator<Range, Domain, Payload>
linear_operator(const LinearOperator<Range, Domain, Payload> &,
                Matrix &&) = delete;
 
template <
  typename Range   = Vector<double>,
  typename Domain  = Range,
  typename Payload = internal::LinearOperatorImplementation::EmptyPayload,
  typename Matrix,
  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
LinearOperator<Range, Domain, Payload>
linear_operator(Matrix &&) = delete;
 
template <
  typename Payload,
  typename Solver,
  typename Preconditioner,
  typename Range  = typename Solver::vector_type,
  typename Domain = Range,
  typename = std::enable_if_t<!std::is_lvalue_reference_v<Preconditioner>>,
  typename =
    std::enable_if_t<!std::is_same_v<Preconditioner, PreconditionIdentity>>,
  typename = std::enable_if_t<
    !std::is_same_v<Preconditioner, LinearOperator<Range, Domain, Payload>>>>
LinearOperator<Domain, Range, Payload>
inverse_operator(const LinearOperator<Range, Domain, Payload> &,
                 Solver &,
                 Preconditioner &&) = delete;
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif