Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
linear_operator.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2014 - 2019 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 
21 #include <deal.II/base/exceptions.h>
22 
23 #include <deal.II/lac/vector_memory.h>
24 
25 #include <array>
26 #include <functional>
27 #include <type_traits>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 // Forward declarations:
32 
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>
51 
52 template <
53  typename Range = Vector<double>,
54  typename Domain = Range,
56  typename OperatorExemplar,
57  typename Matrix>
59 linear_operator(const OperatorExemplar &, const Matrix &);
60 
61 template <
62  typename Range = Vector<double>,
63  typename Domain = Range,
65  typename Matrix>
67 linear_operator(const Matrix &);
68 
69 template <
70  typename Range = Vector<double>,
71  typename Domain = Range,
75 
76 template <typename Range, typename Domain, typename Payload>
79 
80 
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 
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 
648 template <typename Range, typename Domain, typename Payload>
651 {
652  LinearOperator<Domain, Range, Payload> return_op{op.transpose_payload()};
653 
655  return_op.reinit_domain_vector = op.reinit_range_vector;
656 
657  return_op.vmult = op.Tvmult;
658  return_op.vmult_add = op.Tvmult_add;
659  return_op.Tvmult = op.vmult;
660  return_op.Tvmult_add = op.vmult_add;
661 
662  return return_op;
663 }
664 
665 
686 template <typename Payload,
687  typename Solver,
688  typename Preconditioner,
689  typename Range = typename Solver::vector_type,
690  typename Domain = Range>
693  Solver & solver,
694  const Preconditioner & preconditioner)
695 {
697  op.inverse_payload(solver, preconditioner)};
698 
700  return_op.reinit_domain_vector = op.reinit_range_vector;
701 
702  return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
703  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
704  solver.solve(op, v, u, preconditioner);
705  };
706 
707  return_op.vmult_add = [op, &solver, &preconditioner](Range & v,
708  const Domain &u) {
709  GrowingVectorMemory<Range> vector_memory;
710 
711  typename VectorMemory<Range>::Pointer v2(vector_memory);
712  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
713  solver.solve(op, *v2, u, preconditioner);
714  v += *v2;
715  };
716 
717  return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
718  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
719  solver.solve(transpose_operator(op), v, u, preconditioner);
720  };
721 
722  return_op.Tvmult_add = [op, &solver, &preconditioner](Range & v,
723  const Domain &u) {
724  GrowingVectorMemory<Range> vector_memory;
725 
726  typename VectorMemory<Range>::Pointer v2(vector_memory);
727  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
728  solver.solve(transpose_operator(op), *v2, u, preconditioner);
729  v += *v2;
730  };
731 
732  return return_op;
733 }
734 
735 
744 template <typename Payload,
745  typename Solver,
746  typename Range = typename Solver::vector_type,
747  typename Domain = Range>
750  Solver & solver,
751  const LinearOperator<Range, Domain, Payload> &preconditioner)
752 {
754  op.inverse_payload(solver, preconditioner)};
755 
757  return_op.reinit_domain_vector = op.reinit_range_vector;
758 
759  return_op.vmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
760  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
761  solver.solve(op, v, u, preconditioner);
762  };
763 
764  return_op.vmult_add = [op, &solver, preconditioner](Range & v,
765  const Domain &u) {
766  GrowingVectorMemory<Range> vector_memory;
767 
768  typename VectorMemory<Range>::Pointer v2(vector_memory);
769  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
770  solver.solve(op, *v2, u, preconditioner);
771  v += *v2;
772  };
773 
774  return_op.Tvmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
775  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
776  solver.solve(transpose_operator(op), v, u, preconditioner);
777  };
778 
779  return_op.Tvmult_add = [op, &solver, preconditioner](Range & v,
780  const Domain &u) {
781  GrowingVectorMemory<Range> vector_memory;
782 
783  typename VectorMemory<Range>::Pointer v2(vector_memory);
784  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
785  solver.solve(transpose_operator(op), *v2, u, preconditioner);
786  v += *v2;
787  };
788 
789  return return_op;
790 }
791 
792 
802 template <typename Payload,
803  typename Solver,
804  typename Range = typename Solver::vector_type,
805  typename Domain = Range>
808  Solver & solver)
809 {
810  return inverse_operator(op, solver, identity_operator(op));
811 }
812 
813 
822 template <typename Payload,
823  typename Solver,
824  typename Range = typename Solver::vector_type,
825  typename Domain = Range>
828  Solver & solver,
829  const PreconditionIdentity &)
830 {
831  return inverse_operator(op, solver);
832 }
833 
835 
836 
841 
853 template <
854  typename Range,
857 identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
858 {
859  LinearOperator<Range, Range, Payload> return_op{Payload()};
860 
861  return_op.reinit_range_vector = reinit_vector;
862  return_op.reinit_domain_vector = reinit_vector;
863 
864  return_op.vmult = [](Range &v, const Range &u) { v = u; };
865 
866  return_op.vmult_add = [](Range &v, const Range &u) { v += u; };
867 
868  return_op.Tvmult = [](Range &v, const Range &u) { v = u; };
869 
870  return_op.Tvmult_add = [](Range &v, const Range &u) { v += u; };
871 
872  return return_op;
873 }
874 
875 
888 template <typename Range, typename Domain, typename Payload>
891 {
892  auto return_op = identity_operator<Range, Payload>(op.reinit_range_vector);
893  static_cast<Payload &>(return_op) = op.identity_payload();
894 
895  return return_op;
896 }
897 
898 
908 template <typename Range, typename Domain, typename Payload>
911 {
912  LinearOperator<Range, Domain, Payload> return_op{op.null_payload()};
913 
914  return_op.is_null_operator = true;
915 
916  return_op.reinit_range_vector = op.reinit_range_vector;
917  return_op.reinit_domain_vector = op.reinit_domain_vector;
918 
919  return_op.vmult = [](Range &v, const Domain &) { v = 0.; };
920 
921  return_op.vmult_add = [](Range &, const Domain &) {};
922 
923  return_op.Tvmult = [](Domain &v, const Range &) { v = 0.; };
924 
925  return_op.Tvmult_add = [](Domain &, const Range &) {};
926 
927  return return_op;
928 }
929 
930 
943 template <
944  typename Range,
947 mean_value_filter(const std::function<void(Range &, bool)> &reinit_vector)
948 {
949  LinearOperator<Range, Range, Payload> return_op{Payload()};
950 
951  return_op.reinit_range_vector = reinit_vector;
952  return_op.reinit_domain_vector = reinit_vector;
953 
954  return_op.vmult = [](Range &v, const Range &u) {
955  const auto mean = u.mean_value();
956 
957  v = u;
958  v.add(-mean);
959  };
960 
961  return_op.vmult_add = [](Range &v, const Range &u) {
962  const auto mean = u.mean_value();
963 
964  v += u;
965  v.add(-mean);
966  };
967 
968  return_op.Tvmult = return_op.vmult_add;
969  return_op.Tvmult_add = return_op.vmult_add;
970 
971  return return_op;
972 }
973 
974 
987 template <typename Range, typename Domain, typename Payload>
990 {
991  auto return_op = mean_value_filter<Range, Payload>(op.reinit_range_vector);
992  static_cast<Payload &>(return_op) = op.identity_payload();
993 
994  return return_op;
995 }
996 
997 
998 namespace internal
999 {
1000  namespace LinearOperatorImplementation
1001  {
1013  template <typename Vector>
1014  class ReinitHelper
1015  {
1016  public:
1028  template <typename Matrix>
1029  static void
1030  reinit_range_vector(const Matrix &matrix,
1031  Vector & v,
1032  bool omit_zeroing_entries)
1033  {
1034  v.reinit(matrix.m(), omit_zeroing_entries);
1035  }
1036 
1048  template <typename Matrix>
1049  static void
1050  reinit_domain_vector(const Matrix &matrix,
1051  Vector & v,
1052  bool omit_zeroing_entries)
1053  {
1054  v.reinit(matrix.n(), omit_zeroing_entries);
1055  }
1056  };
1057 
1058 
1072  {
1073  public:
1081  template <typename... Args>
1082  EmptyPayload(const Args &...)
1083  {}
1084 
1085 
1089  EmptyPayload
1091  {
1092  return *this;
1093  }
1094 
1095 
1099  EmptyPayload
1101  {
1102  return *this;
1103  }
1104 
1105 
1109  EmptyPayload
1111  {
1112  return *this;
1113  }
1114 
1115 
1119  template <typename Solver, typename Preconditioner>
1120  EmptyPayload
1121  inverse_payload(Solver &, const Preconditioner &) const
1122  {
1123  return *this;
1124  }
1125  };
1126 
1131  inline EmptyPayload
1132  operator+(const EmptyPayload &, const EmptyPayload &)
1133  {
1134  return {};
1135  }
1136 
1141  inline EmptyPayload operator*(const EmptyPayload &, const EmptyPayload &)
1142  {
1143  return {};
1144  }
1145 
1146 
1147 
1148  // A trait class that determines whether type T provides public
1149  // (templated or non-templated) vmult_add member functions
1150  template <typename Range, typename Domain, typename T>
1151  class has_vmult_add_and_Tvmult_add
1152  {
1153  template <typename C>
1154  static std::false_type
1155  test(...);
1156 
1157  template <typename C>
1158  static auto
1159  test(Range *r, Domain *d)
1160  -> decltype(std::declval<C>().vmult_add(*r, *d),
1161  std::declval<C>().Tvmult_add(*d, *r),
1162  std::true_type());
1163 
1164  public:
1165  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
1166  // otherwise it is std::false_type
1167 
1168  using type = decltype(test<T>(nullptr, nullptr));
1169  };
1170 
1171 
1172  // A helper function to apply a given vmult, or Tvmult to a vector with
1173  // intermediate storage
1174  template <typename Function, typename Range, typename Domain>
1175  void
1176  apply_with_intermediate_storage(Function function,
1177  Range & v,
1178  const Domain &u,
1179  bool add)
1180  {
1181  GrowingVectorMemory<Range> vector_memory;
1182 
1183  typename VectorMemory<Range>::Pointer i(vector_memory);
1184  i->reinit(v, /*bool omit_zeroing_entries =*/true);
1185 
1186  function(*i, u);
1187 
1188  if (add)
1189  v += *i;
1190  else
1191  v = *i;
1192  }
1193 
1194 
1195  // A helper class to add a reduced matrix interface to a LinearOperator
1196  // (typically provided by Preconditioner classes)
1197  template <typename Range, typename Domain, typename Payload>
1198  class MatrixInterfaceWithoutVmultAdd
1199  {
1200  public:
1201  template <typename Matrix>
1202  void
1204  const Matrix & matrix)
1205  {
1206  op.vmult = [&matrix](Range &v, const Domain &u) {
1207  if (PointerComparison::equal(&v, &u))
1208  {
1209  // If v and u are the same memory location use intermediate
1210  // storage
1211  apply_with_intermediate_storage(
1212  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1213  v,
1214  u,
1215  /*bool add =*/false);
1216  }
1217  else
1218  {
1219  matrix.vmult(v, u);
1220  }
1221  };
1222 
1223  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1224  // use intermediate storage to implement vmult_add with vmult
1225  apply_with_intermediate_storage(
1226  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1227  v,
1228  u,
1229  /*bool add =*/true);
1230  };
1231 
1232  op.Tvmult = [&matrix](Domain &v, const Range &u) {
1233  if (PointerComparison::equal(&v, &u))
1234  {
1235  // If v and u are the same memory location use intermediate
1236  // storage
1237  apply_with_intermediate_storage(
1238  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1239  v,
1240  u,
1241  /*bool add =*/false);
1242  }
1243  else
1244  {
1245  matrix.Tvmult(v, u);
1246  }
1247  };
1248 
1249  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1250  // use intermediate storage to implement Tvmult_add with Tvmult
1251  apply_with_intermediate_storage(
1252  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1253  v,
1254  u,
1255  /*bool add =*/true);
1256  };
1257  }
1258  };
1259 
1260 
1261  // A helper class to add the full matrix interface to a LinearOperator
1262  template <typename Range, typename Domain, typename Payload>
1263  class MatrixInterfaceWithVmultAdd
1264  {
1265  public:
1266  template <typename Matrix>
1267  void
1269  const Matrix & matrix)
1270  {
1271  // As above ...
1272 
1273  MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>().operator()(
1274  op, matrix);
1275 
1276  // ... but add native vmult_add and Tvmult_add variants:
1277 
1278  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1279  if (PointerComparison::equal(&v, &u))
1280  {
1281  apply_with_intermediate_storage(
1282  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1283  v,
1284  u,
1285  /*bool add =*/true);
1286  }
1287  else
1288  {
1289  matrix.vmult_add(v, u);
1290  }
1291  };
1292 
1293  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1294  if (PointerComparison::equal(&v, &u))
1295  {
1296  apply_with_intermediate_storage(
1297  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1298  v,
1299  u,
1300  /*bool add =*/true);
1301  }
1302  else
1303  {
1304  matrix.Tvmult_add(v, u);
1305  }
1306  };
1307  }
1308  };
1309  } // namespace LinearOperatorImplementation
1310 } // namespace internal
1311 
1312 
1370 template <typename Range, typename Domain, typename Payload, typename Matrix>
1372 linear_operator(const Matrix &matrix)
1373 {
1374  // implement with the more generic variant below...
1375  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1376  matrix);
1377 }
1378 
1379 
1395 template <typename Range,
1396  typename Domain,
1397  typename Payload,
1398  typename OperatorExemplar,
1399  typename Matrix>
1401 linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
1402 {
1404  // Initialize the payload based on the input exemplar matrix
1406  Payload(operator_exemplar, matrix)};
1407 
1408  // Always store a reference to matrix and operator_exemplar in the lambda
1409  // functions. This ensures that a modification of the matrix after the
1410  // creation of a LinearOperator wrapper is respected - further a matrix
1411  // or an operator_exemplar cannot usually be copied...
1412 
1413  return_op.reinit_range_vector =
1414  [&operator_exemplar](Range &v, bool omit_zeroing_entries) {
1416  Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1417  };
1418 
1419  return_op.reinit_domain_vector = [&operator_exemplar](
1420  Domain &v, bool omit_zeroing_entries) {
1422  Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1423  };
1424 
1425  typename std::conditional<
1426  has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
1427  MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
1428  MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>::type()
1429  .
1430  operator()(return_op, matrix);
1431 
1432  return return_op;
1433 }
1434 
1435 
1436 
1454 template <typename Range, typename Domain, typename Payload, typename Matrix>
1457  const Matrix & matrix)
1458 {
1460  // Initialize the payload based on the LinearOperator exemplar
1461  auto return_op = operator_exemplar;
1462 
1463  typename std::conditional<
1464  has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
1465  MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
1466  MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>::type()
1467  .
1468  operator()(return_op, matrix);
1469 
1470  return return_op;
1471 }
1472 
1473 
1475 
1476 #ifndef DOXYGEN
1477 
1478 //
1479 // Ensure that we never capture a reference to a temporary by accident.
1480 // to avoid "stack use after free".
1481 //
1482 
1483 template <
1484  typename Range = Vector<double>,
1485  typename Domain = Range,
1487  typename OperatorExemplar,
1488  typename Matrix,
1489  typename =
1490  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1492 linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1493 
1494 template <
1495  typename Range = Vector<double>,
1496  typename Domain = Range,
1498  typename OperatorExemplar,
1499  typename Matrix,
1500  typename = typename std::enable_if<
1501  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1502  typename = typename std::enable_if<
1503  !std::is_same<OperatorExemplar,
1506 linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1507 
1508 template <
1509  typename Range = Vector<double>,
1510  typename Domain = Range,
1512  typename OperatorExemplar,
1513  typename Matrix,
1514  typename =
1515  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type,
1516  typename = typename std::enable_if<
1517  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1518  typename = typename std::enable_if<
1519  !std::is_same<OperatorExemplar,
1522 linear_operator(OperatorExemplar &&, Matrix &&) = delete;
1523 
1524 template <
1525  typename Range = Vector<double>,
1526  typename Domain = Range,
1528  typename Matrix,
1529  typename =
1530  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1533  Matrix &&) = delete;
1534 
1535 template <
1536  typename Range = Vector<double>,
1537  typename Domain = Range,
1539  typename Matrix,
1540  typename =
1541  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1543 linear_operator(Matrix &&) = delete;
1544 
1545 template <typename Payload,
1546  typename Solver,
1547  typename Preconditioner,
1548  typename Range = typename Solver::vector_type,
1549  typename Domain = Range,
1550  typename = typename std::enable_if<
1551  !std::is_lvalue_reference<Preconditioner>::value>::type,
1552  typename = typename std::enable_if<
1553  !std::is_same<Preconditioner, PreconditionIdentity>::value>::type,
1554  typename = typename std::enable_if<
1555  !std::is_same<Preconditioner,
1559  Solver &,
1560  Preconditioner &&) = delete;
1561 
1562 #endif // DOXYGEN
1563 
1564 DEAL_II_NAMESPACE_CLOSE
1565 
1566 #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< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
LinearOperator< Range, Domain, Payload > & operator*=(const LinearOperator< Domain, Domain, Payload > &second_op)
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator+(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
EmptyPayload inverse_payload(Solver &, const Preconditioner &) const
LinearAlgebra::distributed::Vector< Number > Vector
static void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
LinearOperator< Range, Domain, Payload > & operator=(const LinearOperator< Range, Domain, Payload > &)=default
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &op)
SymmetricTensor< rank_, dim, typename ProductType< Number, OtherNumber >::type > operator-(const SymmetricTensor< rank_, dim, Number > &left, const SymmetricTensor< rank_, dim, OtherNumber > &right)
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
std::function< void(Range &v, const Domain &u)> vmult
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
LinearOperator< Range, Domain, Payload > operator*(typename Range::value_type number, const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Range, Payload > identity_operator(const std::function< void(Range &, bool)> &reinit_vector)
VectorType vector_type
Definition: solver.h:334
std::function< void(Range &v, const Domain &u)> vmult_add
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())
static bool equal(const T *p1, const T *p2)
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
LinearOperator(const Op &op)
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const