Reference documentation for deal.II version Git 44bda21415 2021-03-04 10:32:54 +0100
\(\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 
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 = typename std::enable_if<
243  !std::is_base_of<LinearOperator<Range, Domain, Payload>,
244  Op>::value>::type>
245  LinearOperator(const Op &op)
246  {
247  *this = linear_operator<Range, Domain, Payload, Op>(op);
248  }
249 
254  operator=(const LinearOperator<Range, Domain, Payload> &) = default;
255 
260  template <typename Op,
261  typename = typename std::enable_if<
262  !std::is_base_of<LinearOperator<Range, Domain, Payload>,
263  Op>::value>::type>
265  operator=(const Op &op)
266  {
267  *this = linear_operator<Range, Domain, Payload, Op>(op);
268  return *this;
269  }
270 
275  std::function<void(Range &v, const Domain &u)> vmult;
276 
281  std::function<void(Range &v, const Domain &u)> vmult_add;
282 
287  std::function<void(Domain &v, const Range &u)> Tvmult;
288 
293  std::function<void(Domain &v, const Range &u)> Tvmult_add;
294 
302  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_range_vector;
303 
311  std::function<void(Domain &v, bool omit_zeroing_entries)>
313 
318 
325  {
326  *this = *this + second_op;
327  return *this;
328  }
329 
336  {
337  *this = *this - second_op;
338  return *this;
339  }
340 
347  {
348  *this = *this * second_op;
349  return *this;
350  }
351 
357  operator*=(typename Domain::value_type number)
358  {
359  *this = *this * number;
360  return *this;
361  }
362 
369 
371 };
372 
373 
378 
388 template <typename Range, typename Domain, typename Payload>
392 {
393  if (first_op.is_null_operator)
394  {
395  return second_op;
396  }
397  else if (second_op.is_null_operator)
398  {
399  return first_op;
400  }
401  else
402  {
404  static_cast<const Payload &>(first_op) +
405  static_cast<const Payload &>(second_op)};
406 
407  return_op.reinit_range_vector = first_op.reinit_range_vector;
408  return_op.reinit_domain_vector = first_op.reinit_domain_vector;
409 
410  // ensure to have valid computation objects by catching first_op and
411  // second_op by value
412 
413  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
414  first_op.vmult(v, u);
415  second_op.vmult_add(v, u);
416  };
417 
418  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
419  first_op.vmult_add(v, u);
420  second_op.vmult_add(v, u);
421  };
422 
423  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
424  second_op.Tvmult(v, u);
425  first_op.Tvmult_add(v, u);
426  };
427 
428  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
429  second_op.Tvmult_add(v, u);
430  first_op.Tvmult_add(v, u);
431  };
432 
433  return return_op;
434  }
435 }
436 
437 
447 template <typename Range, typename Domain, typename Payload>
451 {
452  if (first_op.is_null_operator)
453  {
454  return -1. * second_op;
455  }
456  else if (second_op.is_null_operator)
457  {
458  return first_op;
459  }
460  else
461  {
462  // implement with addition and scalar multiplication
463  return first_op + (-1. * second_op);
464  }
465 }
466 
467 
485 template <typename Range, typename Domain, typename Payload>
487 operator*(typename Range::value_type number,
489 {
490  static_assert(
491  std::is_convertible<typename Range::value_type,
492  typename Domain::value_type>::value,
493  "Range and Domain must have implicitly convertible 'value_type's");
494 
495  if (op.is_null_operator)
496  {
497  return op;
498  }
499  else if (number == 0.)
500  {
501  return null_operator(op);
502  }
503  else
504  {
506 
507  // ensure to have valid computation objects by catching number and op by
508  // value
509 
510  return_op.vmult = [number, op](Range &v, const Domain &u) {
511  op.vmult(v, u);
512  v *= number;
513  };
514 
515  return_op.vmult_add = [number, op](Range &v, const Domain &u) {
516  v /= number;
517  op.vmult_add(v, u);
518  v *= number;
519  };
520 
521  return_op.Tvmult = [number, op](Domain &v, const Range &u) {
522  op.Tvmult(v, u);
523  v *= number;
524  };
525 
526  return_op.Tvmult_add = [number, op](Domain &v, const Range &u) {
527  v /= number;
528  op.Tvmult_add(v, u);
529  v *= number;
530  };
531 
532  return return_op;
533  }
534 }
535 
536 
552 template <typename Range, typename Domain, typename Payload>
555  typename Domain::value_type number)
556 {
557  static_assert(
558  std::is_convertible<typename Range::value_type,
559  typename Domain::value_type>::value,
560  "Range and Domain must have implicitly convertible 'value_type's");
561 
562  return number * op;
563 }
564 
566 
567 
572 
582 template <typename Range,
583  typename Intermediate,
584  typename Domain,
585  typename Payload>
589 {
590  if (first_op.is_null_operator || second_op.is_null_operator)
591  {
593  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
594  return_op.reinit_range_vector = first_op.reinit_range_vector;
595  return null_operator(return_op);
596  }
597  else
598  {
600  static_cast<const Payload &>(first_op) *
601  static_cast<const Payload &>(second_op)};
602 
603  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
604  return_op.reinit_range_vector = first_op.reinit_range_vector;
605 
606  // ensure to have valid computation objects by catching first_op and
607  // second_op by value
608 
609  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
610  GrowingVectorMemory<Intermediate> vector_memory;
611 
612  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
613  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
614  second_op.vmult(*i, u);
615  first_op.vmult(v, *i);
616  };
617 
618  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
619  GrowingVectorMemory<Intermediate> vector_memory;
620 
621  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
622  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
623  second_op.vmult(*i, u);
624  first_op.vmult_add(v, *i);
625  };
626 
627  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
628  GrowingVectorMemory<Intermediate> vector_memory;
629 
630  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
631  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
632  first_op.Tvmult(*i, u);
633  second_op.Tvmult(v, *i);
634  };
635 
636  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
637  GrowingVectorMemory<Intermediate> vector_memory;
638 
639  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
640  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
641  first_op.Tvmult(*i, u);
642  second_op.Tvmult_add(v, *i);
643  };
644 
645  return return_op;
646  }
647 }
648 
649 
658 template <typename Range, typename Domain, typename Payload>
661 {
662  LinearOperator<Domain, Range, Payload> return_op{op.transpose_payload()};
663 
665  return_op.reinit_domain_vector = op.reinit_range_vector;
666 
667  return_op.vmult = op.Tvmult;
668  return_op.vmult_add = op.Tvmult_add;
669  return_op.Tvmult = op.vmult;
670  return_op.Tvmult_add = op.vmult_add;
671 
672  return return_op;
673 }
674 
675 
695 template <typename Payload,
696  typename Solver,
697  typename Preconditioner,
698  typename Range = typename Solver::vector_type,
699  typename Domain = Range>
702  Solver & solver,
703  const Preconditioner & preconditioner)
704 {
706  op.inverse_payload(solver, preconditioner)};
707 
709  return_op.reinit_domain_vector = op.reinit_range_vector;
710 
711  return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
712  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
713  solver.solve(op, v, u, preconditioner);
714  };
715 
716  return_op.vmult_add = [op, &solver, &preconditioner](Range & v,
717  const Domain &u) {
718  GrowingVectorMemory<Range> vector_memory;
719 
720  typename VectorMemory<Range>::Pointer v2(vector_memory);
721  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
722  solver.solve(op, *v2, u, preconditioner);
723  v += *v2;
724  };
725 
726  return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
727  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
728  solver.solve(transpose_operator(op), v, u, preconditioner);
729  };
730 
731  return_op.Tvmult_add = [op, &solver, &preconditioner](Range & v,
732  const Domain &u) {
733  GrowingVectorMemory<Range> vector_memory;
734 
735  typename VectorMemory<Range>::Pointer v2(vector_memory);
736  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
737  solver.solve(transpose_operator(op), *v2, u, preconditioner);
738  v += *v2;
739  };
740 
741  return return_op;
742 }
743 
744 
753 template <typename Payload,
754  typename Solver,
755  typename Range = typename Solver::vector_type,
756  typename Domain = Range>
759  Solver & solver,
760  const LinearOperator<Range, Domain, Payload> &preconditioner)
761 {
763  op.inverse_payload(solver, preconditioner)};
764 
766  return_op.reinit_domain_vector = op.reinit_range_vector;
767 
768  return_op.vmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
769  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
770  solver.solve(op, v, u, preconditioner);
771  };
772 
773  return_op.vmult_add = [op, &solver, preconditioner](Range & v,
774  const Domain &u) {
775  GrowingVectorMemory<Range> vector_memory;
776 
777  typename VectorMemory<Range>::Pointer v2(vector_memory);
778  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
779  solver.solve(op, *v2, u, preconditioner);
780  v += *v2;
781  };
782 
783  return_op.Tvmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
784  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
785  solver.solve(transpose_operator(op), v, u, preconditioner);
786  };
787 
788  return_op.Tvmult_add = [op, &solver, preconditioner](Range & v,
789  const Domain &u) {
790  GrowingVectorMemory<Range> vector_memory;
791 
792  typename VectorMemory<Range>::Pointer v2(vector_memory);
793  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
794  solver.solve(transpose_operator(op), *v2, u, preconditioner);
795  v += *v2;
796  };
797 
798  return return_op;
799 }
800 
801 
811 template <typename Payload,
812  typename Solver,
813  typename Range = typename Solver::vector_type,
814  typename Domain = Range>
817  Solver & solver)
818 {
819  return inverse_operator(op, solver, identity_operator(op));
820 }
821 
822 
831 template <typename Payload,
832  typename Solver,
833  typename Range = typename Solver::vector_type,
834  typename Domain = Range>
837  Solver & solver,
838  const PreconditionIdentity &)
839 {
840  return inverse_operator(op, solver);
841 }
842 
844 
845 
850 
862 template <
863  typename Range,
866 identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
867 {
868  LinearOperator<Range, Range, Payload> return_op{Payload()};
869 
870  return_op.reinit_range_vector = reinit_vector;
871  return_op.reinit_domain_vector = reinit_vector;
872 
873  return_op.vmult = [](Range &v, const Range &u) { v = u; };
874 
875  return_op.vmult_add = [](Range &v, const Range &u) { v += u; };
876 
877  return_op.Tvmult = [](Range &v, const Range &u) { v = u; };
878 
879  return_op.Tvmult_add = [](Range &v, const Range &u) { v += u; };
880 
881  return return_op;
882 }
883 
884 
897 template <typename Range, typename Domain, typename Payload>
900 {
901  auto return_op = identity_operator<Range, Payload>(op.reinit_range_vector);
902  static_cast<Payload &>(return_op) = op.identity_payload();
903 
904  return return_op;
905 }
906 
907 
917 template <typename Range, typename Domain, typename Payload>
920 {
921  LinearOperator<Range, Domain, Payload> return_op{op.null_payload()};
922 
923  return_op.is_null_operator = true;
924 
925  return_op.reinit_range_vector = op.reinit_range_vector;
926  return_op.reinit_domain_vector = op.reinit_domain_vector;
927 
928  return_op.vmult = [](Range &v, const Domain &) { v = 0.; };
929 
930  return_op.vmult_add = [](Range &, const Domain &) {};
931 
932  return_op.Tvmult = [](Domain &v, const Range &) { v = 0.; };
933 
934  return_op.Tvmult_add = [](Domain &, const Range &) {};
935 
936  return return_op;
937 }
938 
939 
952 template <
953  typename Range,
956 mean_value_filter(const std::function<void(Range &, bool)> &reinit_vector)
957 {
958  LinearOperator<Range, Range, Payload> return_op{Payload()};
959 
960  return_op.reinit_range_vector = reinit_vector;
961  return_op.reinit_domain_vector = reinit_vector;
962 
963  return_op.vmult = [](Range &v, const Range &u) {
964  const auto mean = u.mean_value();
965 
966  v = u;
967  v.add(-mean);
968  };
969 
970  return_op.vmult_add = [](Range &v, const Range &u) {
971  const auto mean = u.mean_value();
972 
973  v += u;
974  v.add(-mean);
975  };
976 
977  return_op.Tvmult = return_op.vmult_add;
978  return_op.Tvmult_add = return_op.vmult_add;
979 
980  return return_op;
981 }
982 
983 
996 template <typename Range, typename Domain, typename Payload>
999 {
1000  auto return_op = mean_value_filter<Range, Payload>(op.reinit_range_vector);
1001  static_cast<Payload &>(return_op) = op.identity_payload();
1002 
1003  return return_op;
1004 }
1005 
1006 
1007 namespace internal
1008 {
1009  namespace LinearOperatorImplementation
1010  {
1022  template <typename Vector>
1023  class ReinitHelper
1024  {
1025  public:
1037  template <typename Matrix>
1038  static void
1040  Vector & v,
1041  bool omit_zeroing_entries)
1042  {
1043  v.reinit(matrix.m(), omit_zeroing_entries);
1044  }
1045 
1057  template <typename Matrix>
1058  static void
1060  Vector & v,
1061  bool omit_zeroing_entries)
1062  {
1063  v.reinit(matrix.n(), omit_zeroing_entries);
1064  }
1065  };
1066 
1067 
1080  {
1081  public:
1089  template <typename... Args>
1090  EmptyPayload(const Args &...)
1091  {}
1092 
1093 
1097  EmptyPayload
1099  {
1100  return *this;
1101  }
1102 
1103 
1107  EmptyPayload
1109  {
1110  return *this;
1111  }
1112 
1113 
1117  EmptyPayload
1119  {
1120  return *this;
1121  }
1122 
1123 
1127  template <typename Solver, typename Preconditioner>
1128  EmptyPayload
1129  inverse_payload(Solver &, const Preconditioner &) const
1130  {
1131  return *this;
1132  }
1133  };
1134 
1139  inline EmptyPayload
1141  {
1142  return {};
1143  }
1144 
1150  {
1151  return {};
1152  }
1153 
1154 
1155 
1156  // A trait class that determines whether type T provides public
1157  // (templated or non-templated) vmult_add member functions
1158  template <typename Range, typename Domain, typename T>
1160  {
1161  template <typename C>
1162  static std::false_type
1163  test(...);
1164 
1165  template <typename C>
1166  static auto
1167  test(Range *r, Domain *d)
1168  -> decltype(std::declval<C>().vmult_add(*r, *d),
1169  std::declval<C>().Tvmult_add(*d, *r),
1170  std::true_type());
1171 
1172  public:
1173  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
1174  // otherwise it is std::false_type
1175 
1176  using type = decltype(test<T>(nullptr, nullptr));
1177  };
1178 
1179 
1180  // A helper function to apply a given vmult, or Tvmult to a vector with
1181  // intermediate storage
1182  template <typename Function, typename Range, typename Domain>
1183  void
1185  Range & v,
1186  const Domain &u,
1187  bool add)
1188  {
1189  GrowingVectorMemory<Range> vector_memory;
1190 
1191  typename VectorMemory<Range>::Pointer i(vector_memory);
1192  i->reinit(v, /*bool omit_zeroing_entries =*/true);
1193 
1194  function(*i, u);
1195 
1196  if (add)
1197  v += *i;
1198  else
1199  v = *i;
1200  }
1201 
1202 
1203  // A helper class to add a reduced matrix interface to a LinearOperator
1204  // (typically provided by Preconditioner classes)
1205  template <typename Range, typename Domain, typename Payload>
1207  {
1208  public:
1209  template <typename Matrix>
1210  void
1212  const Matrix & matrix)
1213  {
1214  op.vmult = [&matrix](Range &v, const Domain &u) {
1215  if (PointerComparison::equal(&v, &u))
1216  {
1217  // If v and u are the same memory location use intermediate
1218  // storage
1220  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1221  v,
1222  u,
1223  /*bool add =*/false);
1224  }
1225  else
1226  {
1227  matrix.vmult(v, u);
1228  }
1229  };
1230 
1231  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1232  // use intermediate storage to implement vmult_add with vmult
1234  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1235  v,
1236  u,
1237  /*bool add =*/true);
1238  };
1239 
1240  op.Tvmult = [&matrix](Domain &v, const Range &u) {
1241  if (PointerComparison::equal(&v, &u))
1242  {
1243  // If v and u are the same memory location use intermediate
1244  // storage
1246  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1247  v,
1248  u,
1249  /*bool add =*/false);
1250  }
1251  else
1252  {
1253  matrix.Tvmult(v, u);
1254  }
1255  };
1256 
1257  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1258  // use intermediate storage to implement Tvmult_add with Tvmult
1260  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1261  v,
1262  u,
1263  /*bool add =*/true);
1264  };
1265  }
1266  };
1267 
1268 
1269  // A helper class to add the full matrix interface to a LinearOperator
1270  template <typename Range, typename Domain, typename Payload>
1272  {
1273  public:
1274  template <typename Matrix>
1275  void
1277  const Matrix & matrix)
1278  {
1279  // As above ...
1280 
1282  op, matrix);
1283 
1284  // ... but add native vmult_add and Tvmult_add variants:
1285 
1286  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1287  if (PointerComparison::equal(&v, &u))
1288  {
1290  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1291  v,
1292  u,
1293  /*bool add =*/true);
1294  }
1295  else
1296  {
1297  matrix.vmult_add(v, u);
1298  }
1299  };
1300 
1301  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1302  if (PointerComparison::equal(&v, &u))
1303  {
1305  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1306  v,
1307  u,
1308  /*bool add =*/true);
1309  }
1310  else
1311  {
1312  matrix.Tvmult_add(v, u);
1313  }
1314  };
1315  }
1316  };
1317  } // namespace LinearOperatorImplementation
1318 } // namespace internal
1319 
1320 
1377 template <typename Range, typename Domain, typename Payload, typename Matrix>
1379 linear_operator(const Matrix &matrix)
1380 {
1381  // implement with the more generic variant below...
1382  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1383  matrix);
1384 }
1385 
1386 
1401 template <typename Range,
1402  typename Domain,
1403  typename Payload,
1404  typename OperatorExemplar,
1405  typename Matrix>
1407 linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
1408 {
1410  // Initialize the payload based on the input exemplar matrix
1412  Payload(operator_exemplar, matrix)};
1413 
1414  // Always store a reference to matrix and operator_exemplar in the lambda
1415  // functions. This ensures that a modification of the matrix after the
1416  // creation of a LinearOperator wrapper is respected - further a matrix
1417  // or an operator_exemplar cannot usually be copied...
1418 
1419  return_op.reinit_range_vector =
1420  [&operator_exemplar](Range &v, bool omit_zeroing_entries) {
1422  Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1423  };
1424 
1425  return_op.reinit_domain_vector = [&operator_exemplar](
1426  Domain &v, bool omit_zeroing_entries) {
1428  Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1429  };
1430 
1431  typename std::conditional<
1435  .
1436  operator()(return_op, matrix);
1437 
1438  return return_op;
1439 }
1440 
1441 
1442 
1459 template <typename Range, typename Domain, typename Payload, typename Matrix>
1462  const Matrix & matrix)
1463 {
1465  // Initialize the payload based on the LinearOperator exemplar
1466  auto return_op = operator_exemplar;
1467 
1468  typename std::conditional<
1472  .
1473  operator()(return_op, matrix);
1474 
1475  return return_op;
1476 }
1477 
1478 
1480 
1481 #ifndef DOXYGEN
1482 
1483 //
1484 // Ensure that we never capture a reference to a temporary by accident.
1485 // to avoid "stack use after free".
1486 //
1487 
1488 template <
1489  typename Range = Vector<double>,
1490  typename Domain = Range,
1492  typename OperatorExemplar,
1493  typename Matrix,
1494  typename =
1495  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1497 linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1498 
1499 template <
1500  typename Range = Vector<double>,
1501  typename Domain = Range,
1503  typename OperatorExemplar,
1504  typename Matrix,
1505  typename = typename std::enable_if<
1506  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1507  typename = typename std::enable_if<
1508  !std::is_same<OperatorExemplar,
1511 linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1512 
1513 template <
1514  typename Range = Vector<double>,
1515  typename Domain = Range,
1517  typename OperatorExemplar,
1518  typename Matrix,
1519  typename =
1520  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type,
1521  typename = typename std::enable_if<
1522  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1523  typename = typename std::enable_if<
1524  !std::is_same<OperatorExemplar,
1527 linear_operator(OperatorExemplar &&, 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>
1538  Matrix &&) = delete;
1539 
1540 template <
1541  typename Range = Vector<double>,
1542  typename Domain = Range,
1544  typename Matrix,
1545  typename =
1546  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1548 linear_operator(Matrix &&) = delete;
1549 
1550 template <typename Payload,
1551  typename Solver,
1552  typename Preconditioner,
1553  typename Range = typename Solver::vector_type,
1554  typename Domain = Range,
1555  typename = typename std::enable_if<
1556  !std::is_lvalue_reference<Preconditioner>::value>::type,
1557  typename = typename std::enable_if<
1558  !std::is_same<Preconditioner, PreconditionIdentity>::value>::type,
1559  typename = typename std::enable_if<
1560  !std::is_same<Preconditioner,
1564  Solver &,
1565  Preconditioner &&) = delete;
1566 
1567 #endif // DOXYGEN
1568 
1570 
1571 #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:1466
std::function< void(Range &v, const Domain &u)> vmult
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:394
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:393
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 &)