Reference documentation for deal.II version Git 689de043d4 2020-08-10 16:46:15 +0200
\(\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 - 2020 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 
167 template <typename Range, typename Domain, typename Payload>
168 class LinearOperator : public Payload
169 {
170 public:
179  LinearOperator(const Payload &payload = Payload())
180  : Payload(payload)
181  , is_null_operator(false)
182  {
183  vmult = [](Range &, const Domain &) {
184  Assert(false,
185  ExcMessage("Uninitialized LinearOperator<Range, "
186  "Domain>::vmult called"));
187  };
188 
189  vmult_add = [](Range &, const Domain &) {
190  Assert(false,
191  ExcMessage("Uninitialized LinearOperator<Range, "
192  "Domain>::vmult_add called"));
193  };
194 
195  Tvmult = [](Domain &, const Range &) {
196  Assert(false,
197  ExcMessage("Uninitialized LinearOperator<Range, "
198  "Domain>::Tvmult called"));
199  };
200 
201  Tvmult_add = [](Domain &, const Range &) {
202  Assert(false,
203  ExcMessage("Uninitialized LinearOperator<Range, "
204  "Domain>::Tvmult_add called"));
205  };
206 
207  reinit_range_vector = [](Range &, bool) {
208  Assert(false,
209  ExcMessage("Uninitialized LinearOperator<Range, "
210  "Domain>::reinit_range_vector method called"));
211  };
212 
213  reinit_domain_vector = [](Domain &, bool) {
214  Assert(false,
215  ExcMessage("Uninitialized LinearOperator<Range, "
216  "Domain>::reinit_domain_vector method called"));
217  };
218  }
219 
224 
230  template <typename Op,
231  typename = typename std::enable_if<
232  !std::is_base_of<LinearOperator<Range, Domain, Payload>,
233  Op>::value>::type>
234  LinearOperator(const Op &op)
235  {
236  *this = linear_operator<Range, Domain, Payload, Op>(op);
237  }
238 
243  operator=(const LinearOperator<Range, Domain, Payload> &) = default;
244 
249  template <typename Op,
250  typename = typename std::enable_if<
251  !std::is_base_of<LinearOperator<Range, Domain, Payload>,
252  Op>::value>::type>
254  operator=(const Op &op)
255  {
256  *this = linear_operator<Range, Domain, Payload, Op>(op);
257  return *this;
258  }
259 
264  std::function<void(Range &v, const Domain &u)> vmult;
265 
270  std::function<void(Range &v, const Domain &u)> vmult_add;
271 
276  std::function<void(Domain &v, const Range &u)> Tvmult;
277 
282  std::function<void(Domain &v, const Range &u)> Tvmult_add;
283 
291  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_range_vector;
292 
300  std::function<void(Domain &v, bool omit_zeroing_entries)>
302 
307 
314  {
315  *this = *this + second_op;
316  return *this;
317  }
318 
325  {
326  *this = *this - second_op;
327  return *this;
328  }
329 
336  {
337  *this = *this * second_op;
338  return *this;
339  }
340 
346  operator*=(typename Domain::value_type number)
347  {
348  *this = *this * number;
349  return *this;
350  }
351 
358 
360 };
361 
362 
367 
377 template <typename Range, typename Domain, typename Payload>
381 {
382  if (first_op.is_null_operator)
383  {
384  return second_op;
385  }
386  else if (second_op.is_null_operator)
387  {
388  return first_op;
389  }
390  else
391  {
393  static_cast<const Payload &>(first_op) +
394  static_cast<const Payload &>(second_op)};
395 
396  return_op.reinit_range_vector = first_op.reinit_range_vector;
397  return_op.reinit_domain_vector = first_op.reinit_domain_vector;
398 
399  // ensure to have valid computation objects by catching first_op and
400  // second_op by value
401 
402  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
403  first_op.vmult(v, u);
404  second_op.vmult_add(v, u);
405  };
406 
407  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
408  first_op.vmult_add(v, u);
409  second_op.vmult_add(v, u);
410  };
411 
412  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
413  second_op.Tvmult(v, u);
414  first_op.Tvmult_add(v, u);
415  };
416 
417  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
418  second_op.Tvmult_add(v, u);
419  first_op.Tvmult_add(v, u);
420  };
421 
422  return return_op;
423  }
424 }
425 
426 
436 template <typename Range, typename Domain, typename Payload>
440 {
441  if (first_op.is_null_operator)
442  {
443  return -1. * second_op;
444  }
445  else if (second_op.is_null_operator)
446  {
447  return first_op;
448  }
449  else
450  {
451  // implement with addition and scalar multiplication
452  return first_op + (-1. * second_op);
453  }
454 }
455 
456 
474 template <typename Range, typename Domain, typename Payload>
476 operator*(typename Range::value_type number,
478 {
479  static_assert(
480  std::is_convertible<typename Range::value_type,
481  typename Domain::value_type>::value,
482  "Range and Domain must have implicitly convertible 'value_type's");
483 
484  if (op.is_null_operator)
485  {
486  return op;
487  }
488  else if (number == 0.)
489  {
490  return null_operator(op);
491  }
492  else
493  {
495 
496  // ensure to have valid computation objects by catching number and op by
497  // value
498 
499  return_op.vmult = [number, op](Range &v, const Domain &u) {
500  op.vmult(v, u);
501  v *= number;
502  };
503 
504  return_op.vmult_add = [number, op](Range &v, const Domain &u) {
505  v /= number;
506  op.vmult_add(v, u);
507  v *= number;
508  };
509 
510  return_op.Tvmult = [number, op](Domain &v, const Range &u) {
511  op.Tvmult(v, u);
512  v *= number;
513  };
514 
515  return_op.Tvmult_add = [number, op](Domain &v, const Range &u) {
516  v /= number;
517  op.Tvmult_add(v, u);
518  v *= number;
519  };
520 
521  return return_op;
522  }
523 }
524 
525 
541 template <typename Range, typename Domain, typename Payload>
544  typename Domain::value_type number)
545 {
546  static_assert(
547  std::is_convertible<typename Range::value_type,
548  typename Domain::value_type>::value,
549  "Range and Domain must have implicitly convertible 'value_type's");
550 
551  return number * op;
552 }
553 
555 
556 
561 
571 template <typename Range,
572  typename Intermediate,
573  typename Domain,
574  typename Payload>
578 {
579  if (first_op.is_null_operator || second_op.is_null_operator)
580  {
582  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
583  return_op.reinit_range_vector = first_op.reinit_range_vector;
584  return null_operator(return_op);
585  }
586  else
587  {
589  static_cast<const Payload &>(first_op) *
590  static_cast<const Payload &>(second_op)};
591 
592  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
593  return_op.reinit_range_vector = first_op.reinit_range_vector;
594 
595  // ensure to have valid computation objects by catching first_op and
596  // second_op by value
597 
598  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
599  GrowingVectorMemory<Intermediate> vector_memory;
600 
601  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
602  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
603  second_op.vmult(*i, u);
604  first_op.vmult(v, *i);
605  };
606 
607  return_op.vmult_add = [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_add(v, *i);
614  };
615 
616  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
617  GrowingVectorMemory<Intermediate> vector_memory;
618 
619  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
620  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
621  first_op.Tvmult(*i, u);
622  second_op.Tvmult(v, *i);
623  };
624 
625  return_op.Tvmult_add = [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_add(v, *i);
632  };
633 
634  return return_op;
635  }
636 }
637 
638 
647 template <typename Range, typename Domain, typename Payload>
650 {
651  LinearOperator<Domain, Range, Payload> return_op{op.transpose_payload()};
652 
654  return_op.reinit_domain_vector = op.reinit_range_vector;
655 
656  return_op.vmult = op.Tvmult;
657  return_op.vmult_add = op.Tvmult_add;
658  return_op.Tvmult = op.vmult;
659  return_op.Tvmult_add = op.vmult_add;
660 
661  return return_op;
662 }
663 
664 
684 template <typename Payload,
685  typename Solver,
686  typename Preconditioner,
687  typename Range = typename Solver::vector_type,
688  typename Domain = Range>
691  Solver & solver,
692  const Preconditioner & preconditioner)
693 {
695  op.inverse_payload(solver, preconditioner)};
696 
698  return_op.reinit_domain_vector = op.reinit_range_vector;
699 
700  return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
701  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
702  solver.solve(op, v, u, preconditioner);
703  };
704 
705  return_op.vmult_add = [op, &solver, &preconditioner](Range & v,
706  const Domain &u) {
707  GrowingVectorMemory<Range> vector_memory;
708 
709  typename VectorMemory<Range>::Pointer v2(vector_memory);
710  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
711  solver.solve(op, *v2, u, preconditioner);
712  v += *v2;
713  };
714 
715  return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
716  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
717  solver.solve(transpose_operator(op), v, u, preconditioner);
718  };
719 
720  return_op.Tvmult_add = [op, &solver, &preconditioner](Range & v,
721  const Domain &u) {
722  GrowingVectorMemory<Range> vector_memory;
723 
724  typename VectorMemory<Range>::Pointer v2(vector_memory);
725  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
726  solver.solve(transpose_operator(op), *v2, u, preconditioner);
727  v += *v2;
728  };
729 
730  return return_op;
731 }
732 
733 
742 template <typename Payload,
743  typename Solver,
744  typename Range = typename Solver::vector_type,
745  typename Domain = Range>
748  Solver & solver,
749  const LinearOperator<Range, Domain, Payload> &preconditioner)
750 {
752  op.inverse_payload(solver, preconditioner)};
753 
755  return_op.reinit_domain_vector = op.reinit_range_vector;
756 
757  return_op.vmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
758  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
759  solver.solve(op, v, u, preconditioner);
760  };
761 
762  return_op.vmult_add = [op, &solver, preconditioner](Range & v,
763  const Domain &u) {
764  GrowingVectorMemory<Range> vector_memory;
765 
766  typename VectorMemory<Range>::Pointer v2(vector_memory);
767  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
768  solver.solve(op, *v2, u, preconditioner);
769  v += *v2;
770  };
771 
772  return_op.Tvmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
773  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
774  solver.solve(transpose_operator(op), v, u, preconditioner);
775  };
776 
777  return_op.Tvmult_add = [op, &solver, preconditioner](Range & v,
778  const Domain &u) {
779  GrowingVectorMemory<Range> vector_memory;
780 
781  typename VectorMemory<Range>::Pointer v2(vector_memory);
782  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
783  solver.solve(transpose_operator(op), *v2, u, preconditioner);
784  v += *v2;
785  };
786 
787  return return_op;
788 }
789 
790 
800 template <typename Payload,
801  typename Solver,
802  typename Range = typename Solver::vector_type,
803  typename Domain = Range>
806  Solver & solver)
807 {
808  return inverse_operator(op, solver, identity_operator(op));
809 }
810 
811 
820 template <typename Payload,
821  typename Solver,
822  typename Range = typename Solver::vector_type,
823  typename Domain = Range>
826  Solver & solver,
827  const PreconditionIdentity &)
828 {
829  return inverse_operator(op, solver);
830 }
831 
833 
834 
839 
851 template <
852  typename Range,
855 identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
856 {
857  LinearOperator<Range, Range, Payload> return_op{Payload()};
858 
859  return_op.reinit_range_vector = reinit_vector;
860  return_op.reinit_domain_vector = reinit_vector;
861 
862  return_op.vmult = [](Range &v, const Range &u) { v = u; };
863 
864  return_op.vmult_add = [](Range &v, const Range &u) { v += u; };
865 
866  return_op.Tvmult = [](Range &v, const Range &u) { v = u; };
867 
868  return_op.Tvmult_add = [](Range &v, const Range &u) { v += u; };
869 
870  return return_op;
871 }
872 
873 
886 template <typename Range, typename Domain, typename Payload>
889 {
890  auto return_op = identity_operator<Range, Payload>(op.reinit_range_vector);
891  static_cast<Payload &>(return_op) = op.identity_payload();
892 
893  return return_op;
894 }
895 
896 
906 template <typename Range, typename Domain, typename Payload>
909 {
910  LinearOperator<Range, Domain, Payload> return_op{op.null_payload()};
911 
912  return_op.is_null_operator = true;
913 
914  return_op.reinit_range_vector = op.reinit_range_vector;
915  return_op.reinit_domain_vector = op.reinit_domain_vector;
916 
917  return_op.vmult = [](Range &v, const Domain &) { v = 0.; };
918 
919  return_op.vmult_add = [](Range &, const Domain &) {};
920 
921  return_op.Tvmult = [](Domain &v, const Range &) { v = 0.; };
922 
923  return_op.Tvmult_add = [](Domain &, const Range &) {};
924 
925  return return_op;
926 }
927 
928 
941 template <
942  typename Range,
945 mean_value_filter(const std::function<void(Range &, bool)> &reinit_vector)
946 {
947  LinearOperator<Range, Range, Payload> return_op{Payload()};
948 
949  return_op.reinit_range_vector = reinit_vector;
950  return_op.reinit_domain_vector = reinit_vector;
951 
952  return_op.vmult = [](Range &v, const Range &u) {
953  const auto mean = u.mean_value();
954 
955  v = u;
956  v.add(-mean);
957  };
958 
959  return_op.vmult_add = [](Range &v, const Range &u) {
960  const auto mean = u.mean_value();
961 
962  v += u;
963  v.add(-mean);
964  };
965 
966  return_op.Tvmult = return_op.vmult_add;
967  return_op.Tvmult_add = return_op.vmult_add;
968 
969  return return_op;
970 }
971 
972 
985 template <typename Range, typename Domain, typename Payload>
988 {
989  auto return_op = mean_value_filter<Range, Payload>(op.reinit_range_vector);
990  static_cast<Payload &>(return_op) = op.identity_payload();
991 
992  return return_op;
993 }
994 
995 
996 namespace internal
997 {
998  namespace LinearOperatorImplementation
999  {
1011  template <typename Vector>
1012  class ReinitHelper
1013  {
1014  public:
1026  template <typename Matrix>
1027  static void
1029  Vector & v,
1030  bool omit_zeroing_entries)
1031  {
1032  v.reinit(matrix.m(), omit_zeroing_entries);
1033  }
1034 
1046  template <typename Matrix>
1047  static void
1049  Vector & v,
1050  bool omit_zeroing_entries)
1051  {
1052  v.reinit(matrix.n(), omit_zeroing_entries);
1053  }
1054  };
1055 
1056 
1069  {
1070  public:
1078  template <typename... Args>
1079  EmptyPayload(const Args &...)
1080  {}
1081 
1082 
1086  EmptyPayload
1088  {
1089  return *this;
1090  }
1091 
1092 
1096  EmptyPayload
1098  {
1099  return *this;
1100  }
1101 
1102 
1106  EmptyPayload
1108  {
1109  return *this;
1110  }
1111 
1112 
1116  template <typename Solver, typename Preconditioner>
1117  EmptyPayload
1118  inverse_payload(Solver &, const Preconditioner &) const
1119  {
1120  return *this;
1121  }
1122  };
1123 
1128  inline EmptyPayload
1130  {
1131  return {};
1132  }
1133 
1139  {
1140  return {};
1141  }
1142 
1143 
1144 
1145  // A trait class that determines whether type T provides public
1146  // (templated or non-templated) vmult_add member functions
1147  template <typename Range, typename Domain, typename T>
1149  {
1150  template <typename C>
1151  static std::false_type
1152  test(...);
1153 
1154  template <typename C>
1155  static auto
1156  test(Range *r, Domain *d)
1157  -> decltype(std::declval<C>().vmult_add(*r, *d),
1158  std::declval<C>().Tvmult_add(*d, *r),
1159  std::true_type());
1160 
1161  public:
1162  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
1163  // otherwise it is std::false_type
1164 
1165  using type = decltype(test<T>(nullptr, nullptr));
1166  };
1167 
1168 
1169  // A helper function to apply a given vmult, or Tvmult to a vector with
1170  // intermediate storage
1171  template <typename Function, typename Range, typename Domain>
1172  void
1174  Range & v,
1175  const Domain &u,
1176  bool add)
1177  {
1178  GrowingVectorMemory<Range> vector_memory;
1179 
1180  typename VectorMemory<Range>::Pointer i(vector_memory);
1181  i->reinit(v, /*bool omit_zeroing_entries =*/true);
1182 
1183  function(*i, u);
1184 
1185  if (add)
1186  v += *i;
1187  else
1188  v = *i;
1189  }
1190 
1191 
1192  // A helper class to add a reduced matrix interface to a LinearOperator
1193  // (typically provided by Preconditioner classes)
1194  template <typename Range, typename Domain, typename Payload>
1196  {
1197  public:
1198  template <typename Matrix>
1199  void
1201  const Matrix & matrix)
1202  {
1203  op.vmult = [&matrix](Range &v, const Domain &u) {
1204  if (PointerComparison::equal(&v, &u))
1205  {
1206  // If v and u are the same memory location use intermediate
1207  // storage
1209  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1210  v,
1211  u,
1212  /*bool add =*/false);
1213  }
1214  else
1215  {
1216  matrix.vmult(v, u);
1217  }
1218  };
1219 
1220  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1221  // use intermediate storage to implement vmult_add with vmult
1223  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1224  v,
1225  u,
1226  /*bool add =*/true);
1227  };
1228 
1229  op.Tvmult = [&matrix](Domain &v, const Range &u) {
1230  if (PointerComparison::equal(&v, &u))
1231  {
1232  // If v and u are the same memory location use intermediate
1233  // storage
1235  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1236  v,
1237  u,
1238  /*bool add =*/false);
1239  }
1240  else
1241  {
1242  matrix.Tvmult(v, u);
1243  }
1244  };
1245 
1246  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1247  // use intermediate storage to implement Tvmult_add with Tvmult
1249  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1250  v,
1251  u,
1252  /*bool add =*/true);
1253  };
1254  }
1255  };
1256 
1257 
1258  // A helper class to add the full matrix interface to a LinearOperator
1259  template <typename Range, typename Domain, typename Payload>
1261  {
1262  public:
1263  template <typename Matrix>
1264  void
1266  const Matrix & matrix)
1267  {
1268  // As above ...
1269 
1271  op, matrix);
1272 
1273  // ... but add native vmult_add and Tvmult_add variants:
1274 
1275  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1276  if (PointerComparison::equal(&v, &u))
1277  {
1279  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1280  v,
1281  u,
1282  /*bool add =*/true);
1283  }
1284  else
1285  {
1286  matrix.vmult_add(v, u);
1287  }
1288  };
1289 
1290  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1291  if (PointerComparison::equal(&v, &u))
1292  {
1294  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1295  v,
1296  u,
1297  /*bool add =*/true);
1298  }
1299  else
1300  {
1301  matrix.Tvmult_add(v, u);
1302  }
1303  };
1304  }
1305  };
1306  } // namespace LinearOperatorImplementation
1307 } // namespace internal
1308 
1309 
1366 template <typename Range, typename Domain, typename Payload, typename Matrix>
1368 linear_operator(const Matrix &matrix)
1369 {
1370  // implement with the more generic variant below...
1371  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1372  matrix);
1373 }
1374 
1375 
1390 template <typename Range,
1391  typename Domain,
1392  typename Payload,
1393  typename OperatorExemplar,
1394  typename Matrix>
1396 linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
1397 {
1399  // Initialize the payload based on the input exemplar matrix
1401  Payload(operator_exemplar, matrix)};
1402 
1403  // Always store a reference to matrix and operator_exemplar in the lambda
1404  // functions. This ensures that a modification of the matrix after the
1405  // creation of a LinearOperator wrapper is respected - further a matrix
1406  // or an operator_exemplar cannot usually be copied...
1407 
1408  return_op.reinit_range_vector =
1409  [&operator_exemplar](Range &v, bool omit_zeroing_entries) {
1411  Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1412  };
1413 
1414  return_op.reinit_domain_vector = [&operator_exemplar](
1415  Domain &v, bool omit_zeroing_entries) {
1417  Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1418  };
1419 
1420  typename std::conditional<
1424  .
1425  operator()(return_op, matrix);
1426 
1427  return return_op;
1428 }
1429 
1430 
1431 
1448 template <typename Range, typename Domain, typename Payload, typename Matrix>
1451  const Matrix & matrix)
1452 {
1454  // Initialize the payload based on the LinearOperator exemplar
1455  auto return_op = operator_exemplar;
1456 
1457  typename std::conditional<
1461  .
1462  operator()(return_op, matrix);
1463 
1464  return return_op;
1465 }
1466 
1467 
1469 
1470 #ifndef DOXYGEN
1471 
1472 //
1473 // Ensure that we never capture a reference to a temporary by accident.
1474 // to avoid "stack use after free".
1475 //
1476 
1477 template <
1478  typename Range = Vector<double>,
1479  typename Domain = Range,
1481  typename OperatorExemplar,
1482  typename Matrix,
1483  typename =
1484  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1486 linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1487 
1488 template <
1489  typename Range = Vector<double>,
1490  typename Domain = Range,
1492  typename OperatorExemplar,
1493  typename Matrix,
1494  typename = typename std::enable_if<
1495  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1496  typename = typename std::enable_if<
1497  !std::is_same<OperatorExemplar,
1500 linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1501 
1502 template <
1503  typename Range = Vector<double>,
1504  typename Domain = Range,
1506  typename OperatorExemplar,
1507  typename Matrix,
1508  typename =
1509  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type,
1510  typename = typename std::enable_if<
1511  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1512  typename = typename std::enable_if<
1513  !std::is_same<OperatorExemplar,
1516 linear_operator(OperatorExemplar &&, Matrix &&) = delete;
1517 
1518 template <
1519  typename Range = Vector<double>,
1520  typename Domain = Range,
1522  typename Matrix,
1523  typename =
1524  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1527  Matrix &&) = delete;
1528 
1529 template <
1530  typename Range = Vector<double>,
1531  typename Domain = Range,
1533  typename Matrix,
1534  typename =
1535  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1537 linear_operator(Matrix &&) = delete;
1538 
1539 template <typename Payload,
1540  typename Solver,
1541  typename Preconditioner,
1542  typename Range = typename Solver::vector_type,
1543  typename Domain = Range,
1544  typename = typename std::enable_if<
1545  !std::is_lvalue_reference<Preconditioner>::value>::type,
1546  typename = typename std::enable_if<
1547  !std::is_same<Preconditioner, PreconditionIdentity>::value>::type,
1548  typename = typename std::enable_if<
1549  !std::is_same<Preconditioner,
1553  Solver &,
1554  Preconditioner &&) = delete;
1555 
1556 #endif // DOXYGEN
1557 
1559 
1560 #endif
Contents is actually a matrix.
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< Domain, Domain, Payload > &second_op)
void operator()(LinearOperator< Range, Domain, Payload > &op, const Matrix &matrix)
EmptyPayload inverse_payload(Solver &, const Preconditioner &) const
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)
void apply_with_intermediate_storage(Function function, Range &v, const Domain &u, bool add)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1411
std::function< void(Range &v, const Domain &u)> vmult
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
EmptyPayload operator*(const EmptyPayload &, const EmptyPayload &)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
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 > &second_op)
LinearOperator(const Payload &payload=Payload())
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
LinearOperator< Range, Domain, Payload > identity_operator(const LinearOperator< Range, Domain, Payload > &)
void operator()(LinearOperator< Range, Domain, Payload > &op, const Matrix &matrix)
static bool equal(const T *p1, const T *p2)
LinearOperator(const Op &op)
EmptyPayload operator+(const EmptyPayload &, const EmptyPayload &)