Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
linear_operator.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2014 - 2022 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
33namespace internal
34{
35 namespace LinearOperatorImplementation
36 {
37 class EmptyPayload;
38 }
39} // namespace internal
40
41template <typename Number>
42class Vector;
43
45
46template <typename Range = Vector<double>,
47 typename Domain = Range,
48 typename Payload =
49 internal::LinearOperatorImplementation::EmptyPayload>
50class LinearOperator;
51#endif
52
53template <
54 typename Range = Vector<double>,
55 typename Domain = Range,
57 typename OperatorExemplar,
58 typename Matrix>
60linear_operator(const OperatorExemplar &, const Matrix &);
61
62template <
63 typename Range = Vector<double>,
64 typename Domain = Range,
66 typename Matrix>
68linear_operator(const Matrix &);
69
70template <
71 typename Range = Vector<double>,
72 typename Domain = Range,
76
77template <typename Range, typename Domain, typename Payload>
80
81
178template <typename Range, typename Domain, typename Payload>
179class LinearOperator : public Payload
180{
181public:
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 <
242 typename Op,
243 typename = std::enable_if_t<
244 !std::is_base_of<LinearOperator<Range, Domain, Payload>, Op>::value>>
245 LinearOperator(const Op &op)
246 {
247 *this = linear_operator<Range, Domain, Payload, Op>(op);
248 }
249
255
260 template <
261 typename Op,
262 typename = std::enable_if_t<
263 !std::is_base_of<LinearOperator<Range, Domain, Payload>, Op>::value>>
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
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
388template <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
447template <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
485template <typename Range, typename Domain, typename Payload>
487operator*(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
552template <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
582template <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) {
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) {
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) {
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) {
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
658template <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
695template <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
753template <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
811template <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
831template <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
862template <
863 typename Range,
866identity_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
897template <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
917template <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
952template <
953 typename Range,
956mean_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
996template <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
1007namespace internal
1008{
1009 namespace LinearOperatorImplementation
1010 {
1022 template <typename Vector>
1024 {
1025 public:
1037 template <typename Matrix>
1038 static void
1039 reinit_range_vector(const Matrix &matrix,
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
1059 reinit_domain_vector(const Matrix &matrix,
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
1099 {
1100 return *this;
1101 }
1102
1103
1109 {
1110 return *this;
1111 }
1112
1113
1119 {
1120 return *this;
1121 }
1122
1123
1127 template <typename Solver, typename Preconditioner>
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
1378template <typename Range, typename Domain, typename Payload, typename Matrix>
1380linear_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
1402template <typename Range,
1403 typename Domain,
1404 typename Payload,
1405 typename OperatorExemplar,
1406 typename Matrix>
1408linear_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
1460template <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
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
1489template <
1490 typename Range = Vector<double>,
1491 typename Domain = Range,
1493 typename OperatorExemplar,
1494 typename Matrix,
1495 typename = std::enable_if_t<!std::is_lvalue_reference<Matrix>::value>>
1497linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1498
1499template <
1500 typename Range = Vector<double>,
1501 typename Domain = Range,
1503 typename OperatorExemplar,
1504 typename Matrix,
1505 typename =
1506 std::enable_if_t<!std::is_lvalue_reference<OperatorExemplar>::value>,
1507 typename = std::enable_if_t<
1508 !std::is_same<OperatorExemplar,
1511linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1512
1513template <
1514 typename Range = Vector<double>,
1515 typename Domain = Range,
1517 typename OperatorExemplar,
1518 typename Matrix,
1519 typename = std::enable_if_t<!std::is_lvalue_reference<Matrix>::value>,
1520 typename =
1521 std::enable_if_t<!std::is_lvalue_reference<OperatorExemplar>::value>,
1522 typename = std::enable_if_t<
1523 !std::is_same<OperatorExemplar,
1526linear_operator(OperatorExemplar &&, Matrix &&) = delete;
1527
1528template <
1529 typename Range = Vector<double>,
1530 typename Domain = Range,
1532 typename Matrix,
1533 typename = std::enable_if_t<!std::is_lvalue_reference<Matrix>::value>>
1536 Matrix &&) = delete;
1537
1538template <
1539 typename Range = Vector<double>,
1540 typename Domain = Range,
1542 typename Matrix,
1543 typename = std::enable_if_t<!std::is_lvalue_reference<Matrix>::value>>
1545linear_operator(Matrix &&) = delete;
1546
1547template <
1548 typename Payload,
1549 typename Solver,
1550 typename Preconditioner,
1551 typename Range = typename Solver::vector_type,
1552 typename Domain = Range,
1553 typename = std::enable_if_t<!std::is_lvalue_reference<Preconditioner>::value>,
1554 typename = std::enable_if_t<
1555 !std::is_same<Preconditioner, PreconditionIdentity>::value>,
1556 typename = std::enable_if_t<
1557 !std::is_same<Preconditioner,
1561 Solver &,
1562 Preconditioner &&) = delete;
1563
1564#endif // DOXYGEN
1565
1567
1568#endif
LinearOperator< Range, Domain, Payload > operator*=(typename Domain::value_type number)
LinearOperator< Range, Domain, Payload > & operator=(const LinearOperator< Range, Domain, Payload > &)=default
std::function< void(Range &v, const Domain &u)> vmult_add
LinearOperator< Range, Domain, Payload > & operator=(const Op &op)
std::function< void(Domain &v, const Range &u)> Tvmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
LinearOperator(const Payload &payload=Payload())
std::function< void(Range &v, const Domain &u)> vmult
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
LinearOperator< Range, Domain, Payload > & operator+=(const LinearOperator< Range, Domain, Payload > &second_op)
std::function< void(Domain &v, const Range &u)> Tvmult_add
LinearOperator(const Op &op)
LinearOperator< Range, Domain, Payload > & operator*=(const LinearOperator< Domain, Domain, Payload > &second_op)
LinearOperator(const LinearOperator< Range, Domain, Payload > &)=default
LinearOperator< Range, Domain, Payload > & operator-=(const LinearOperator< Range, Domain, Payload > &second_op)
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
EmptyPayload inverse_payload(Solver &, const Preconditioner &) const
void operator()(LinearOperator< Range, Domain, Payload > &op, const Matrix &matrix)
void operator()(LinearOperator< Range, Domain, Payload > &op, const Matrix &matrix)
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static auto test(Range *r, Domain *d) -> decltype(std::declval< C >().vmult_add(*r, *d), std::declval< C >().Tvmult_add(*d, *r), std::true_type())
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver)
LinearOperator< Range, Domain, Payload > operator*(const LinearOperator< Range, Intermediate, Payload > &first_op, const LinearOperator< Intermediate, Domain, Payload > &second_op)
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
LinearOperator< Range, Domain, Payload > operator-(const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
LinearOperator< Range, Range, Payload > identity_operator(const std::function< void(Range &, bool)> &reinit_vector)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > identity_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const PreconditionIdentity &)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const LinearOperator< Range, Domain, Payload > &preconditioner)
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &)
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > mean_value_filter(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > operator*(typename Range::value_type number, const LinearOperator< Range, Domain, Payload > &op)
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 > &op, typename Domain::value_type number)
LinearOperator< Range, Domain, Payload > linear_operator(const LinearOperator< Range, Domain, Payload > &operator_exemplar, const Matrix &matrix)
LinearOperator< Range, Range, Payload > mean_value_filter(const std::function< void(Range &, bool)> &reinit_vector)
LinearOperator< Range, Domain, Payload > identity_operator(const LinearOperator< Range, Domain, Payload > &)
LinearOperator< Range, Domain, Payload > operator+(const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
EmptyPayload operator+(const EmptyPayload &, const EmptyPayload &)
void apply_with_intermediate_storage(Function function, Range &v, const Domain &u, bool add)
EmptyPayload operator*(const EmptyPayload &, const EmptyPayload &)
static bool equal(const T *p1, const T *p2)