16 #ifndef dealii_linear_operator_h
17 #define dealii_linear_operator_h
27 #include <type_traits>
35 namespace LinearOperatorImplementation
41 template <
typename Number>
46 template <
typename Range = Vector<
double>,
47 typename Domain = Range,
49 internal::LinearOperatorImplementation::EmptyPayload>
55 typename Domain = Range,
57 typename OperatorExemplar,
64 typename Domain = Range,
72 typename Domain = Range,
77 template <
typename Range,
typename Domain,
typename Payload>
178 template <
typename Range,
typename Domain,
typename Payload>
194 vmult = [](Range &,
const Domain &) {
196 ExcMessage(
"Uninitialized LinearOperator<Range, "
197 "Domain>::vmult called"));
200 vmult_add = [](Range &,
const Domain &) {
202 ExcMessage(
"Uninitialized LinearOperator<Range, "
203 "Domain>::vmult_add called"));
206 Tvmult = [](Domain &,
const Range &) {
208 ExcMessage(
"Uninitialized LinearOperator<Range, "
209 "Domain>::Tvmult called"));
214 ExcMessage(
"Uninitialized LinearOperator<Range, "
215 "Domain>::Tvmult_add called"));
220 ExcMessage(
"Uninitialized LinearOperator<Range, "
221 "Domain>::reinit_range_vector method called"));
226 ExcMessage(
"Uninitialized LinearOperator<Range, "
227 "Domain>::reinit_domain_vector method called"));
241 template <
typename Op,
242 typename = std::enable_if_t<
243 !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
246 *
this = linear_operator<Range, Domain, Payload, Op>(op);
259 template <
typename Op,
260 typename = std::enable_if_t<
261 !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
265 *
this = linear_operator<Range, Domain, Payload, Op>(op);
273 std::function<void(Range &v,
const Domain &u)>
vmult;
279 std::function<void(Range &v,
const Domain &u)>
vmult_add;
285 std::function<void(Domain &v,
const Range &u)>
Tvmult;
309 std::function<void(Domain &v,
bool omit_zeroing_entries)>
324 *
this = *
this + second_op;
335 *
this = *
this - second_op;
346 *
this = *
this * second_op;
357 *
this = *
this * number;
386 template <
typename Range,
typename Domain,
typename Payload>
402 static_cast<const Payload &
>(first_op) +
403 static_cast<const Payload &
>(second_op)};
411 return_op.vmult = [first_op, second_op](Range &v,
const Domain &u) {
412 first_op.
vmult(v, u);
416 return_op.vmult_add = [first_op, second_op](Range &v,
const Domain &u) {
421 return_op.Tvmult = [first_op, second_op](Domain &v,
const Range &u) {
426 return_op.Tvmult_add = [first_op, second_op](Domain &v,
const Range &u) {
445 template <
typename Range,
typename Domain,
typename Payload>
452 return -1. * second_op;
461 return first_op + (-1. * second_op);
483 template <
typename Range,
typename Domain,
typename Payload>
489 std::is_convertible<
typename Range::value_type,
490 typename Domain::value_type>::value,
491 "Range and Domain must have implicitly convertible 'value_type's");
497 else if (number == 0.)
508 return_op.
vmult = [number, op](Range &v,
const Domain &u) {
513 return_op.
vmult_add = [number, op](Range &v,
const Domain &u) {
519 return_op.
Tvmult = [number, op](Domain &v,
const Range &u) {
524 return_op.
Tvmult_add = [number, op](Domain &v,
const Range &u) {
550 template <
typename Range,
typename Domain,
typename Payload>
553 typename Domain::value_type number)
556 std::is_convertible<
typename Range::value_type,
557 typename Domain::value_type>::value,
558 "Range and Domain must have implicitly convertible 'value_type's");
580 template <
typename Range,
581 typename Intermediate,
598 static_cast<const Payload &
>(first_op) *
599 static_cast<const Payload &
>(second_op)};
607 return_op.vmult = [first_op, second_op](Range &v,
const Domain &u) {
612 second_op.
vmult(*i, u);
613 first_op.
vmult(v, *i);
616 return_op.vmult_add = [first_op, second_op](Range &v,
const Domain &u) {
621 second_op.
vmult(*i, u);
625 return_op.Tvmult = [first_op, second_op](Domain &v,
const Range &u) {
634 return_op.Tvmult_add = [first_op, second_op](Domain &v,
const Range &u) {
656 template <
typename Range,
typename Domain,
typename Payload>
665 return_op.vmult = op.
Tvmult;
667 return_op.Tvmult = op.
vmult;
695 typename Preconditioner,
696 typename Range =
typename Solver::vector_type,
697 typename Domain = Range>
701 const Preconditioner &preconditioner)
704 op.inverse_payload(solver, preconditioner)};
709 return_op.vmult = [op, &solver, &preconditioner](Range &v,
const Domain &u) {
711 solver.solve(op, v, u, preconditioner);
714 return_op.vmult_add = [op, &solver, &preconditioner](Range &v,
720 solver.solve(op, *v2, u, preconditioner);
724 return_op.Tvmult = [op, &solver, &preconditioner](Range &v,
const Domain &u) {
729 return_op.Tvmult_add = [op, &solver, &preconditioner](Range &v,
753 typename Range =
typename Solver::vector_type,
754 typename Domain = Range>
761 op.inverse_payload(solver, preconditioner)};
766 return_op.vmult = [op, &solver, preconditioner](Range &v,
const Domain &u) {
768 solver.solve(op, v, u, preconditioner);
771 return_op.vmult_add = [op, &solver, preconditioner](Range &v,
777 solver.solve(op, *v2, u, preconditioner);
781 return_op.Tvmult = [op, &solver, preconditioner](Range &v,
const Domain &u) {
786 return_op.Tvmult_add = [op, &solver, preconditioner](Range &v,
811 typename Range =
typename Solver::vector_type,
812 typename Domain = Range>
831 typename Range =
typename Solver::vector_type,
832 typename Domain = Range>
868 return_op.reinit_range_vector = reinit_vector;
869 return_op.reinit_domain_vector = reinit_vector;
871 return_op.vmult = [](Range &v,
const Range &u) { v = u; };
873 return_op.vmult_add = [](Range &v,
const Range &u) { v += u; };
875 return_op.Tvmult = [](Range &v,
const Range &u) { v = u; };
877 return_op.Tvmult_add = [](Range &v,
const Range &u) { v += u; };
895 template <
typename Range,
typename Domain,
typename Payload>
900 static_cast<Payload &
>(return_op) = op.identity_payload();
915 template <
typename Range,
typename Domain,
typename Payload>
926 return_op.vmult = [](Range &v,
const Domain &) { v = 0.; };
928 return_op.vmult_add = [](Range &,
const Domain &) {};
930 return_op.Tvmult = [](Domain &v,
const Range &) { v = 0.; };
932 return_op.Tvmult_add = [](Domain &,
const Range &) {};
958 return_op.reinit_range_vector = reinit_vector;
959 return_op.reinit_domain_vector = reinit_vector;
961 return_op.vmult = [](Range &v,
const Range &u) {
962 const auto mean = u.mean_value();
968 return_op.vmult_add = [](Range &v,
const Range &u) {
969 const auto mean = u.mean_value();
975 return_op.Tvmult = return_op.vmult_add;
976 return_op.Tvmult_add = return_op.vmult_add;
994 template <
typename Range,
typename Domain,
typename Payload>
999 static_cast<Payload &
>(return_op) = op.identity_payload();
1007 namespace LinearOperatorImplementation
1020 template <
typename Vector>
1035 template <
typename Matrix>
1039 bool omit_zeroing_entries)
1055 template <
typename Matrix>
1059 bool omit_zeroing_entries)
1087 template <
typename... Args>
1125 template <
typename Solver,
typename Preconditioner>
1157 template <
typename Range,
typename Domain,
typename T>
1160 template <
typename C>
1161 static std::false_type
1164 template <
typename C>
1167 -> decltype(std::declval<C>().vmult_add(*r, *
d),
1168 std::declval<C>().Tvmult_add(*
d, *r),
1175 using type = decltype(test<T>(
nullptr,
nullptr));
1181 template <
typename Function,
typename Range,
typename Domain>
1204 template <
typename Range,
typename Domain,
typename Payload>
1208 template <
typename Matrix>
1269 template <
typename Range,
typename Domain,
typename Payload>
1273 template <
typename Matrix>
1376 template <
typename Range,
typename Domain,
typename Payload,
typename Matrix>
1381 return linear_operator<Range, Domain, Payload, Matrix, Matrix>(
matrix,
1400 template <
typename Range,
1403 typename OperatorExemplar,
1418 return_op.reinit_range_vector =
1419 [&operator_exemplar](Range &v,
bool omit_zeroing_entries) {
1421 Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1424 return_op.reinit_domain_vector = [&operator_exemplar](
1425 Domain &v,
bool omit_zeroing_entries) {
1427 Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1430 typename std::conditional<
1435 operator()(return_op,
matrix);
1458 template <
typename Range,
typename Domain,
typename Payload,
typename Matrix>
1465 auto return_op = operator_exemplar;
1467 typename std::conditional<
1472 operator()(return_op,
matrix);
1489 typename Domain = Range,
1491 typename OperatorExemplar,
1493 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1499 typename Domain = Range,
1501 typename OperatorExemplar,
1503 typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
1504 typename = std::enable_if_t<
1505 !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
1511 typename Domain = Range,
1513 typename OperatorExemplar,
1515 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>,
1516 typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
1517 typename = std::enable_if_t<
1518 !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
1524 typename Domain = Range,
1527 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1530 Matrix &&) =
delete;
1534 typename Domain = Range,
1537 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1544 typename Preconditioner,
1545 typename Range =
typename Solver::vector_type,
1546 typename Domain = Range,
1547 typename = std::enable_if_t<!std::is_lvalue_reference_v<Preconditioner>>,
1549 std::enable_if_t<!std::is_same_v<Preconditioner, PreconditionIdentity>>,
1550 typename = std::enable_if_t<
1551 !std::is_same_v<Preconditioner, LinearOperator<Range, Domain, Payload>>>>
1555 Preconditioner &&) =
delete;
LinearOperator< Range, Domain, Payload > & operator*=(const LinearOperator< Domain, Domain, Payload > &second_op)
LinearOperator< Range, Domain, Payload > operator*=(typename Domain::value_type number)
LinearOperator< Range, Domain, Payload > & operator-=(const LinearOperator< Range, Domain, Payload > &second_op)
std::function< void(Range &v, const Domain &u)> vmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
LinearOperator(const Payload &payload=Payload())
std::function< void(Range &v, const Domain &u)> vmult
LinearOperator< Range, Domain, Payload > & operator+=(const LinearOperator< Range, Domain, Payload > &second_op)
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
LinearOperator(const Op &op)
LinearOperator< Range, Domain, Payload > & operator=(const Op &op)
LinearOperator< Range, Domain, Payload > & operator=(const LinearOperator< Range, Domain, Payload > &)=default
LinearOperator(const LinearOperator< Range, Domain, Payload > &)=default
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
EmptyPayload(const Args &...)
EmptyPayload transpose_payload() const
EmptyPayload identity_payload() const
EmptyPayload null_payload() const
EmptyPayload inverse_payload(Solver &, const Preconditioner &) const
void operator()(LinearOperator< Range, Domain, Payload > &op, const Matrix &matrix)
void operator()(LinearOperator< Range, Domain, Payload > &op, const Matrix &matrix)
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static std::false_type test(...)
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())
decltype(test< T >(nullptr, nullptr)) type
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver)
LinearOperator< Range, Domain, Payload > operator*(const LinearOperator< Range, Intermediate, Payload > &first_op, const LinearOperator< Intermediate, Domain, Payload > &second_op)
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
LinearOperator< Range, Domain, Payload > operator-(const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
LinearOperator< Range, Range, Payload > identity_operator(const std::function< void(Range &, bool)> &reinit_vector)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > identity_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const PreconditionIdentity &)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const LinearOperator< Range, Domain, Payload > &preconditioner)
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > mean_value_filter(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > operator*(typename Range::value_type number, const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
LinearOperator< Range, Domain, Payload > operator*(const LinearOperator< Range, Domain, Payload > &op, typename Domain::value_type number)
LinearOperator< Range, Domain, Payload > linear_operator(const LinearOperator< Range, Domain, Payload > &operator_exemplar, const Matrix &matrix)
LinearOperator< Range, Range, Payload > mean_value_filter(const std::function< void(Range &, bool)> &reinit_vector)
LinearOperator< Range, Domain, Payload > operator+(const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
EmptyPayload operator+(const EmptyPayload &, const EmptyPayload &)
void apply_with_intermediate_storage(Function function, Range &v, const Domain &u, bool add)
EmptyPayload operator*(const EmptyPayload &, const EmptyPayload &)
static bool equal(const T *p1, const T *p2)