Reference documentation for deal.II version GIT c00406db16 2023-09-28 18:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
linear_operator.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2014 - 2022 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_linear_operator_h
17 #define dealii_linear_operator_h
18 
19 #include <deal.II/base/config.h>
20 
22 
24 
25 #include <array>
26 #include <functional>
27 #include <type_traits>
28 
30 
31 // Forward declarations:
32 #ifndef DOXYGEN
33 namespace internal
34 {
35  namespace LinearOperatorImplementation
36  {
37  class EmptyPayload;
38  }
39 } // namespace internal
40 
41 template <typename Number>
42 class Vector;
43 
45 
46 template <typename Range = Vector<double>,
47  typename Domain = Range,
48  typename Payload =
49  internal::LinearOperatorImplementation::EmptyPayload>
50 class LinearOperator;
51 #endif
52 
53 template <
54  typename Range = Vector<double>,
55  typename Domain = Range,
57  typename OperatorExemplar,
58  typename Matrix>
60 linear_operator(const OperatorExemplar &, const Matrix &);
61 
62 template <
63  typename Range = Vector<double>,
64  typename Domain = Range,
66  typename Matrix>
68 linear_operator(const Matrix &);
69 
70 template <
71  typename Range = Vector<double>,
72  typename Domain = Range,
76 
77 template <typename Range, typename Domain, typename Payload>
80 
81 
178 template <typename Range, typename Domain, typename Payload>
179 class LinearOperator : public Payload
180 {
181 public:
190  LinearOperator(const Payload &payload = Payload())
191  : Payload(payload)
192  , is_null_operator(false)
193  {
194  vmult = [](Range &, const Domain &) {
195  Assert(false,
196  ExcMessage("Uninitialized LinearOperator<Range, "
197  "Domain>::vmult called"));
198  };
199 
200  vmult_add = [](Range &, const Domain &) {
201  Assert(false,
202  ExcMessage("Uninitialized LinearOperator<Range, "
203  "Domain>::vmult_add called"));
204  };
205 
206  Tvmult = [](Domain &, const Range &) {
207  Assert(false,
208  ExcMessage("Uninitialized LinearOperator<Range, "
209  "Domain>::Tvmult called"));
210  };
211 
212  Tvmult_add = [](Domain &, const Range &) {
213  Assert(false,
214  ExcMessage("Uninitialized LinearOperator<Range, "
215  "Domain>::Tvmult_add called"));
216  };
217 
218  reinit_range_vector = [](Range &, bool) {
219  Assert(false,
220  ExcMessage("Uninitialized LinearOperator<Range, "
221  "Domain>::reinit_range_vector method called"));
222  };
223 
224  reinit_domain_vector = [](Domain &, bool) {
225  Assert(false,
226  ExcMessage("Uninitialized LinearOperator<Range, "
227  "Domain>::reinit_domain_vector method called"));
228  };
229  }
230 
235 
241  template <typename Op,
242  typename = std::enable_if_t<
243  !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
244  LinearOperator(const Op &op)
245  {
246  *this = linear_operator<Range, Domain, Payload, Op>(op);
247  }
248 
254 
259  template <typename Op,
260  typename = std::enable_if_t<
261  !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
263  operator=(const Op &op)
264  {
265  *this = linear_operator<Range, Domain, Payload, Op>(op);
266  return *this;
267  }
268 
273  std::function<void(Range &v, const Domain &u)> vmult;
274 
279  std::function<void(Range &v, const Domain &u)> vmult_add;
280 
285  std::function<void(Domain &v, const Range &u)> Tvmult;
286 
291  std::function<void(Domain &v, const Range &u)> Tvmult_add;
292 
300  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_range_vector;
301 
309  std::function<void(Domain &v, bool omit_zeroing_entries)>
311 
323  {
324  *this = *this + second_op;
325  return *this;
326  }
327 
334  {
335  *this = *this - second_op;
336  return *this;
337  }
338 
345  {
346  *this = *this * second_op;
347  return *this;
348  }
349 
355  operator*=(typename Domain::value_type number)
356  {
357  *this = *this * number;
358  return *this;
359  }
360 
367 
369 };
370 
371 
386 template <typename Range, typename Domain, typename Payload>
390 {
391  if (first_op.is_null_operator)
392  {
393  return second_op;
394  }
395  else if (second_op.is_null_operator)
396  {
397  return first_op;
398  }
399  else
400  {
402  static_cast<const Payload &>(first_op) +
403  static_cast<const Payload &>(second_op)};
404 
405  return_op.reinit_range_vector = first_op.reinit_range_vector;
406  return_op.reinit_domain_vector = first_op.reinit_domain_vector;
407 
408  // ensure to have valid computation objects by catching first_op and
409  // second_op by value
410 
411  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
412  first_op.vmult(v, u);
413  second_op.vmult_add(v, u);
414  };
415 
416  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
417  first_op.vmult_add(v, u);
418  second_op.vmult_add(v, u);
419  };
420 
421  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
422  second_op.Tvmult(v, u);
423  first_op.Tvmult_add(v, u);
424  };
425 
426  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
427  second_op.Tvmult_add(v, u);
428  first_op.Tvmult_add(v, u);
429  };
430 
431  return return_op;
432  }
433 }
434 
435 
445 template <typename Range, typename Domain, typename Payload>
449 {
450  if (first_op.is_null_operator)
451  {
452  return -1. * second_op;
453  }
454  else if (second_op.is_null_operator)
455  {
456  return first_op;
457  }
458  else
459  {
460  // implement with addition and scalar multiplication
461  return first_op + (-1. * second_op);
462  }
463 }
464 
465 
483 template <typename Range, typename Domain, typename Payload>
485 operator*(typename Range::value_type number,
487 {
488  static_assert(
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");
492 
493  if (op.is_null_operator)
494  {
495  return op;
496  }
497  else if (number == 0.)
498  {
499  return null_operator(op);
500  }
501  else
502  {
504 
505  // ensure to have valid computation objects by catching number and op by
506  // value
507 
508  return_op.vmult = [number, op](Range &v, const Domain &u) {
509  op.vmult(v, u);
510  v *= number;
511  };
512 
513  return_op.vmult_add = [number, op](Range &v, const Domain &u) {
514  v /= number;
515  op.vmult_add(v, u);
516  v *= number;
517  };
518 
519  return_op.Tvmult = [number, op](Domain &v, const Range &u) {
520  op.Tvmult(v, u);
521  v *= number;
522  };
523 
524  return_op.Tvmult_add = [number, op](Domain &v, const Range &u) {
525  v /= number;
526  op.Tvmult_add(v, u);
527  v *= number;
528  };
529 
530  return return_op;
531  }
532 }
533 
534 
550 template <typename Range, typename Domain, typename Payload>
553  typename Domain::value_type number)
554 {
555  static_assert(
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");
559 
560  return number * op;
561 }
562 
580 template <typename Range,
581  typename Intermediate,
582  typename Domain,
583  typename Payload>
587 {
588  if (first_op.is_null_operator || second_op.is_null_operator)
589  {
591  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
592  return_op.reinit_range_vector = first_op.reinit_range_vector;
593  return null_operator(return_op);
594  }
595  else
596  {
598  static_cast<const Payload &>(first_op) *
599  static_cast<const Payload &>(second_op)};
600 
601  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
602  return_op.reinit_range_vector = first_op.reinit_range_vector;
603 
604  // ensure to have valid computation objects by catching first_op and
605  // second_op by value
606 
607  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
608  GrowingVectorMemory<Intermediate> vector_memory;
609 
610  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
611  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
612  second_op.vmult(*i, u);
613  first_op.vmult(v, *i);
614  };
615 
616  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
617  GrowingVectorMemory<Intermediate> vector_memory;
618 
619  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
620  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
621  second_op.vmult(*i, u);
622  first_op.vmult_add(v, *i);
623  };
624 
625  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
626  GrowingVectorMemory<Intermediate> vector_memory;
627 
628  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
629  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
630  first_op.Tvmult(*i, u);
631  second_op.Tvmult(v, *i);
632  };
633 
634  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
635  GrowingVectorMemory<Intermediate> vector_memory;
636 
637  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
638  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
639  first_op.Tvmult(*i, u);
640  second_op.Tvmult_add(v, *i);
641  };
642 
643  return return_op;
644  }
645 }
646 
647 
656 template <typename Range, typename Domain, typename Payload>
659 {
660  LinearOperator<Domain, Range, Payload> return_op{op.transpose_payload()};
661 
663  return_op.reinit_domain_vector = op.reinit_range_vector;
664 
665  return_op.vmult = op.Tvmult;
666  return_op.vmult_add = op.Tvmult_add;
667  return_op.Tvmult = op.vmult;
668  return_op.Tvmult_add = op.vmult_add;
669 
670  return return_op;
671 }
672 
673 
693 template <typename Payload,
694  typename Solver,
695  typename Preconditioner,
696  typename Range = typename Solver::vector_type,
697  typename Domain = Range>
700  Solver &solver,
701  const Preconditioner &preconditioner)
702 {
704  op.inverse_payload(solver, preconditioner)};
705 
707  return_op.reinit_domain_vector = op.reinit_range_vector;
708 
709  return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
710  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
711  solver.solve(op, v, u, preconditioner);
712  };
713 
714  return_op.vmult_add = [op, &solver, &preconditioner](Range &v,
715  const Domain &u) {
716  GrowingVectorMemory<Range> vector_memory;
717 
718  typename VectorMemory<Range>::Pointer v2(vector_memory);
719  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
720  solver.solve(op, *v2, u, preconditioner);
721  v += *v2;
722  };
723 
724  return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
725  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
726  solver.solve(transpose_operator(op), v, u, preconditioner);
727  };
728 
729  return_op.Tvmult_add = [op, &solver, &preconditioner](Range &v,
730  const Domain &u) {
731  GrowingVectorMemory<Range> vector_memory;
732 
733  typename VectorMemory<Range>::Pointer v2(vector_memory);
734  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
735  solver.solve(transpose_operator(op), *v2, u, preconditioner);
736  v += *v2;
737  };
738 
739  return return_op;
740 }
741 
742 
751 template <typename Payload,
752  typename Solver,
753  typename Range = typename Solver::vector_type,
754  typename Domain = Range>
757  Solver &solver,
758  const LinearOperator<Range, Domain, Payload> &preconditioner)
759 {
761  op.inverse_payload(solver, preconditioner)};
762 
764  return_op.reinit_domain_vector = op.reinit_range_vector;
765 
766  return_op.vmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
767  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
768  solver.solve(op, v, u, preconditioner);
769  };
770 
771  return_op.vmult_add = [op, &solver, preconditioner](Range &v,
772  const Domain &u) {
773  GrowingVectorMemory<Range> vector_memory;
774 
775  typename VectorMemory<Range>::Pointer v2(vector_memory);
776  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
777  solver.solve(op, *v2, u, preconditioner);
778  v += *v2;
779  };
780 
781  return_op.Tvmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
782  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
783  solver.solve(transpose_operator(op), v, u, preconditioner);
784  };
785 
786  return_op.Tvmult_add = [op, &solver, preconditioner](Range &v,
787  const Domain &u) {
788  GrowingVectorMemory<Range> vector_memory;
789 
790  typename VectorMemory<Range>::Pointer v2(vector_memory);
791  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
792  solver.solve(transpose_operator(op), *v2, u, preconditioner);
793  v += *v2;
794  };
795 
796  return return_op;
797 }
798 
799 
809 template <typename Payload,
810  typename Solver,
811  typename Range = typename Solver::vector_type,
812  typename Domain = Range>
815  Solver &solver)
816 {
817  return inverse_operator(op, solver, identity_operator(op));
818 }
819 
820 
829 template <typename Payload,
830  typename Solver,
831  typename Range = typename Solver::vector_type,
832  typename Domain = Range>
835  Solver &solver,
836  const PreconditionIdentity &)
837 {
838  return inverse_operator(op, solver);
839 }
840 
860 template <
861  typename Range,
864 identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
865 {
867 
868  return_op.reinit_range_vector = reinit_vector;
869  return_op.reinit_domain_vector = reinit_vector;
870 
871  return_op.vmult = [](Range &v, const Range &u) { v = u; };
872 
873  return_op.vmult_add = [](Range &v, const Range &u) { v += u; };
874 
875  return_op.Tvmult = [](Range &v, const Range &u) { v = u; };
876 
877  return_op.Tvmult_add = [](Range &v, const Range &u) { v += u; };
878 
879  return return_op;
880 }
881 
882 
895 template <typename Range, typename Domain, typename Payload>
898 {
899  auto return_op = identity_operator<Range, Payload>(op.reinit_range_vector);
900  static_cast<Payload &>(return_op) = op.identity_payload();
901 
902  return return_op;
903 }
904 
905 
915 template <typename Range, typename Domain, typename Payload>
918 {
919  LinearOperator<Range, Domain, Payload> return_op{op.null_payload()};
920 
921  return_op.is_null_operator = true;
922 
923  return_op.reinit_range_vector = op.reinit_range_vector;
924  return_op.reinit_domain_vector = op.reinit_domain_vector;
925 
926  return_op.vmult = [](Range &v, const Domain &) { v = 0.; };
927 
928  return_op.vmult_add = [](Range &, const Domain &) {};
929 
930  return_op.Tvmult = [](Domain &v, const Range &) { v = 0.; };
931 
932  return_op.Tvmult_add = [](Domain &, const Range &) {};
933 
934  return return_op;
935 }
936 
937 
950 template <
951  typename Range,
954 mean_value_filter(const std::function<void(Range &, bool)> &reinit_vector)
955 {
957 
958  return_op.reinit_range_vector = reinit_vector;
959  return_op.reinit_domain_vector = reinit_vector;
960 
961  return_op.vmult = [](Range &v, const Range &u) {
962  const auto mean = u.mean_value();
963 
964  v = u;
965  v.add(-mean);
966  };
967 
968  return_op.vmult_add = [](Range &v, const Range &u) {
969  const auto mean = u.mean_value();
970 
971  v += u;
972  v.add(-mean);
973  };
974 
975  return_op.Tvmult = return_op.vmult_add;
976  return_op.Tvmult_add = return_op.vmult_add;
977 
978  return return_op;
979 }
980 
981 
994 template <typename Range, typename Domain, typename Payload>
997 {
998  auto return_op = mean_value_filter<Range, Payload>(op.reinit_range_vector);
999  static_cast<Payload &>(return_op) = op.identity_payload();
1000 
1001  return return_op;
1002 }
1003 
1004 
1005 namespace internal
1006 {
1007  namespace LinearOperatorImplementation
1008  {
1020  template <typename Vector>
1022  {
1023  public:
1035  template <typename Matrix>
1036  static void
1038  Vector &v,
1039  bool omit_zeroing_entries)
1040  {
1041  v.reinit(matrix.m(), omit_zeroing_entries);
1042  }
1043 
1055  template <typename Matrix>
1056  static void
1058  Vector &v,
1059  bool omit_zeroing_entries)
1060  {
1061  v.reinit(matrix.n(), omit_zeroing_entries);
1062  }
1063  };
1064 
1065 
1078  {
1079  public:
1087  template <typename... Args>
1088  EmptyPayload(const Args &...)
1089  {}
1090 
1091 
1095  EmptyPayload
1097  {
1098  return *this;
1099  }
1100 
1101 
1105  EmptyPayload
1107  {
1108  return *this;
1109  }
1110 
1111 
1115  EmptyPayload
1117  {
1118  return *this;
1119  }
1120 
1121 
1125  template <typename Solver, typename Preconditioner>
1126  EmptyPayload
1127  inverse_payload(Solver &, const Preconditioner &) const
1128  {
1129  return *this;
1130  }
1131  };
1132 
1137  inline EmptyPayload
1139  {
1140  return {};
1141  }
1142 
1147  inline EmptyPayload
1149  {
1150  return {};
1151  }
1152 
1153 
1154 
1155  // A trait class that determines whether type T provides public
1156  // (templated or non-templated) vmult_add member functions
1157  template <typename Range, typename Domain, typename T>
1159  {
1160  template <typename C>
1161  static std::false_type
1162  test(...);
1163 
1164  template <typename C>
1165  static auto
1166  test(Range *r, Domain *d)
1167  -> decltype(std::declval<C>().vmult_add(*r, *d),
1168  std::declval<C>().Tvmult_add(*d, *r),
1169  std::true_type());
1170 
1171  public:
1172  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
1173  // otherwise it is std::false_type
1174 
1175  using type = decltype(test<T>(nullptr, nullptr));
1176  };
1177 
1178 
1179  // A helper function to apply a given vmult, or Tvmult to a vector with
1180  // intermediate storage
1181  template <typename Function, typename Range, typename Domain>
1182  void
1184  Range &v,
1185  const Domain &u,
1186  bool add)
1187  {
1188  GrowingVectorMemory<Range> vector_memory;
1189 
1190  typename VectorMemory<Range>::Pointer i(vector_memory);
1191  i->reinit(v, /*bool omit_zeroing_entries =*/true);
1192 
1193  function(*i, u);
1194 
1195  if (add)
1196  v += *i;
1197  else
1198  v = *i;
1199  }
1200 
1201 
1202  // A helper class to add a reduced matrix interface to a LinearOperator
1203  // (typically provided by Preconditioner classes)
1204  template <typename Range, typename Domain, typename Payload>
1206  {
1207  public:
1208  template <typename Matrix>
1209  void
1211  const Matrix &matrix)
1212  {
1213  op.vmult = [&matrix](Range &v, const Domain &u) {
1214  if (PointerComparison::equal(&v, &u))
1215  {
1216  // If v and u are the same memory location use intermediate
1217  // storage
1219  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1220  v,
1221  u,
1222  /*bool add =*/false);
1223  }
1224  else
1225  {
1226  matrix.vmult(v, u);
1227  }
1228  };
1229 
1230  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1231  // use intermediate storage to implement vmult_add with vmult
1233  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1234  v,
1235  u,
1236  /*bool add =*/true);
1237  };
1238 
1239  op.Tvmult = [&matrix](Domain &v, const Range &u) {
1240  if (PointerComparison::equal(&v, &u))
1241  {
1242  // If v and u are the same memory location use intermediate
1243  // storage
1245  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1246  v,
1247  u,
1248  /*bool add =*/false);
1249  }
1250  else
1251  {
1252  matrix.Tvmult(v, u);
1253  }
1254  };
1255 
1256  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1257  // use intermediate storage to implement Tvmult_add with Tvmult
1259  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1260  v,
1261  u,
1262  /*bool add =*/true);
1263  };
1264  }
1265  };
1266 
1267 
1268  // A helper class to add the full matrix interface to a LinearOperator
1269  template <typename Range, typename Domain, typename Payload>
1271  {
1272  public:
1273  template <typename Matrix>
1274  void
1276  const Matrix &matrix)
1277  {
1278  // As above ...
1279 
1281  op, matrix);
1282 
1283  // ... but add native vmult_add and Tvmult_add variants:
1284 
1285  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1286  if (PointerComparison::equal(&v, &u))
1287  {
1289  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1290  v,
1291  u,
1292  /*bool add =*/true);
1293  }
1294  else
1295  {
1296  matrix.vmult_add(v, u);
1297  }
1298  };
1299 
1300  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1301  if (PointerComparison::equal(&v, &u))
1302  {
1304  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1305  v,
1306  u,
1307  /*bool add =*/true);
1308  }
1309  else
1310  {
1311  matrix.Tvmult_add(v, u);
1312  }
1313  };
1314  }
1315  };
1316  } // namespace LinearOperatorImplementation
1317 } // namespace internal
1318 
1319 
1376 template <typename Range, typename Domain, typename Payload, typename Matrix>
1378 linear_operator(const Matrix &matrix)
1379 {
1380  // implement with the more generic variant below...
1381  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1382  matrix);
1383 }
1384 
1385 
1400 template <typename Range,
1401  typename Domain,
1402  typename Payload,
1403  typename OperatorExemplar,
1404  typename Matrix>
1406 linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
1407 {
1409  // Initialize the payload based on the input exemplar matrix
1411  Payload(operator_exemplar, matrix)};
1412 
1413  // Always store a reference to matrix and operator_exemplar in the lambda
1414  // functions. This ensures that a modification of the matrix after the
1415  // creation of a LinearOperator wrapper is respected - further a matrix
1416  // or an operator_exemplar cannot usually be copied...
1417 
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);
1422  };
1423 
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);
1428  };
1429 
1430  typename std::conditional<
1434  .
1435  operator()(return_op, matrix);
1436 
1437  return return_op;
1438 }
1439 
1440 
1441 
1458 template <typename Range, typename Domain, typename Payload, typename Matrix>
1461  const Matrix &matrix)
1462 {
1464  // Initialize the payload based on the LinearOperator exemplar
1465  auto return_op = operator_exemplar;
1466 
1467  typename std::conditional<
1471  .
1472  operator()(return_op, matrix);
1473 
1474  return return_op;
1475 }
1476 
1477 
1480 #ifndef DOXYGEN
1481 
1482 //
1483 // Ensure that we never capture a reference to a temporary by accident.
1484 // to avoid "stack use after free".
1485 //
1486 
1487 template <
1488  typename Range = Vector<double>,
1489  typename Domain = Range,
1491  typename OperatorExemplar,
1492  typename Matrix,
1493  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1495 linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1496 
1497 template <
1498  typename Range = Vector<double>,
1499  typename Domain = Range,
1501  typename OperatorExemplar,
1502  typename Matrix,
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>>>>
1507 linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1508 
1509 template <
1510  typename Range = Vector<double>,
1511  typename Domain = Range,
1513  typename OperatorExemplar,
1514  typename Matrix,
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>>>>
1520 linear_operator(OperatorExemplar &&, Matrix &&) = delete;
1521 
1522 template <
1523  typename Range = Vector<double>,
1524  typename Domain = Range,
1526  typename Matrix,
1527  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1530  Matrix &&) = delete;
1531 
1532 template <
1533  typename Range = Vector<double>,
1534  typename Domain = Range,
1536  typename Matrix,
1537  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1539 linear_operator(Matrix &&) = delete;
1540 
1541 template <
1542  typename Payload,
1543  typename Solver,
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>>,
1548  typename =
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>>>>
1554  Solver &,
1555  Preconditioner &&) = delete;
1556 
1557 #endif // DOXYGEN
1558 
1560 
1561 #endif
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
Definition: vector.h:110
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
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 auto test(Range *r, Domain *d) -> decltype(std::declval< C >().vmult_add(*r, *d), std::declval< C >().Tvmult_add(*d, *r), std::true_type())
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1616
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)