16 #ifndef dealii_linear_operator_h 17 #define dealii_linear_operator_h 19 #include <deal.II/base/config.h> 21 #include <deal.II/base/exceptions.h> 23 #include <deal.II/lac/vector_memory.h> 27 #include <type_traits> 29 DEAL_II_NAMESPACE_OPEN
35 namespace LinearOperatorImplementation
41 template <
typename Number>
44 template <
typename Range = Vector<
double>,
45 typename Domain = Range,
47 internal::LinearOperatorImplementation::EmptyPayload>
52 typename Domain = Range,
54 typename OperatorExemplar,
61 typename Domain = Range,
69 typename Domain = Range,
158 template <
typename Range,
typename Domain,
typename Payload>
172 , is_null_operator(false)
174 vmult = [](Range &,
const Domain &) {
176 ExcMessage(
"Uninitialized LinearOperator<Range, " 177 "Domain>::vmult called"));
180 vmult_add = [](Range &,
const Domain &) {
182 ExcMessage(
"Uninitialized LinearOperator<Range, " 183 "Domain>::vmult_add called"));
186 Tvmult = [](Domain &,
const Range &) {
188 ExcMessage(
"Uninitialized LinearOperator<Range, " 189 "Domain>::Tvmult called"));
192 Tvmult_add = [](Domain &,
const Range &) {
194 ExcMessage(
"Uninitialized LinearOperator<Range, " 195 "Domain>::Tvmult_add called"));
198 reinit_range_vector = [](Range &, bool) {
200 ExcMessage(
"Uninitialized LinearOperator<Range, " 201 "Domain>::reinit_range_vector method called"));
204 reinit_domain_vector = [](Domain &, bool) {
206 ExcMessage(
"Uninitialized LinearOperator<Range, " 207 "Domain>::reinit_domain_vector method called"));
221 template <
typename Op,
222 typename =
typename std::enable_if<
223 !std::is_base_of<LinearOperator<Range, Domain, Payload>,
227 *
this = linear_operator<Range, Domain, Payload, Op>(op);
240 template <
typename Op,
241 typename =
typename std::enable_if<
242 !std::is_base_of<LinearOperator<Range, Domain, Payload>,
247 *
this = linear_operator<Range, Domain, Payload, Op>(op);
255 std::function<void(Range &v, const Domain &u)>
vmult;
261 std::function<void(Range &v, const Domain &u)>
vmult_add;
267 std::function<void(Domain &v, const Range &u)>
Tvmult;
291 std::function<void(Domain &v, bool omit_zeroing_entries)>
306 *
this = *
this + second_op;
317 *
this = *
this - second_op;
328 *
this = *
this * second_op;
339 *
this = *
this * number;
368 template <
typename Range,
typename Domain,
typename Payload>
384 static_cast<const Payload &>(first_op) +
385 static_cast<const Payload &>(second_op));
393 return_op.
vmult = [first_op, second_op](Range &v,
const Domain &u) {
394 first_op.
vmult(v, u);
398 return_op.
vmult_add = [first_op, second_op](Range &v,
const Domain &u) {
403 return_op.
Tvmult = [first_op, second_op](Domain &v,
const Range &u) {
408 return_op.
Tvmult_add = [first_op, second_op](Domain &v,
const Range &u) {
427 template <
typename Range,
typename Domain,
typename Payload>
434 return -1. * second_op;
443 return first_op + (-1. * second_op);
465 template <
typename Range,
typename Domain,
typename Payload>
471 std::is_convertible<
typename Range::value_type,
472 typename Domain::value_type>::value,
473 "Range and Domain must have implicitly convertible 'value_type's");
479 else if (number == 0.)
490 return_op.
vmult = [number, op](Range &v,
const Domain &u) {
495 return_op.
vmult_add = [number, op](Range &v,
const Domain &u) {
501 return_op.
Tvmult = [number, op](Domain &v,
const Range &u) {
506 return_op.
Tvmult_add = [number, op](Domain &v,
const Range &u) {
532 template <
typename Range,
typename Domain,
typename Payload>
535 typename Domain::value_type number)
538 std::is_convertible<
typename Range::value_type,
539 typename Domain::value_type>::value,
540 "Range and Domain must have implicitly convertible 'value_type's");
562 template <
typename Range,
563 typename Intermediate,
580 static_cast<const Payload &>(first_op) *
581 static_cast<const Payload &>(second_op));
589 return_op.
vmult = [first_op, second_op](Range &v,
const Domain &u) {
594 second_op.
vmult(*i, u);
595 first_op.
vmult(v, *i);
598 return_op.
vmult_add = [first_op, second_op](Range &v,
const Domain &u) {
603 second_op.
vmult(*i, u);
607 return_op.
Tvmult = [first_op, second_op](Domain &v,
const Range &u) {
616 return_op.
Tvmult_add = [first_op, second_op](Domain &v,
const Range &u) {
638 template <
typename Range,
typename Domain,
typename Payload>
647 return_op.vmult = op.
Tvmult;
649 return_op.Tvmult = op.
vmult;
676 template <
typename Payload,
678 typename Preconditioner,
680 typename Domain = Range>
684 const Preconditioner & preconditioner)
687 op.inverse_payload(solver, preconditioner));
692 return_op.vmult = [op, &solver, &preconditioner](Range &v,
const Domain &u) {
694 solver.solve(op, v, u, preconditioner);
697 return_op.vmult_add = [op, &solver, &preconditioner](Range & v,
703 solver.solve(op, *v2, u, preconditioner);
707 return_op.Tvmult = [op, &solver, &preconditioner](Range &v,
const Domain &u) {
712 return_op.Tvmult_add = [op, &solver, &preconditioner](Range & v,
755 return_op.
vmult = [](Range &v,
const Range &u) { v = u; };
757 return_op.
vmult_add = [](Range &v,
const Range &u) { v += u; };
759 return_op.
Tvmult = [](Range &v,
const Range &u) { v = u; };
761 return_op.
Tvmult_add = [](Range &v,
const Range &u) { v += u; };
779 template <
typename Range,
typename Domain,
typename Payload>
784 static_cast<Payload &
>(return_op) = op.identity_payload();
799 template <
typename Range,
typename Domain,
typename Payload>
810 return_op.vmult = [](Range &v,
const Domain &) { v = 0.; };
812 return_op.vmult_add = [](Range &,
const Domain &) {};
814 return_op.Tvmult = [](Domain &v,
const Range &) { v = 0.; };
816 return_op.Tvmult_add = [](Domain &,
const Range &) {};
845 return_op.
vmult = [](Range &v,
const Range &u) {
846 const auto mean = u.mean_value();
852 return_op.
vmult_add = [](Range &v,
const Range &u) {
853 const auto mean = u.mean_value();
878 template <
typename Range,
typename Domain,
typename Payload>
883 static_cast<Payload &
>(return_op) = op.identity_payload();
891 namespace LinearOperatorImplementation
904 template <
typename Vector>
919 template <
typename Matrix>
923 bool omit_zeroing_entries)
925 v.reinit(matrix.m(), omit_zeroing_entries);
939 template <
typename Matrix>
943 bool omit_zeroing_entries)
945 v.reinit(matrix.n(), omit_zeroing_entries);
972 template <
typename... Args>
1010 template <
typename Solver,
typename Preconditioner>
1041 template <
typename Range,
typename Domain,
typename T>
1042 class has_vmult_add_and_Tvmult_add
1044 template <
typename C>
1045 static std::false_type
1048 template <
typename C>
1050 test(Range *r, Domain *d)
1051 -> decltype(std::declval<C>().vmult_add(*r, *d),
1052 std::declval<C>().Tvmult_add(*d, *r),
1059 using type = decltype(test<T>(
nullptr,
nullptr));
1065 template <
typename Function,
typename Range,
typename Domain>
1067 apply_with_intermediate_storage(
Function function,
1088 template <
typename Range,
typename Domain,
typename Payload>
1089 class MatrixInterfaceWithoutVmultAdd
1092 template <
typename Matrix>
1095 const Matrix & matrix)
1097 op.
vmult = [&matrix](Range &v,
const Domain &u) {
1102 apply_with_intermediate_storage(
1103 [&matrix](Range &b,
const Domain &a) { matrix.vmult(b, a); },
1114 op.
vmult_add = [&matrix](Range &v,
const Domain &u) {
1116 apply_with_intermediate_storage(
1117 [&matrix](Range &b,
const Domain &a) { matrix.vmult(b, a); },
1123 op.
Tvmult = [&matrix](Domain &v,
const Range &u) {
1128 apply_with_intermediate_storage(
1129 [&matrix](Domain &b,
const Range &a) { matrix.Tvmult(b, a); },
1136 matrix.Tvmult(v, u);
1140 op.
Tvmult_add = [&matrix](Domain &v,
const Range &u) {
1142 apply_with_intermediate_storage(
1143 [&matrix](Domain &b,
const Range &a) { matrix.Tvmult(b, a); },
1153 template <
typename Range,
typename Domain,
typename Payload>
1154 class MatrixInterfaceWithVmultAdd
1157 template <
typename Matrix>
1160 const Matrix & matrix)
1164 MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>().
operator()(
1169 op.
vmult_add = [&matrix](Range &v,
const Domain &u) {
1172 apply_with_intermediate_storage(
1173 [&matrix](Range &b,
const Domain &a) { matrix.vmult(b, a); },
1180 matrix.vmult_add(v, u);
1184 op.
Tvmult_add = [&matrix](Domain &v,
const Range &u) {
1187 apply_with_intermediate_storage(
1188 [&matrix](Domain &b,
const Range &a) { matrix.Tvmult(b, a); },
1195 matrix.Tvmult_add(v, u);
1261 template <
typename Range,
typename Domain,
typename Payload,
typename Matrix>
1266 return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1286 template <
typename Range,
1289 typename OperatorExemplar,
1297 Payload(operator_exemplar, matrix));
1305 [&operator_exemplar](Range &v,
bool omit_zeroing_entries) {
1307 Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1311 Domain &v,
bool omit_zeroing_entries) {
1313 Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1316 typename std::conditional<
1317 has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
1318 MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
1319 MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>::type()
1321 operator()(return_op, matrix);
1345 template <
typename Range,
typename Domain,
typename Payload,
typename Matrix>
1348 const Matrix & matrix)
1352 auto return_op = operator_exemplar;
1354 typename std::conditional<
1355 has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
1356 MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
1357 MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>::type()
1359 operator()(return_op, matrix);
1368 DEAL_II_NAMESPACE_CLOSE
EmptyPayload inverse_payload(Solver &, const Preconditioner &) const
EmptyPayload transpose_payload() const
LinearOperator< Range, Range, Payload > mean_value_filter(const std::function< void(Range &, bool)> &reinit_vector)
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
LinearOperator< Range, Domain, Payload > operator*=(typename Domain::value_type number)
LinearOperator< Range, Domain, Payload > operator+(const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
LinearOperator< Range, Domain, Payload > & operator*=(const LinearOperator< Domain, Domain, Payload > &second_op)
static void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &)
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
EmptyPayload null_payload() const
static::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
std::function< void(Range &v, const Domain &u)> vmult
LinearOperator< Range, Domain, Payload > & operator=(const Op &op)
LinearOperator< Range, Domain, Payload > & operator-=(const LinearOperator< Range, Domain, Payload > &second_op)
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
LinearOperator< Range, Domain, Payload > operator*(typename Range::value_type number, const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Range, Payload > identity_operator(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)
std::function< void(Range &v, const Domain &u)> vmult_add
EmptyPayload(const Args &...)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
EmptyPayload identity_payload() const
LinearOperator< Range, Domain, Payload > & operator+=(const LinearOperator< Range, Domain, Payload > &second_op)
LinearOperator(const Payload &payload=Payload())
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
static bool equal(const T *p1, const T *p2)
LinearOperator(const Op &op)