Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50: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 
200 template <typename Range, typename Domain, typename Payload>
201 class LinearOperator : public Payload
202 {
203 public:
212  LinearOperator(const Payload &payload = Payload())
213  : Payload(payload)
214  , is_null_operator(false)
215  {
216  vmult = [](Range &, const Domain &) {
217  Assert(false,
218  ExcMessage("Uninitialized LinearOperator<Range, "
219  "Domain>::vmult called"));
220  };
221 
222  vmult_add = [](Range &, const Domain &) {
223  Assert(false,
224  ExcMessage("Uninitialized LinearOperator<Range, "
225  "Domain>::vmult_add called"));
226  };
227 
228  Tvmult = [](Domain &, const Range &) {
229  Assert(false,
230  ExcMessage("Uninitialized LinearOperator<Range, "
231  "Domain>::Tvmult called"));
232  };
233 
234  Tvmult_add = [](Domain &, const Range &) {
235  Assert(false,
236  ExcMessage("Uninitialized LinearOperator<Range, "
237  "Domain>::Tvmult_add called"));
238  };
239 
240  reinit_range_vector = [](Range &, bool) {
241  Assert(false,
242  ExcMessage("Uninitialized LinearOperator<Range, "
243  "Domain>::reinit_range_vector method called"));
244  };
245 
246  reinit_domain_vector = [](Domain &, bool) {
247  Assert(false,
248  ExcMessage("Uninitialized LinearOperator<Range, "
249  "Domain>::reinit_domain_vector method called"));
250  };
251  }
252 
257 
263  template <typename Op,
264  typename = std::enable_if_t<
265  !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
266  LinearOperator(const Op &op)
267  {
268  *this = linear_operator<Range, Domain, Payload, Op>(op);
269  }
270 
276 
281  template <typename Op,
282  typename = std::enable_if_t<
283  !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
285  operator=(const Op &op)
286  {
287  *this = linear_operator<Range, Domain, Payload, Op>(op);
288  return *this;
289  }
290 
295  std::function<void(Range &v, const Domain &u)> vmult;
296 
301  std::function<void(Range &v, const Domain &u)> vmult_add;
302 
307  std::function<void(Domain &v, const Range &u)> Tvmult;
308 
313  std::function<void(Domain &v, const Range &u)> Tvmult_add;
314 
322  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_range_vector;
323 
331  std::function<void(Domain &v, bool omit_zeroing_entries)>
333 
345  {
346  *this = *this + second_op;
347  return *this;
348  }
349 
356  {
357  *this = *this - second_op;
358  return *this;
359  }
360 
367  {
368  *this = *this * second_op;
369  return *this;
370  }
371 
377  operator*=(typename Domain::value_type number)
378  {
379  *this = *this * number;
380  return *this;
381  }
382 
389 
391 };
392 
393 
408 template <typename Range, typename Domain, typename Payload>
412 {
413  if (first_op.is_null_operator)
414  {
415  return second_op;
416  }
417  else if (second_op.is_null_operator)
418  {
419  return first_op;
420  }
421  else
422  {
424  static_cast<const Payload &>(first_op) +
425  static_cast<const Payload &>(second_op)};
426 
427  return_op.reinit_range_vector = first_op.reinit_range_vector;
428  return_op.reinit_domain_vector = first_op.reinit_domain_vector;
429 
430  // ensure to have valid computation objects by catching first_op and
431  // second_op by value
432 
433  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
434  first_op.vmult(v, u);
435  second_op.vmult_add(v, u);
436  };
437 
438  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
439  first_op.vmult_add(v, u);
440  second_op.vmult_add(v, u);
441  };
442 
443  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
444  second_op.Tvmult(v, u);
445  first_op.Tvmult_add(v, u);
446  };
447 
448  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
449  second_op.Tvmult_add(v, u);
450  first_op.Tvmult_add(v, u);
451  };
452 
453  return return_op;
454  }
455 }
456 
457 
467 template <typename Range, typename Domain, typename Payload>
471 {
472  if (first_op.is_null_operator)
473  {
474  return -1. * second_op;
475  }
476  else if (second_op.is_null_operator)
477  {
478  return first_op;
479  }
480  else
481  {
482  // implement with addition and scalar multiplication
483  return first_op + (-1. * second_op);
484  }
485 }
486 
487 
505 template <typename Range, typename Domain, typename Payload>
507 operator*(typename Range::value_type number,
509 {
510  static_assert(
511  std::is_convertible<typename Range::value_type,
512  typename Domain::value_type>::value,
513  "Range and Domain must have implicitly convertible 'value_type's");
514 
515  if (op.is_null_operator)
516  {
517  return op;
518  }
519  else if (number == 0.)
520  {
521  return null_operator(op);
522  }
523  else
524  {
526 
527  // ensure to have valid computation objects by catching number and op by
528  // value
529 
530  return_op.vmult = [number, op](Range &v, const Domain &u) {
531  op.vmult(v, u);
532  v *= number;
533  };
534 
535  return_op.vmult_add = [number, op](Range &v, const Domain &u) {
536  v /= number;
537  op.vmult_add(v, u);
538  v *= number;
539  };
540 
541  return_op.Tvmult = [number, op](Domain &v, const Range &u) {
542  op.Tvmult(v, u);
543  v *= number;
544  };
545 
546  return_op.Tvmult_add = [number, op](Domain &v, const Range &u) {
547  v /= number;
548  op.Tvmult_add(v, u);
549  v *= number;
550  };
551 
552  return return_op;
553  }
554 }
555 
556 
572 template <typename Range, typename Domain, typename Payload>
575  typename Domain::value_type number)
576 {
577  static_assert(
578  std::is_convertible<typename Range::value_type,
579  typename Domain::value_type>::value,
580  "Range and Domain must have implicitly convertible 'value_type's");
581 
582  return number * op;
583 }
584 
602 template <typename Range,
603  typename Intermediate,
604  typename Domain,
605  typename Payload>
609 {
610  if (first_op.is_null_operator || second_op.is_null_operator)
611  {
613  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
614  return_op.reinit_range_vector = first_op.reinit_range_vector;
615  return null_operator(return_op);
616  }
617  else
618  {
620  static_cast<const Payload &>(first_op) *
621  static_cast<const Payload &>(second_op)};
622 
623  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
624  return_op.reinit_range_vector = first_op.reinit_range_vector;
625 
626  // ensure to have valid computation objects by catching first_op and
627  // second_op by value
628 
629  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
630  GrowingVectorMemory<Intermediate> vector_memory;
631 
632  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
633  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
634  second_op.vmult(*i, u);
635  first_op.vmult(v, *i);
636  };
637 
638  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
639  GrowingVectorMemory<Intermediate> vector_memory;
640 
641  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
642  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
643  second_op.vmult(*i, u);
644  first_op.vmult_add(v, *i);
645  };
646 
647  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
648  GrowingVectorMemory<Intermediate> vector_memory;
649 
650  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
651  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
652  first_op.Tvmult(*i, u);
653  second_op.Tvmult(v, *i);
654  };
655 
656  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
657  GrowingVectorMemory<Intermediate> vector_memory;
658 
659  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
660  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
661  first_op.Tvmult(*i, u);
662  second_op.Tvmult_add(v, *i);
663  };
664 
665  return return_op;
666  }
667 }
668 
669 
678 template <typename Range, typename Domain, typename Payload>
681 {
682  LinearOperator<Domain, Range, Payload> return_op{op.transpose_payload()};
683 
685  return_op.reinit_domain_vector = op.reinit_range_vector;
686 
687  return_op.vmult = op.Tvmult;
688  return_op.vmult_add = op.Tvmult_add;
689  return_op.Tvmult = op.vmult;
690  return_op.Tvmult_add = op.vmult_add;
691 
692  return return_op;
693 }
694 
695 
715 template <typename Payload,
716  typename Solver,
717  typename Preconditioner,
718  typename Range = typename Solver::vector_type,
719  typename Domain = Range>
722  Solver &solver,
723  const Preconditioner &preconditioner)
724 {
726  op.inverse_payload(solver, preconditioner)};
727 
729  return_op.reinit_domain_vector = op.reinit_range_vector;
730 
731  return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
732  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
733  solver.solve(op, v, u, preconditioner);
734  };
735 
736  return_op.vmult_add = [op, &solver, &preconditioner](Range &v,
737  const Domain &u) {
738  GrowingVectorMemory<Range> vector_memory;
739 
740  typename VectorMemory<Range>::Pointer v2(vector_memory);
741  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
742  solver.solve(op, *v2, u, preconditioner);
743  v += *v2;
744  };
745 
746  return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
747  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
748  solver.solve(transpose_operator(op), v, u, preconditioner);
749  };
750 
751  return_op.Tvmult_add = [op, &solver, &preconditioner](Range &v,
752  const Domain &u) {
753  GrowingVectorMemory<Range> vector_memory;
754 
755  typename VectorMemory<Range>::Pointer v2(vector_memory);
756  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
757  solver.solve(transpose_operator(op), *v2, u, preconditioner);
758  v += *v2;
759  };
760 
761  return return_op;
762 }
763 
764 
773 template <typename Payload,
774  typename Solver,
775  typename Range = typename Solver::vector_type,
776  typename Domain = Range>
779  Solver &solver,
780  const LinearOperator<Range, Domain, Payload> &preconditioner)
781 {
783  op.inverse_payload(solver, preconditioner)};
784 
786  return_op.reinit_domain_vector = op.reinit_range_vector;
787 
788  return_op.vmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
789  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
790  solver.solve(op, v, u, preconditioner);
791  };
792 
793  return_op.vmult_add = [op, &solver, preconditioner](Range &v,
794  const Domain &u) {
795  GrowingVectorMemory<Range> vector_memory;
796 
797  typename VectorMemory<Range>::Pointer v2(vector_memory);
798  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
799  solver.solve(op, *v2, u, preconditioner);
800  v += *v2;
801  };
802 
803  return_op.Tvmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
804  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
805  solver.solve(transpose_operator(op), v, u, preconditioner);
806  };
807 
808  return_op.Tvmult_add = [op, &solver, preconditioner](Range &v,
809  const Domain &u) {
810  GrowingVectorMemory<Range> vector_memory;
811 
812  typename VectorMemory<Range>::Pointer v2(vector_memory);
813  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
814  solver.solve(transpose_operator(op), *v2, u, preconditioner);
815  v += *v2;
816  };
817 
818  return return_op;
819 }
820 
821 
831 template <typename Payload,
832  typename Solver,
833  typename Range = typename Solver::vector_type,
834  typename Domain = Range>
837  Solver &solver)
838 {
839  return inverse_operator(op, solver, identity_operator(op));
840 }
841 
842 
851 template <typename Payload,
852  typename Solver,
853  typename Range = typename Solver::vector_type,
854  typename Domain = Range>
857  Solver &solver,
858  const PreconditionIdentity &)
859 {
860  return inverse_operator(op, solver);
861 }
862 
882 template <
883  typename Range,
886 identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
887 {
889 
890  return_op.reinit_range_vector = reinit_vector;
891  return_op.reinit_domain_vector = reinit_vector;
892 
893  return_op.vmult = [](Range &v, const Range &u) { v = u; };
894 
895  return_op.vmult_add = [](Range &v, const Range &u) { v += u; };
896 
897  return_op.Tvmult = [](Range &v, const Range &u) { v = u; };
898 
899  return_op.Tvmult_add = [](Range &v, const Range &u) { v += u; };
900 
901  return return_op;
902 }
903 
904 
917 template <typename Range, typename Domain, typename Payload>
920 {
921  auto return_op = identity_operator<Range, Payload>(op.reinit_range_vector);
922  static_cast<Payload &>(return_op) = op.identity_payload();
923 
924  return return_op;
925 }
926 
927 
937 template <typename Range, typename Domain, typename Payload>
940 {
941  LinearOperator<Range, Domain, Payload> return_op{op.null_payload()};
942 
943  return_op.is_null_operator = true;
944 
945  return_op.reinit_range_vector = op.reinit_range_vector;
946  return_op.reinit_domain_vector = op.reinit_domain_vector;
947 
948  return_op.vmult = [](Range &v, const Domain &) { v = 0.; };
949 
950  return_op.vmult_add = [](Range &, const Domain &) {};
951 
952  return_op.Tvmult = [](Domain &v, const Range &) { v = 0.; };
953 
954  return_op.Tvmult_add = [](Domain &, const Range &) {};
955 
956  return return_op;
957 }
958 
959 
972 template <
973  typename Range,
976 mean_value_filter(const std::function<void(Range &, bool)> &reinit_vector)
977 {
979 
980  return_op.reinit_range_vector = reinit_vector;
981  return_op.reinit_domain_vector = reinit_vector;
982 
983  return_op.vmult = [](Range &v, const Range &u) {
984  const auto mean = u.mean_value();
985 
986  v = u;
987  v.add(-mean);
988  };
989 
990  return_op.vmult_add = [](Range &v, const Range &u) {
991  const auto mean = u.mean_value();
992 
993  v += u;
994  v.add(-mean);
995  };
996 
997  return_op.Tvmult = return_op.vmult_add;
998  return_op.Tvmult_add = return_op.vmult_add;
999 
1000  return return_op;
1001 }
1002 
1003 
1016 template <typename Range, typename Domain, typename Payload>
1019 {
1020  auto return_op = mean_value_filter<Range, Payload>(op.reinit_range_vector);
1021  static_cast<Payload &>(return_op) = op.identity_payload();
1022 
1023  return return_op;
1024 }
1025 
1026 
1027 namespace internal
1028 {
1029  namespace LinearOperatorImplementation
1030  {
1042  template <typename Vector>
1044  {
1045  public:
1057  template <typename Matrix>
1058  static void
1060  Vector &v,
1061  bool omit_zeroing_entries)
1062  {
1063  v.reinit(matrix.m(), omit_zeroing_entries);
1064  }
1065 
1077  template <typename Matrix>
1078  static void
1080  Vector &v,
1081  bool omit_zeroing_entries)
1082  {
1083  v.reinit(matrix.n(), omit_zeroing_entries);
1084  }
1085  };
1086 
1087 
1100  {
1101  public:
1109  template <typename... Args>
1110  EmptyPayload(const Args &...)
1111  {}
1112 
1113 
1117  EmptyPayload
1119  {
1120  return *this;
1121  }
1122 
1123 
1127  EmptyPayload
1129  {
1130  return *this;
1131  }
1132 
1133 
1137  EmptyPayload
1139  {
1140  return *this;
1141  }
1142 
1143 
1147  template <typename Solver, typename Preconditioner>
1148  EmptyPayload
1149  inverse_payload(Solver &, const Preconditioner &) const
1150  {
1151  return *this;
1152  }
1153  };
1154 
1159  inline EmptyPayload
1161  {
1162  return {};
1163  }
1164 
1169  inline EmptyPayload
1171  {
1172  return {};
1173  }
1174 
1175 
1176 
1177  // A trait class that determines whether type T provides public
1178  // (templated or non-templated) vmult_add member functions
1179  template <typename Range, typename Domain, typename T>
1181  {
1182  template <typename C>
1183  static std::false_type
1184  test(...);
1185 
1186  template <typename C>
1187  static auto
1188  test(Range *r, Domain *d)
1189  -> decltype(std::declval<C>().vmult_add(*r, *d),
1190  std::declval<C>().Tvmult_add(*d, *r),
1191  std::true_type());
1192 
1193  public:
1194  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
1195  // otherwise it is std::false_type
1196 
1197  using type = decltype(test<T>(nullptr, nullptr));
1198  };
1199 
1200 
1201  // A helper function to apply a given vmult, or Tvmult to a vector with
1202  // intermediate storage
1203  template <typename Function, typename Range, typename Domain>
1204  void
1206  Range &v,
1207  const Domain &u,
1208  bool add)
1209  {
1210  GrowingVectorMemory<Range> vector_memory;
1211 
1212  typename VectorMemory<Range>::Pointer i(vector_memory);
1213  i->reinit(v, /*bool omit_zeroing_entries =*/true);
1214 
1215  function(*i, u);
1216 
1217  if (add)
1218  v += *i;
1219  else
1220  v = *i;
1221  }
1222 
1223 
1224  // A helper class to add a reduced matrix interface to a LinearOperator
1225  // (typically provided by Preconditioner classes)
1226  template <typename Range, typename Domain, typename Payload>
1228  {
1229  public:
1230  template <typename Matrix>
1231  void
1233  const Matrix &matrix)
1234  {
1235  op.vmult = [&matrix](Range &v, const Domain &u) {
1236  if (PointerComparison::equal(&v, &u))
1237  {
1238  // If v and u are the same memory location use intermediate
1239  // storage
1241  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1242  v,
1243  u,
1244  /*bool add =*/false);
1245  }
1246  else
1247  {
1248  matrix.vmult(v, u);
1249  }
1250  };
1251 
1252  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1253  // use intermediate storage to implement vmult_add with vmult
1255  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1256  v,
1257  u,
1258  /*bool add =*/true);
1259  };
1260 
1261  op.Tvmult = [&matrix](Domain &v, const Range &u) {
1262  if (PointerComparison::equal(&v, &u))
1263  {
1264  // If v and u are the same memory location use intermediate
1265  // storage
1267  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1268  v,
1269  u,
1270  /*bool add =*/false);
1271  }
1272  else
1273  {
1274  matrix.Tvmult(v, u);
1275  }
1276  };
1277 
1278  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1279  // use intermediate storage to implement Tvmult_add with Tvmult
1281  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1282  v,
1283  u,
1284  /*bool add =*/true);
1285  };
1286  }
1287  };
1288 
1289 
1290  // A helper class to add the full matrix interface to a LinearOperator
1291  template <typename Range, typename Domain, typename Payload>
1293  {
1294  public:
1295  template <typename Matrix>
1296  void
1298  const Matrix &matrix)
1299  {
1300  // As above ...
1301 
1303  op, matrix);
1304 
1305  // ... but add native vmult_add and Tvmult_add variants:
1306 
1307  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1308  if (PointerComparison::equal(&v, &u))
1309  {
1311  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1312  v,
1313  u,
1314  /*bool add =*/true);
1315  }
1316  else
1317  {
1318  matrix.vmult_add(v, u);
1319  }
1320  };
1321 
1322  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1323  if (PointerComparison::equal(&v, &u))
1324  {
1326  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1327  v,
1328  u,
1329  /*bool add =*/true);
1330  }
1331  else
1332  {
1333  matrix.Tvmult_add(v, u);
1334  }
1335  };
1336  }
1337  };
1338  } // namespace LinearOperatorImplementation
1339 } // namespace internal
1340 
1341 
1398 template <typename Range, typename Domain, typename Payload, typename Matrix>
1400 linear_operator(const Matrix &matrix)
1401 {
1402  // implement with the more generic variant below...
1403  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1404  matrix);
1405 }
1406 
1407 
1422 template <typename Range,
1423  typename Domain,
1424  typename Payload,
1425  typename OperatorExemplar,
1426  typename Matrix>
1428 linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
1429 {
1431  // Initialize the payload based on the input exemplar matrix
1433  Payload(operator_exemplar, matrix)};
1434 
1435  // Always store a reference to matrix and operator_exemplar in the lambda
1436  // functions. This ensures that a modification of the matrix after the
1437  // creation of a LinearOperator wrapper is respected - further a matrix
1438  // or an operator_exemplar cannot usually be copied...
1439 
1440  return_op.reinit_range_vector =
1441  [&operator_exemplar](Range &v, bool omit_zeroing_entries) {
1443  Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1444  };
1445 
1446  return_op.reinit_domain_vector = [&operator_exemplar](
1447  Domain &v, bool omit_zeroing_entries) {
1449  Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1450  };
1451 
1452  std::conditional_t<
1456  .
1457  operator()(return_op, matrix);
1458 
1459  return return_op;
1460 }
1461 
1462 
1463 
1480 template <typename Range, typename Domain, typename Payload, typename Matrix>
1483  const Matrix &matrix)
1484 {
1486  // Initialize the payload based on the LinearOperator exemplar
1487  auto return_op = operator_exemplar;
1488 
1489  std::conditional_t<
1493  .
1494  operator()(return_op, matrix);
1495 
1496  return return_op;
1497 }
1498 
1499 
1502 #ifndef DOXYGEN
1503 
1504 //
1505 // Ensure that we never capture a reference to a temporary by accident.
1506 // to avoid "stack use after free".
1507 //
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>>>
1517 linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1518 
1519 template <
1520  typename Range = Vector<double>,
1521  typename Domain = Range,
1523  typename OperatorExemplar,
1524  typename Matrix,
1525  typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
1526  typename = std::enable_if_t<
1527  !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
1529 linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1530 
1531 template <
1532  typename Range = Vector<double>,
1533  typename Domain = Range,
1535  typename OperatorExemplar,
1536  typename Matrix,
1537  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>,
1538  typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
1539  typename = std::enable_if_t<
1540  !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
1542 linear_operator(OperatorExemplar &&, Matrix &&) = delete;
1543 
1544 template <
1545  typename Range = Vector<double>,
1546  typename Domain = Range,
1548  typename Matrix,
1549  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1552  Matrix &&) = delete;
1553 
1554 template <
1555  typename Range = Vector<double>,
1556  typename Domain = Range,
1558  typename Matrix,
1559  typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1561 linear_operator(Matrix &&) = delete;
1562 
1563 template <
1564  typename Payload,
1565  typename Solver,
1566  typename Preconditioner,
1567  typename Range = typename Solver::vector_type,
1568  typename Domain = Range,
1569  typename = std::enable_if_t<!std::is_lvalue_reference_v<Preconditioner>>,
1570  typename =
1571  std::enable_if_t<!std::is_same_v<Preconditioner, PreconditionIdentity>>,
1572  typename = std::enable_if_t<
1573  !std::is_same_v<Preconditioner, LinearOperator<Range, Domain, Payload>>>>
1576  Solver &,
1577  Preconditioner &&) = delete;
1578 
1579 #endif // DOXYGEN
1580 
1582 
1583 #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:1631
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)