Reference documentation for deal.II version Git 5f32ffafa8 2019-10-18 00:26:50 -0600
\(\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 #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 
168 template <typename Range, typename Domain, typename Payload>
169 class LinearOperator : public Payload
170 {
171 public:
180  LinearOperator(const Payload &payload = Payload())
181  : Payload(payload)
182  , is_null_operator(false)
183  {
184  vmult = [](Range &, const Domain &) {
185  Assert(false,
186  ExcMessage("Uninitialized LinearOperator<Range, "
187  "Domain>::vmult called"));
188  };
189 
190  vmult_add = [](Range &, const Domain &) {
191  Assert(false,
192  ExcMessage("Uninitialized LinearOperator<Range, "
193  "Domain>::vmult_add called"));
194  };
195 
196  Tvmult = [](Domain &, const Range &) {
197  Assert(false,
198  ExcMessage("Uninitialized LinearOperator<Range, "
199  "Domain>::Tvmult called"));
200  };
201 
202  Tvmult_add = [](Domain &, const Range &) {
203  Assert(false,
204  ExcMessage("Uninitialized LinearOperator<Range, "
205  "Domain>::Tvmult_add called"));
206  };
207 
208  reinit_range_vector = [](Range &, bool) {
209  Assert(false,
210  ExcMessage("Uninitialized LinearOperator<Range, "
211  "Domain>::reinit_range_vector method called"));
212  };
213 
214  reinit_domain_vector = [](Domain &, bool) {
215  Assert(false,
216  ExcMessage("Uninitialized LinearOperator<Range, "
217  "Domain>::reinit_domain_vector method called"));
218  };
219  }
220 
225 
231  template <typename Op,
232  typename = typename std::enable_if<
233  !std::is_base_of<LinearOperator<Range, Domain, Payload>,
234  Op>::value>::type>
235  LinearOperator(const Op &op)
236  {
237  *this = linear_operator<Range, Domain, Payload, Op>(op);
238  }
239 
244  operator=(const LinearOperator<Range, Domain, Payload> &) = default;
245 
250  template <typename Op,
251  typename = typename std::enable_if<
252  !std::is_base_of<LinearOperator<Range, Domain, Payload>,
253  Op>::value>::type>
255  operator=(const Op &op)
256  {
257  *this = linear_operator<Range, Domain, Payload, Op>(op);
258  return *this;
259  }
260 
265  std::function<void(Range &v, const Domain &u)> vmult;
266 
271  std::function<void(Range &v, const Domain &u)> vmult_add;
272 
277  std::function<void(Domain &v, const Range &u)> Tvmult;
278 
283  std::function<void(Domain &v, const Range &u)> Tvmult_add;
284 
292  std::function<void(Range &v, bool omit_zeroing_entries)> reinit_range_vector;
293 
301  std::function<void(Domain &v, bool omit_zeroing_entries)>
303 
308 
315  {
316  *this = *this + second_op;
317  return *this;
318  }
319 
326  {
327  *this = *this - second_op;
328  return *this;
329  }
330 
337  {
338  *this = *this * second_op;
339  return *this;
340  }
341 
347  operator*=(typename Domain::value_type number)
348  {
349  *this = *this * number;
350  return *this;
351  }
352 
359 
361 };
362 
363 
368 
378 template <typename Range, typename Domain, typename Payload>
382 {
383  if (first_op.is_null_operator)
384  {
385  return second_op;
386  }
387  else if (second_op.is_null_operator)
388  {
389  return first_op;
390  }
391  else
392  {
394  static_cast<const Payload &>(first_op) +
395  static_cast<const Payload &>(second_op)};
396 
397  return_op.reinit_range_vector = first_op.reinit_range_vector;
398  return_op.reinit_domain_vector = first_op.reinit_domain_vector;
399 
400  // ensure to have valid computation objects by catching first_op and
401  // second_op by value
402 
403  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
404  first_op.vmult(v, u);
405  second_op.vmult_add(v, u);
406  };
407 
408  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
409  first_op.vmult_add(v, u);
410  second_op.vmult_add(v, u);
411  };
412 
413  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
414  second_op.Tvmult(v, u);
415  first_op.Tvmult_add(v, u);
416  };
417 
418  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
419  second_op.Tvmult_add(v, u);
420  first_op.Tvmult_add(v, u);
421  };
422 
423  return return_op;
424  }
425 }
426 
427 
437 template <typename Range, typename Domain, typename Payload>
441 {
442  if (first_op.is_null_operator)
443  {
444  return -1. * second_op;
445  }
446  else if (second_op.is_null_operator)
447  {
448  return first_op;
449  }
450  else
451  {
452  // implement with addition and scalar multiplication
453  return first_op + (-1. * second_op);
454  }
455 }
456 
457 
475 template <typename Range, typename Domain, typename Payload>
477 operator*(typename Range::value_type number,
479 {
480  static_assert(
481  std::is_convertible<typename Range::value_type,
482  typename Domain::value_type>::value,
483  "Range and Domain must have implicitly convertible 'value_type's");
484 
485  if (op.is_null_operator)
486  {
487  return op;
488  }
489  else if (number == 0.)
490  {
491  return null_operator(op);
492  }
493  else
494  {
496 
497  // ensure to have valid computation objects by catching number and op by
498  // value
499 
500  return_op.vmult = [number, op](Range &v, const Domain &u) {
501  op.vmult(v, u);
502  v *= number;
503  };
504 
505  return_op.vmult_add = [number, op](Range &v, const Domain &u) {
506  v /= number;
507  op.vmult_add(v, u);
508  v *= number;
509  };
510 
511  return_op.Tvmult = [number, op](Domain &v, const Range &u) {
512  op.Tvmult(v, u);
513  v *= number;
514  };
515 
516  return_op.Tvmult_add = [number, op](Domain &v, const Range &u) {
517  v /= number;
518  op.Tvmult_add(v, u);
519  v *= number;
520  };
521 
522  return return_op;
523  }
524 }
525 
526 
542 template <typename Range, typename Domain, typename Payload>
545  typename Domain::value_type number)
546 {
547  static_assert(
548  std::is_convertible<typename Range::value_type,
549  typename Domain::value_type>::value,
550  "Range and Domain must have implicitly convertible 'value_type's");
551 
552  return number * op;
553 }
554 
556 
557 
562 
572 template <typename Range,
573  typename Intermediate,
574  typename Domain,
575  typename Payload>
579 {
580  if (first_op.is_null_operator || second_op.is_null_operator)
581  {
583  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
584  return_op.reinit_range_vector = first_op.reinit_range_vector;
585  return null_operator(return_op);
586  }
587  else
588  {
590  static_cast<const Payload &>(first_op) *
591  static_cast<const Payload &>(second_op)};
592 
593  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
594  return_op.reinit_range_vector = first_op.reinit_range_vector;
595 
596  // ensure to have valid computation objects by catching first_op and
597  // second_op by value
598 
599  return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
600  GrowingVectorMemory<Intermediate> vector_memory;
601 
602  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
603  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
604  second_op.vmult(*i, u);
605  first_op.vmult(v, *i);
606  };
607 
608  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
609  GrowingVectorMemory<Intermediate> vector_memory;
610 
611  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
612  second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
613  second_op.vmult(*i, u);
614  first_op.vmult_add(v, *i);
615  };
616 
617  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
618  GrowingVectorMemory<Intermediate> vector_memory;
619 
620  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
621  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
622  first_op.Tvmult(*i, u);
623  second_op.Tvmult(v, *i);
624  };
625 
626  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
627  GrowingVectorMemory<Intermediate> vector_memory;
628 
629  typename VectorMemory<Intermediate>::Pointer i(vector_memory);
630  first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
631  first_op.Tvmult(*i, u);
632  second_op.Tvmult_add(v, *i);
633  };
634 
635  return return_op;
636  }
637 }
638 
639 
649 template <typename Range, typename Domain, typename Payload>
652 {
653  LinearOperator<Domain, Range, Payload> return_op{op.transpose_payload()};
654 
656  return_op.reinit_domain_vector = op.reinit_range_vector;
657 
658  return_op.vmult = op.Tvmult;
659  return_op.vmult_add = op.Tvmult_add;
660  return_op.Tvmult = op.vmult;
661  return_op.Tvmult_add = op.vmult_add;
662 
663  return return_op;
664 }
665 
666 
687 template <typename Payload,
688  typename Solver,
689  typename Preconditioner,
690  typename Range = typename Solver::vector_type,
691  typename Domain = Range>
694  Solver & solver,
695  const Preconditioner & preconditioner)
696 {
698  op.inverse_payload(solver, preconditioner)};
699 
701  return_op.reinit_domain_vector = op.reinit_range_vector;
702 
703  return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
704  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
705  solver.solve(op, v, u, preconditioner);
706  };
707 
708  return_op.vmult_add = [op, &solver, &preconditioner](Range & v,
709  const Domain &u) {
710  GrowingVectorMemory<Range> vector_memory;
711 
712  typename VectorMemory<Range>::Pointer v2(vector_memory);
713  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
714  solver.solve(op, *v2, u, preconditioner);
715  v += *v2;
716  };
717 
718  return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
719  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
720  solver.solve(transpose_operator(op), v, u, preconditioner);
721  };
722 
723  return_op.Tvmult_add = [op, &solver, &preconditioner](Range & v,
724  const Domain &u) {
725  GrowingVectorMemory<Range> vector_memory;
726 
727  typename VectorMemory<Range>::Pointer v2(vector_memory);
728  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
729  solver.solve(transpose_operator(op), *v2, u, preconditioner);
730  v += *v2;
731  };
732 
733  return return_op;
734 }
735 
736 
745 template <typename Payload,
746  typename Solver,
747  typename Range = typename Solver::vector_type,
748  typename Domain = Range>
751  Solver & solver,
752  const LinearOperator<Range, Domain, Payload> &preconditioner)
753 {
755  op.inverse_payload(solver, preconditioner)};
756 
758  return_op.reinit_domain_vector = op.reinit_range_vector;
759 
760  return_op.vmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
761  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
762  solver.solve(op, v, u, preconditioner);
763  };
764 
765  return_op.vmult_add = [op, &solver, preconditioner](Range & v,
766  const Domain &u) {
767  GrowingVectorMemory<Range> vector_memory;
768 
769  typename VectorMemory<Range>::Pointer v2(vector_memory);
770  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
771  solver.solve(op, *v2, u, preconditioner);
772  v += *v2;
773  };
774 
775  return_op.Tvmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
776  op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
777  solver.solve(transpose_operator(op), v, u, preconditioner);
778  };
779 
780  return_op.Tvmult_add = [op, &solver, preconditioner](Range & v,
781  const Domain &u) {
782  GrowingVectorMemory<Range> vector_memory;
783 
784  typename VectorMemory<Range>::Pointer v2(vector_memory);
785  op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
786  solver.solve(transpose_operator(op), *v2, u, preconditioner);
787  v += *v2;
788  };
789 
790  return return_op;
791 }
792 
793 
803 template <typename Payload,
804  typename Solver,
805  typename Range = typename Solver::vector_type,
806  typename Domain = Range>
809  Solver & solver)
810 {
811  return inverse_operator(op, solver, identity_operator(op));
812 }
813 
814 
823 template <typename Payload,
824  typename Solver,
825  typename Range = typename Solver::vector_type,
826  typename Domain = Range>
829  Solver & solver,
830  const PreconditionIdentity &)
831 {
832  return inverse_operator(op, solver);
833 }
834 
836 
837 
842 
854 template <
855  typename Range,
858 identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
859 {
860  LinearOperator<Range, Range, Payload> return_op{Payload()};
861 
862  return_op.reinit_range_vector = reinit_vector;
863  return_op.reinit_domain_vector = reinit_vector;
864 
865  return_op.vmult = [](Range &v, const Range &u) { v = u; };
866 
867  return_op.vmult_add = [](Range &v, const Range &u) { v += u; };
868 
869  return_op.Tvmult = [](Range &v, const Range &u) { v = u; };
870 
871  return_op.Tvmult_add = [](Range &v, const Range &u) { v += u; };
872 
873  return return_op;
874 }
875 
876 
889 template <typename Range, typename Domain, typename Payload>
892 {
893  auto return_op = identity_operator<Range, Payload>(op.reinit_range_vector);
894  static_cast<Payload &>(return_op) = op.identity_payload();
895 
896  return return_op;
897 }
898 
899 
909 template <typename Range, typename Domain, typename Payload>
912 {
913  LinearOperator<Range, Domain, Payload> return_op{op.null_payload()};
914 
915  return_op.is_null_operator = true;
916 
917  return_op.reinit_range_vector = op.reinit_range_vector;
918  return_op.reinit_domain_vector = op.reinit_domain_vector;
919 
920  return_op.vmult = [](Range &v, const Domain &) { v = 0.; };
921 
922  return_op.vmult_add = [](Range &, const Domain &) {};
923 
924  return_op.Tvmult = [](Domain &v, const Range &) { v = 0.; };
925 
926  return_op.Tvmult_add = [](Domain &, const Range &) {};
927 
928  return return_op;
929 }
930 
931 
944 template <
945  typename Range,
948 mean_value_filter(const std::function<void(Range &, bool)> &reinit_vector)
949 {
950  LinearOperator<Range, Range, Payload> return_op{Payload()};
951 
952  return_op.reinit_range_vector = reinit_vector;
953  return_op.reinit_domain_vector = reinit_vector;
954 
955  return_op.vmult = [](Range &v, const Range &u) {
956  const auto mean = u.mean_value();
957 
958  v = u;
959  v.add(-mean);
960  };
961 
962  return_op.vmult_add = [](Range &v, const Range &u) {
963  const auto mean = u.mean_value();
964 
965  v += u;
966  v.add(-mean);
967  };
968 
969  return_op.Tvmult = return_op.vmult_add;
970  return_op.Tvmult_add = return_op.vmult_add;
971 
972  return return_op;
973 }
974 
975 
988 template <typename Range, typename Domain, typename Payload>
991 {
992  auto return_op = mean_value_filter<Range, Payload>(op.reinit_range_vector);
993  static_cast<Payload &>(return_op) = op.identity_payload();
994 
995  return return_op;
996 }
997 
998 
999 namespace internal
1000 {
1001  namespace LinearOperatorImplementation
1002  {
1014  template <typename Vector>
1015  class ReinitHelper
1016  {
1017  public:
1029  template <typename Matrix>
1030  static void
1031  reinit_range_vector(const Matrix &matrix,
1032  Vector & v,
1033  bool omit_zeroing_entries)
1034  {
1035  v.reinit(matrix.m(), omit_zeroing_entries);
1036  }
1037 
1049  template <typename Matrix>
1050  static void
1051  reinit_domain_vector(const Matrix &matrix,
1052  Vector & v,
1053  bool omit_zeroing_entries)
1054  {
1055  v.reinit(matrix.n(), omit_zeroing_entries);
1056  }
1057  };
1058 
1059 
1073  {
1074  public:
1082  template <typename... Args>
1083  EmptyPayload(const Args &...)
1084  {}
1085 
1086 
1090  EmptyPayload
1092  {
1093  return *this;
1094  }
1095 
1096 
1100  EmptyPayload
1102  {
1103  return *this;
1104  }
1105 
1106 
1110  EmptyPayload
1112  {
1113  return *this;
1114  }
1115 
1116 
1120  template <typename Solver, typename Preconditioner>
1121  EmptyPayload
1122  inverse_payload(Solver &, const Preconditioner &) const
1123  {
1124  return *this;
1125  }
1126  };
1127 
1132  inline EmptyPayload
1133  operator+(const EmptyPayload &, const EmptyPayload &)
1134  {
1135  return {};
1136  }
1137 
1142  inline EmptyPayload operator*(const EmptyPayload &, const EmptyPayload &)
1143  {
1144  return {};
1145  }
1146 
1147 
1148 
1149  // A trait class that determines whether type T provides public
1150  // (templated or non-templated) vmult_add member functions
1151  template <typename Range, typename Domain, typename T>
1152  class has_vmult_add_and_Tvmult_add
1153  {
1154  template <typename C>
1155  static std::false_type
1156  test(...);
1157 
1158  template <typename C>
1159  static auto
1160  test(Range *r, Domain *d)
1161  -> decltype(std::declval<C>().vmult_add(*r, *d),
1162  std::declval<C>().Tvmult_add(*d, *r),
1163  std::true_type());
1164 
1165  public:
1166  // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
1167  // otherwise it is std::false_type
1168 
1169  using type = decltype(test<T>(nullptr, nullptr));
1170  };
1171 
1172 
1173  // A helper function to apply a given vmult, or Tvmult to a vector with
1174  // intermediate storage
1175  template <typename Function, typename Range, typename Domain>
1176  void
1177  apply_with_intermediate_storage(Function function,
1178  Range & v,
1179  const Domain &u,
1180  bool add)
1181  {
1182  GrowingVectorMemory<Range> vector_memory;
1183 
1184  typename VectorMemory<Range>::Pointer i(vector_memory);
1185  i->reinit(v, /*bool omit_zeroing_entries =*/true);
1186 
1187  function(*i, u);
1188 
1189  if (add)
1190  v += *i;
1191  else
1192  v = *i;
1193  }
1194 
1195 
1196  // A helper class to add a reduced matrix interface to a LinearOperator
1197  // (typically provided by Preconditioner classes)
1198  template <typename Range, typename Domain, typename Payload>
1199  class MatrixInterfaceWithoutVmultAdd
1200  {
1201  public:
1202  template <typename Matrix>
1203  void
1205  const Matrix & matrix)
1206  {
1207  op.vmult = [&matrix](Range &v, const Domain &u) {
1208  if (PointerComparison::equal(&v, &u))
1209  {
1210  // If v and u are the same memory location use intermediate
1211  // storage
1212  apply_with_intermediate_storage(
1213  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1214  v,
1215  u,
1216  /*bool add =*/false);
1217  }
1218  else
1219  {
1220  matrix.vmult(v, u);
1221  }
1222  };
1223 
1224  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1225  // use intermediate storage to implement vmult_add with vmult
1226  apply_with_intermediate_storage(
1227  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1228  v,
1229  u,
1230  /*bool add =*/true);
1231  };
1232 
1233  op.Tvmult = [&matrix](Domain &v, const Range &u) {
1234  if (PointerComparison::equal(&v, &u))
1235  {
1236  // If v and u are the same memory location use intermediate
1237  // storage
1238  apply_with_intermediate_storage(
1239  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1240  v,
1241  u,
1242  /*bool add =*/false);
1243  }
1244  else
1245  {
1246  matrix.Tvmult(v, u);
1247  }
1248  };
1249 
1250  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1251  // use intermediate storage to implement Tvmult_add with Tvmult
1252  apply_with_intermediate_storage(
1253  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1254  v,
1255  u,
1256  /*bool add =*/true);
1257  };
1258  }
1259  };
1260 
1261 
1262  // A helper class to add the full matrix interface to a LinearOperator
1263  template <typename Range, typename Domain, typename Payload>
1264  class MatrixInterfaceWithVmultAdd
1265  {
1266  public:
1267  template <typename Matrix>
1268  void
1270  const Matrix & matrix)
1271  {
1272  // As above ...
1273 
1274  MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>().operator()(
1275  op, matrix);
1276 
1277  // ... but add native vmult_add and Tvmult_add variants:
1278 
1279  op.vmult_add = [&matrix](Range &v, const Domain &u) {
1280  if (PointerComparison::equal(&v, &u))
1281  {
1282  apply_with_intermediate_storage(
1283  [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1284  v,
1285  u,
1286  /*bool add =*/true);
1287  }
1288  else
1289  {
1290  matrix.vmult_add(v, u);
1291  }
1292  };
1293 
1294  op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1295  if (PointerComparison::equal(&v, &u))
1296  {
1297  apply_with_intermediate_storage(
1298  [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1299  v,
1300  u,
1301  /*bool add =*/true);
1302  }
1303  else
1304  {
1305  matrix.Tvmult_add(v, u);
1306  }
1307  };
1308  }
1309  };
1310  } // namespace LinearOperatorImplementation
1311 } // namespace internal
1312 
1313 
1371 template <typename Range, typename Domain, typename Payload, typename Matrix>
1373 linear_operator(const Matrix &matrix)
1374 {
1375  // implement with the more generic variant below...
1376  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1377  matrix);
1378 }
1379 
1380 
1396 template <typename Range,
1397  typename Domain,
1398  typename Payload,
1399  typename OperatorExemplar,
1400  typename Matrix>
1402 linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
1403 {
1405  // Initialize the payload based on the input exemplar matrix
1407  Payload(operator_exemplar, matrix)};
1408 
1409  // Always store a reference to matrix and operator_exemplar in the lambda
1410  // functions. This ensures that a modification of the matrix after the
1411  // creation of a LinearOperator wrapper is respected - further a matrix
1412  // or an operator_exemplar cannot usually be copied...
1413 
1414  return_op.reinit_range_vector =
1415  [&operator_exemplar](Range &v, bool omit_zeroing_entries) {
1417  Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1418  };
1419 
1420  return_op.reinit_domain_vector = [&operator_exemplar](
1421  Domain &v, bool omit_zeroing_entries) {
1423  Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1424  };
1425 
1426  typename std::conditional<
1427  has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
1428  MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
1429  MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>::type()
1430  .
1431  operator()(return_op, matrix);
1432 
1433  return return_op;
1434 }
1435 
1436 
1437 
1455 template <typename Range, typename Domain, typename Payload, typename Matrix>
1458  const Matrix & matrix)
1459 {
1461  // Initialize the payload based on the LinearOperator exemplar
1462  auto return_op = operator_exemplar;
1463 
1464  typename std::conditional<
1465  has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
1466  MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
1467  MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>::type()
1468  .
1469  operator()(return_op, matrix);
1470 
1471  return return_op;
1472 }
1473 
1474 
1476 
1477 #ifndef DOXYGEN
1478 
1479 //
1480 // Ensure that we never capture a reference to a temporary by accident.
1481 // to avoid "stack use after free".
1482 //
1483 
1484 template <
1485  typename Range = Vector<double>,
1486  typename Domain = Range,
1488  typename OperatorExemplar,
1489  typename Matrix,
1490  typename =
1491  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1493 linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1494 
1495 template <
1496  typename Range = Vector<double>,
1497  typename Domain = Range,
1499  typename OperatorExemplar,
1500  typename Matrix,
1501  typename = typename std::enable_if<
1502  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1503  typename = typename std::enable_if<
1504  !std::is_same<OperatorExemplar,
1507 linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1508 
1509 template <
1510  typename Range = Vector<double>,
1511  typename Domain = Range,
1513  typename OperatorExemplar,
1514  typename Matrix,
1515  typename =
1516  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type,
1517  typename = typename std::enable_if<
1518  !std::is_lvalue_reference<OperatorExemplar>::value>::type,
1519  typename = typename std::enable_if<
1520  !std::is_same<OperatorExemplar,
1523 linear_operator(OperatorExemplar &&, Matrix &&) = delete;
1524 
1525 template <
1526  typename Range = Vector<double>,
1527  typename Domain = Range,
1529  typename Matrix,
1530  typename =
1531  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1534  Matrix &&) = delete;
1535 
1536 template <
1537  typename Range = Vector<double>,
1538  typename Domain = Range,
1540  typename Matrix,
1541  typename =
1542  typename std::enable_if<!std::is_lvalue_reference<Matrix>::value>::type>
1544 linear_operator(Matrix &&) = delete;
1545 
1546 template <typename Payload,
1547  typename Solver,
1548  typename Preconditioner,
1549  typename Range = typename Solver::vector_type,
1550  typename Domain = Range,
1551  typename = typename std::enable_if<
1552  !std::is_lvalue_reference<Preconditioner>::value>::type,
1553  typename = typename std::enable_if<
1554  !std::is_same<Preconditioner, PreconditionIdentity>::value>::type,
1555  typename = typename std::enable_if<
1556  !std::is_same<Preconditioner,
1560  Solver &,
1561  Preconditioner &&) = delete;
1562 
1563 #endif // DOXYGEN
1564 
1565 DEAL_II_NAMESPACE_CLOSE
1566 
1567 #endif
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)
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< 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)
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, Domain, Payload > operator-(const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
VectorType vector_type
Definition: solver.h:337
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())
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
LinearOperator< Range, Domain, Payload > identity_operator(const LinearOperator< Range, Domain, Payload > &)
static bool equal(const T *p1, const T *p2)
LinearOperator(const Op &op)