Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\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 - 2021 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 
1149  inline EmptyPayload
1151  {
1152  return {};
1153  }
1154 
1155 
1156 
1157  // A trait class that determines whether type T provides public
1158  // (templated or non-templated) vmult_add member functions
1159  template <typename Range, typename Domain, typename T>
1161  {
1162  template <typename C>
1163  static std::false_type
1164  test(...);
1165 
1166  template <typename C>
1167  static auto
1168  test(Range *r, Domain *d)
1169  -> decltype(std::declval<C>().vmult_add(*r, *d),
1170  std::declval<C>().Tvmult_add(*d, *r),
1171  std::true_type());
1172 
1173  public:
1174  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
1175  // otherwise it is std::false_type
1176 
1177  using type = decltype(test<T>(nullptr, nullptr));
1178  };
1179 
1180 
1181  // A helper function to apply a given vmult, or Tvmult to a vector with
1182  // intermediate storage
1183  template <typename Function, typename Range, typename Domain>
1184  void
1186  Range & v,
1187  const Domain &u,
1188  bool add)
1189  {
1190  GrowingVectorMemory<Range> vector_memory;
1191 
1192  typename VectorMemory<Range>::Pointer i(vector_memory);
1193  i->reinit(v, /*bool omit_zeroing_entries =*/true);
1194 
1195  function(*i, u);
1196 
1197  if (add)
1198  v += *i;
1199  else
1200  v = *i;
1201  }
1202 
1203 
1204  // A helper class to add a reduced matrix interface to a LinearOperator
1205  // (typically provided by Preconditioner classes)
1206  template <typename Range, typename Domain, typename Payload>
1208  {
1209  public:
1210  template <typename Matrix>
1211  void
1213  const Matrix & matrix)
1214  {
1215  op.vmult = [&matrix](Range &v, const Domain &u) {
1216  if (PointerComparison::equal(&v, &u))
1217  {
1218  // If v and u are the same memory location use intermediate
1219  // storage
1221  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1222  v,
1223  u,
1224  /*bool add =*/false);
1225  }
1226  else
1227  {
1228  matrix.vmult(v, u);
1229  }
1230  };
1231 
1232  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1233  // use intermediate storage to implement vmult_add with vmult
1235  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1236  v,
1237  u,
1238  /*bool add =*/true);
1239  };
1240 
1241  op.Tvmult = [&matrix](Domain &v, const Range &u) {
1242  if (PointerComparison::equal(&v, &u))
1243  {
1244  // If v and u are the same memory location use intermediate
1245  // storage
1247  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1248  v,
1249  u,
1250  /*bool add =*/false);
1251  }
1252  else
1253  {
1254  matrix.Tvmult(v, u);
1255  }
1256  };
1257 
1258  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1259  // use intermediate storage to implement Tvmult_add with Tvmult
1261  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1262  v,
1263  u,
1264  /*bool add =*/true);
1265  };
1266  }
1267  };
1268 
1269 
1270  // A helper class to add the full matrix interface to a LinearOperator
1271  template <typename Range, typename Domain, typename Payload>
1273  {
1274  public:
1275  template <typename Matrix>
1276  void
1278  const Matrix & matrix)
1279  {
1280  // As above ...
1281 
1283  op, matrix);
1284 
1285  // ... but add native vmult_add and Tvmult_add variants:
1286 
1287  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1288  if (PointerComparison::equal(&v, &u))
1289  {
1291  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1292  v,
1293  u,
1294  /*bool add =*/true);
1295  }
1296  else
1297  {
1298  matrix.vmult_add(v, u);
1299  }
1300  };
1301 
1302  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1303  if (PointerComparison::equal(&v, &u))
1304  {
1306  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1307  v,
1308  u,
1309  /*bool add =*/true);
1310  }
1311  else
1312  {
1313  matrix.Tvmult_add(v, u);
1314  }
1315  };
1316  }
1317  };
1318  } // namespace LinearOperatorImplementation
1319 } // namespace internal
1320 
1321 
1378 template <typename Range, typename Domain, typename Payload, typename Matrix>
1380 linear_operator(const Matrix &matrix)
1381 {
1382  // implement with the more generic variant below...
1383  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1384  matrix);
1385 }
1386 
1387 
1402 template <typename Range,
1403  typename Domain,
1404  typename Payload,
1405  typename OperatorExemplar,
1406  typename Matrix>
1408 linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
1409 {
1411  // Initialize the payload based on the input exemplar matrix
1413  Payload(operator_exemplar, matrix)};
1414 
1415  // Always store a reference to matrix and operator_exemplar in the lambda
1416  // functions. This ensures that a modification of the matrix after the
1417  // creation of a LinearOperator wrapper is respected - further a matrix
1418  // or an operator_exemplar cannot usually be copied...
1419 
1420  return_op.reinit_range_vector =
1421  [&operator_exemplar](Range &v, bool omit_zeroing_entries) {
1423  Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1424  };
1425 
1426  return_op.reinit_domain_vector = [&operator_exemplar](
1427  Domain &v, bool omit_zeroing_entries) {
1429  Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1430  };
1431 
1432  typename std::conditional<
1436  .
1437  operator()(return_op, matrix);
1438 
1439  return return_op;
1440 }
1441 
1442 
1443 
1460 template <typename Range, typename Domain, typename Payload, typename Matrix>
1463  const Matrix & matrix)
1464 {
1466  // Initialize the payload based on the LinearOperator exemplar
1467  auto return_op = operator_exemplar;
1468 
1469  typename std::conditional<
1473  .
1474  operator()(return_op, matrix);
1475 
1476  return return_op;
1477 }
1478 
1479 
1481 
1482 #ifndef DOXYGEN
1483 
1484 //
1485 // Ensure that we never capture a reference to a temporary by accident.
1486 // to avoid "stack use after free".
1487 //
1488 
1489 template <
1490  typename Range = Vector<double>,
1491  typename Domain = Range,
1493  typename OperatorExemplar,
1494  typename Matrix,
1495  typename =
1496  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1498 linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1499 
1500 template <
1501  typename Range = Vector<double>,
1502  typename Domain = Range,
1504  typename OperatorExemplar,
1505  typename Matrix,
1506  typename = typename std::enable_if<
1507  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1508  typename = typename std::enable_if<
1509  !std::is_same<OperatorExemplar,
1512 linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1513 
1514 template <
1515  typename Range = Vector<double>,
1516  typename Domain = Range,
1518  typename OperatorExemplar,
1519  typename Matrix,
1520  typename =
1521  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type,
1522  typename = typename std::enable_if<
1523  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1524  typename = typename std::enable_if<
1525  !std::is_same<OperatorExemplar,
1528 linear_operator(OperatorExemplar &&, Matrix &&) = delete;
1529 
1530 template <
1531  typename Range = Vector<double>,
1532  typename Domain = Range,
1534  typename Matrix,
1535  typename =
1536  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1539  Matrix &&) = delete;
1540 
1541 template <
1542  typename Range = Vector<double>,
1543  typename Domain = Range,
1545  typename Matrix,
1546  typename =
1547  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1549 linear_operator(Matrix &&) = delete;
1550 
1551 template <typename Payload,
1552  typename Solver,
1553  typename Preconditioner,
1554  typename Range = typename Solver::vector_type,
1555  typename Domain = Range,
1556  typename = typename std::enable_if<
1557  !std::is_lvalue_reference<Preconditioner>::value>::type,
1558  typename = typename std::enable_if<
1559  !std::is_same<Preconditioner, PreconditionIdentity>::value>::type,
1560  typename = typename std::enable_if<
1561  !std::is_same<Preconditioner,
1565  Solver &,
1566  Preconditioner &&) = delete;
1567 
1568 #endif // DOXYGEN
1569 
1571 
1572 #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:1461
std::function< void(Range &v, const Domain &u)> vmult
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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:400
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 &)