Reference documentation for deal.II version GIT relicensing-218-g1f873f3929 2024-03-28 15:00:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2015 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_linear_operator_h
16#define dealii_linear_operator_h
17
18#include <deal.II/base/config.h>
19
21
23
24#include <array>
25#include <functional>
26#include <type_traits>
27
29
30// Forward declarations:
31#ifndef DOXYGEN
32namespace internal
33{
34 namespace LinearOperatorImplementation
35 {
36 class EmptyPayload;
37 }
38} // namespace internal
39
40template <typename Number>
41class Vector;
42
44
45template <typename Range = Vector<double>,
46 typename Domain = Range,
47 typename Payload =
48 internal::LinearOperatorImplementation::EmptyPayload>
49class LinearOperator;
50#endif
51
52template <
53 typename Range = Vector<double>,
54 typename Domain = Range,
56 typename OperatorExemplar,
57 typename Matrix>
59linear_operator(const OperatorExemplar &, const Matrix &);
60
61template <
62 typename Range = Vector<double>,
63 typename Domain = Range,
65 typename Matrix>
67linear_operator(const Matrix &);
68
69template <
70 typename Range = Vector<double>,
71 typename Domain = Range,
75
76template <typename Range, typename Domain, typename Payload>
79
80
199template <typename Range, typename Domain, typename Payload>
200class LinearOperator : public Payload
201{
202public:
211 LinearOperator(const Payload &payload = Payload())
212 : Payload(payload)
213 , is_null_operator(false)
214 {
215 vmult = [](Range &, const Domain &) {
216 Assert(false,
217 ExcMessage("Uninitialized LinearOperator<Range, "
218 "Domain>::vmult called"));
219 };
220
221 vmult_add = [](Range &, const Domain &) {
222 Assert(false,
223 ExcMessage("Uninitialized LinearOperator<Range, "
224 "Domain>::vmult_add called"));
225 };
226
227 Tvmult = [](Domain &, const Range &) {
228 Assert(false,
229 ExcMessage("Uninitialized LinearOperator<Range, "
230 "Domain>::Tvmult called"));
231 };
232
233 Tvmult_add = [](Domain &, const Range &) {
234 Assert(false,
235 ExcMessage("Uninitialized LinearOperator<Range, "
236 "Domain>::Tvmult_add called"));
237 };
238
239 reinit_range_vector = [](Range &, bool) {
240 Assert(false,
241 ExcMessage("Uninitialized LinearOperator<Range, "
242 "Domain>::reinit_range_vector method called"));
243 };
244
245 reinit_domain_vector = [](Domain &, bool) {
246 Assert(false,
247 ExcMessage("Uninitialized LinearOperator<Range, "
248 "Domain>::reinit_domain_vector method called"));
249 };
250 }
251
256
262 template <typename Op,
263 typename = std::enable_if_t<
264 !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
265 LinearOperator(const Op &op)
266 {
267 *this = linear_operator<Range, Domain, Payload, Op>(op);
268 }
269
275
280 template <typename Op,
281 typename = std::enable_if_t<
282 !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
284 operator=(const Op &op)
285 {
286 *this = linear_operator<Range, Domain, Payload, Op>(op);
287 return *this;
288 }
289
294 std::function<void(Range &v, const Domain &u)> vmult;
295
300 std::function<void(Range &v, const Domain &u)> vmult_add;
301
306 std::function<void(Domain &v, const Range &u)> Tvmult;
307
312 std::function<void(Domain &v, const Range &u)> Tvmult_add;
313
321 std::function<void(Range &v, bool omit_zeroing_entries)> reinit_range_vector;
322
330 std::function<void(Domain &v, bool omit_zeroing_entries)>
332
344 {
345 *this = *this + second_op;
346 return *this;
347 }
348
355 {
356 *this = *this - second_op;
357 return *this;
358 }
359
366 {
367 *this = *this * second_op;
368 return *this;
369 }
370
376 operator*=(typename Domain::value_type number)
377 {
378 *this = *this * number;
379 return *this;
380 }
381
388
390};
391
392
407template <typename Range, typename Domain, typename Payload>
411{
412 if (first_op.is_null_operator)
413 {
414 return second_op;
415 }
416 else if (second_op.is_null_operator)
417 {
418 return first_op;
419 }
420 else
421 {
423 static_cast<const Payload &>(first_op) +
424 static_cast<const Payload &>(second_op)};
425
426 return_op.reinit_range_vector = first_op.reinit_range_vector;
427 return_op.reinit_domain_vector = first_op.reinit_domain_vector;
428
429 // ensure to have valid computation objects by catching first_op and
430 // second_op by value
431
432 return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
433 first_op.vmult(v, u);
434 second_op.vmult_add(v, u);
435 };
436
437 return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
438 first_op.vmult_add(v, u);
439 second_op.vmult_add(v, u);
440 };
441
442 return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
443 second_op.Tvmult(v, u);
444 first_op.Tvmult_add(v, u);
445 };
446
447 return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
448 second_op.Tvmult_add(v, u);
449 first_op.Tvmult_add(v, u);
450 };
451
452 return return_op;
453 }
454}
455
456
466template <typename Range, typename Domain, typename Payload>
470{
471 if (first_op.is_null_operator)
472 {
473 return -1. * second_op;
474 }
475 else if (second_op.is_null_operator)
476 {
477 return first_op;
478 }
479 else
480 {
481 // implement with addition and scalar multiplication
482 return first_op + (-1. * second_op);
483 }
484}
485
486
504template <typename Range, typename Domain, typename Payload>
506operator*(typename Range::value_type number,
508{
509 static_assert(
510 std::is_convertible<typename Range::value_type,
511 typename Domain::value_type>::value,
512 "Range and Domain must have implicitly convertible 'value_type's");
513
514 if (op.is_null_operator)
515 {
516 return op;
517 }
518 else if (number == 0.)
519 {
520 return null_operator(op);
521 }
522 else
523 {
525
526 // ensure to have valid computation objects by catching number and op by
527 // value
528
529 return_op.vmult = [number, op](Range &v, const Domain &u) {
530 op.vmult(v, u);
531 v *= number;
532 };
533
534 return_op.vmult_add = [number, op](Range &v, const Domain &u) {
535 v /= number;
536 op.vmult_add(v, u);
537 v *= number;
538 };
539
540 return_op.Tvmult = [number, op](Domain &v, const Range &u) {
541 op.Tvmult(v, u);
542 v *= number;
543 };
544
545 return_op.Tvmult_add = [number, op](Domain &v, const Range &u) {
546 v /= number;
547 op.Tvmult_add(v, u);
548 v *= number;
549 };
550
551 return return_op;
552 }
553}
554
555
571template <typename Range, typename Domain, typename Payload>
574 typename Domain::value_type number)
575{
576 static_assert(
577 std::is_convertible<typename Range::value_type,
578 typename Domain::value_type>::value,
579 "Range and Domain must have implicitly convertible 'value_type's");
580
581 return number * op;
582}
583
601template <typename Range,
602 typename Intermediate,
603 typename Domain,
604 typename Payload>
608{
609 if (first_op.is_null_operator || second_op.is_null_operator)
610 {
612 return_op.reinit_domain_vector = second_op.reinit_domain_vector;
613 return_op.reinit_range_vector = first_op.reinit_range_vector;
614 return null_operator(return_op);
615 }
616 else
617 {
619 static_cast<const Payload &>(first_op) *
620 static_cast<const Payload &>(second_op)};
621
622 return_op.reinit_domain_vector = second_op.reinit_domain_vector;
623 return_op.reinit_range_vector = first_op.reinit_range_vector;
624
625 // ensure to have valid computation objects by catching first_op and
626 // second_op by value
627
628 return_op.vmult = [first_op, second_op](Range &v, const Domain &u) {
630
631 typename VectorMemory<Intermediate>::Pointer i(vector_memory);
632 second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
633 second_op.vmult(*i, u);
634 first_op.vmult(v, *i);
635 };
636
637 return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) {
639
640 typename VectorMemory<Intermediate>::Pointer i(vector_memory);
641 second_op.reinit_range_vector(*i, /*bool omit_zeroing_entries =*/true);
642 second_op.vmult(*i, u);
643 first_op.vmult_add(v, *i);
644 };
645
646 return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) {
648
649 typename VectorMemory<Intermediate>::Pointer i(vector_memory);
650 first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
651 first_op.Tvmult(*i, u);
652 second_op.Tvmult(v, *i);
653 };
654
655 return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) {
657
658 typename VectorMemory<Intermediate>::Pointer i(vector_memory);
659 first_op.reinit_domain_vector(*i, /*bool omit_zeroing_entries =*/true);
660 first_op.Tvmult(*i, u);
661 second_op.Tvmult_add(v, *i);
662 };
663
664 return return_op;
665 }
666}
667
668
677template <typename Range, typename Domain, typename Payload>
680{
681 LinearOperator<Domain, Range, Payload> return_op{op.transpose_payload()};
682
684 return_op.reinit_domain_vector = op.reinit_range_vector;
685
686 return_op.vmult = op.Tvmult;
687 return_op.vmult_add = op.Tvmult_add;
688 return_op.Tvmult = op.vmult;
689 return_op.Tvmult_add = op.vmult_add;
690
691 return return_op;
692}
693
694
714template <typename Payload,
715 typename Solver,
716 typename Preconditioner,
717 typename Range = typename Solver::vector_type,
718 typename Domain = Range>
721 Solver &solver,
722 const Preconditioner &preconditioner)
723{
725 op.inverse_payload(solver, preconditioner)};
726
728 return_op.reinit_domain_vector = op.reinit_range_vector;
729
730 return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
731 op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
732 solver.solve(op, v, u, preconditioner);
733 };
734
735 return_op.vmult_add = [op, &solver, &preconditioner](Range &v,
736 const Domain &u) {
737 GrowingVectorMemory<Range> vector_memory;
738
739 typename VectorMemory<Range>::Pointer v2(vector_memory);
740 op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
741 solver.solve(op, *v2, u, preconditioner);
742 v += *v2;
743 };
744
745 return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u) {
746 op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
747 solver.solve(transpose_operator(op), v, u, preconditioner);
748 };
749
750 return_op.Tvmult_add = [op, &solver, &preconditioner](Range &v,
751 const Domain &u) {
752 GrowingVectorMemory<Range> vector_memory;
753
754 typename VectorMemory<Range>::Pointer v2(vector_memory);
755 op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
756 solver.solve(transpose_operator(op), *v2, u, preconditioner);
757 v += *v2;
758 };
759
760 return return_op;
761}
762
763
772template <typename Payload,
773 typename Solver,
774 typename Range = typename Solver::vector_type,
775 typename Domain = Range>
778 Solver &solver,
779 const LinearOperator<Range, Domain, Payload> &preconditioner)
780{
782 op.inverse_payload(solver, preconditioner)};
783
785 return_op.reinit_domain_vector = op.reinit_range_vector;
786
787 return_op.vmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
788 op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
789 solver.solve(op, v, u, preconditioner);
790 };
791
792 return_op.vmult_add = [op, &solver, preconditioner](Range &v,
793 const Domain &u) {
794 GrowingVectorMemory<Range> vector_memory;
795
796 typename VectorMemory<Range>::Pointer v2(vector_memory);
797 op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
798 solver.solve(op, *v2, u, preconditioner);
799 v += *v2;
800 };
801
802 return_op.Tvmult = [op, &solver, preconditioner](Range &v, const Domain &u) {
803 op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/false);
804 solver.solve(transpose_operator(op), v, u, preconditioner);
805 };
806
807 return_op.Tvmult_add = [op, &solver, preconditioner](Range &v,
808 const Domain &u) {
809 GrowingVectorMemory<Range> vector_memory;
810
811 typename VectorMemory<Range>::Pointer v2(vector_memory);
812 op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
813 solver.solve(transpose_operator(op), *v2, u, preconditioner);
814 v += *v2;
815 };
816
817 return return_op;
818}
819
820
830template <typename Payload,
831 typename Solver,
832 typename Range = typename Solver::vector_type,
833 typename Domain = Range>
836 Solver &solver)
837{
838 return inverse_operator(op, solver, identity_operator(op));
839}
840
841
850template <typename Payload,
851 typename Solver,
852 typename Range = typename Solver::vector_type,
853 typename Domain = Range>
856 Solver &solver,
857 const PreconditionIdentity &)
858{
859 return inverse_operator(op, solver);
860}
861
881template <
882 typename Range,
885identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
886{
887 LinearOperator<Range, Range, Payload> return_op{Payload()};
888
889 return_op.reinit_range_vector = reinit_vector;
890 return_op.reinit_domain_vector = reinit_vector;
891
892 return_op.vmult = [](Range &v, const Range &u) { v = u; };
893
894 return_op.vmult_add = [](Range &v, const Range &u) { v += u; };
895
896 return_op.Tvmult = [](Range &v, const Range &u) { v = u; };
897
898 return_op.Tvmult_add = [](Range &v, const Range &u) { v += u; };
899
900 return return_op;
901}
902
903
916template <typename Range, typename Domain, typename Payload>
919{
920 auto return_op = identity_operator<Range, Payload>(op.reinit_range_vector);
921 static_cast<Payload &>(return_op) = op.identity_payload();
922
923 return return_op;
924}
925
926
936template <typename Range, typename Domain, typename Payload>
939{
940 LinearOperator<Range, Domain, Payload> return_op{op.null_payload()};
941
942 return_op.is_null_operator = true;
943
944 return_op.reinit_range_vector = op.reinit_range_vector;
945 return_op.reinit_domain_vector = op.reinit_domain_vector;
946
947 return_op.vmult = [](Range &v, const Domain &) { v = 0.; };
948
949 return_op.vmult_add = [](Range &, const Domain &) {};
950
951 return_op.Tvmult = [](Domain &v, const Range &) { v = 0.; };
952
953 return_op.Tvmult_add = [](Domain &, const Range &) {};
954
955 return return_op;
956}
957
958
971template <
972 typename Range,
975mean_value_filter(const std::function<void(Range &, bool)> &reinit_vector)
976{
977 LinearOperator<Range, Range, Payload> return_op{Payload()};
978
979 return_op.reinit_range_vector = reinit_vector;
980 return_op.reinit_domain_vector = reinit_vector;
981
982 return_op.vmult = [](Range &v, const Range &u) {
983 const auto mean = u.mean_value();
984
985 v = u;
986 v.add(-mean);
987 };
988
989 return_op.vmult_add = [](Range &v, const Range &u) {
990 const auto mean = u.mean_value();
991
992 v += u;
993 v.add(-mean);
994 };
995
996 return_op.Tvmult = return_op.vmult_add;
997 return_op.Tvmult_add = return_op.vmult_add;
998
999 return return_op;
1000}
1001
1002
1015template <typename Range, typename Domain, typename Payload>
1018{
1019 auto return_op = mean_value_filter<Range, Payload>(op.reinit_range_vector);
1020 static_cast<Payload &>(return_op) = op.identity_payload();
1021
1022 return return_op;
1023}
1024
1025
1026namespace internal
1027{
1028 namespace LinearOperatorImplementation
1029 {
1041 template <typename Vector>
1043 {
1044 public:
1056 template <typename Matrix>
1057 static void
1058 reinit_range_vector(const Matrix &matrix,
1059 Vector &v,
1060 bool omit_zeroing_entries)
1061 {
1062 v.reinit(matrix.m(), omit_zeroing_entries);
1063 }
1064
1076 template <typename Matrix>
1077 static void
1078 reinit_domain_vector(const Matrix &matrix,
1079 Vector &v,
1080 bool omit_zeroing_entries)
1081 {
1082 v.reinit(matrix.n(), omit_zeroing_entries);
1083 }
1084 };
1085
1086
1099 {
1100 public:
1108 template <typename... Args>
1109 EmptyPayload(const Args &...)
1110 {}
1111
1112
1118 {
1119 return *this;
1120 }
1121
1122
1128 {
1129 return *this;
1130 }
1131
1132
1138 {
1139 return *this;
1140 }
1141
1142
1146 template <typename Solver, typename Preconditioner>
1148 inverse_payload(Solver &, const Preconditioner &) const
1149 {
1150 return *this;
1151 }
1152 };
1153
1158 inline EmptyPayload
1160 {
1161 return {};
1162 }
1163
1168 inline EmptyPayload
1170 {
1171 return {};
1172 }
1173
1174
1175
1176 // A trait class that determines whether type T provides public
1177 // (templated or non-templated) vmult_add member functions
1178 template <typename Range, typename Domain, typename T>
1180 {
1181 template <typename C>
1182 static std::false_type
1183 test(...);
1184
1185 template <typename C>
1186 static auto
1187 test(Range *r, Domain *d)
1188 -> decltype(std::declval<C>().vmult_add(*r, *d),
1189 std::declval<C>().Tvmult_add(*d, *r),
1190 std::true_type());
1191
1192 public:
1193 // type is std::true_type if Matrix provides vmult_add and Tvmult_add,
1194 // otherwise it is std::false_type
1195
1196 using type = decltype(test<T>(nullptr, nullptr));
1197 };
1198
1199
1200 // A helper function to apply a given vmult, or Tvmult to a vector with
1201 // intermediate storage
1202 template <typename Function, typename Range, typename Domain>
1203 void
1205 Range &v,
1206 const Domain &u,
1207 bool add)
1208 {
1209 GrowingVectorMemory<Range> vector_memory;
1210
1211 typename VectorMemory<Range>::Pointer i(vector_memory);
1212 i->reinit(v, /*bool omit_zeroing_entries =*/true);
1213
1214 function(*i, u);
1215
1216 if (add)
1217 v += *i;
1218 else
1219 v = *i;
1220 }
1221
1222
1223 // A helper class to add a reduced matrix interface to a LinearOperator
1224 // (typically provided by Preconditioner classes)
1225 template <typename Range, typename Domain, typename Payload>
1227 {
1228 public:
1229 template <typename Matrix>
1230 void
1232 const Matrix &matrix)
1233 {
1234 op.vmult = [&matrix](Range &v, const Domain &u) {
1235 if (PointerComparison::equal(&v, &u))
1236 {
1237 // If v and u are the same memory location use intermediate
1238 // storage
1240 [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1241 v,
1242 u,
1243 /*bool add =*/false);
1244 }
1245 else
1246 {
1247 matrix.vmult(v, u);
1248 }
1249 };
1250
1251 op.vmult_add = [&matrix](Range &v, const Domain &u) {
1252 // use intermediate storage to implement vmult_add with vmult
1254 [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1255 v,
1256 u,
1257 /*bool add =*/true);
1258 };
1259
1260 op.Tvmult = [&matrix](Domain &v, const Range &u) {
1261 if (PointerComparison::equal(&v, &u))
1262 {
1263 // If v and u are the same memory location use intermediate
1264 // storage
1266 [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1267 v,
1268 u,
1269 /*bool add =*/false);
1270 }
1271 else
1272 {
1273 matrix.Tvmult(v, u);
1274 }
1275 };
1276
1277 op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1278 // use intermediate storage to implement Tvmult_add with Tvmult
1280 [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1281 v,
1282 u,
1283 /*bool add =*/true);
1284 };
1285 }
1286 };
1287
1288
1289 // A helper class to add the full matrix interface to a LinearOperator
1290 template <typename Range, typename Domain, typename Payload>
1292 {
1293 public:
1294 template <typename Matrix>
1295 void
1297 const Matrix &matrix)
1298 {
1299 // As above ...
1300
1302 op, matrix);
1303
1304 // ... but add native vmult_add and Tvmult_add variants:
1305
1306 op.vmult_add = [&matrix](Range &v, const Domain &u) {
1307 if (PointerComparison::equal(&v, &u))
1308 {
1310 [&matrix](Range &b, const Domain &a) { matrix.vmult(b, a); },
1311 v,
1312 u,
1313 /*bool add =*/true);
1314 }
1315 else
1316 {
1317 matrix.vmult_add(v, u);
1318 }
1319 };
1320
1321 op.Tvmult_add = [&matrix](Domain &v, const Range &u) {
1322 if (PointerComparison::equal(&v, &u))
1323 {
1325 [&matrix](Domain &b, const Range &a) { matrix.Tvmult(b, a); },
1326 v,
1327 u,
1328 /*bool add =*/true);
1329 }
1330 else
1331 {
1332 matrix.Tvmult_add(v, u);
1333 }
1334 };
1335 }
1336 };
1337 } // namespace LinearOperatorImplementation
1338} // namespace internal
1339
1340
1397template <typename Range, typename Domain, typename Payload, typename Matrix>
1399linear_operator(const Matrix &matrix)
1400{
1401 // implement with the more generic variant below...
1402 return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix,
1403 matrix);
1404}
1405
1406
1421template <typename Range,
1422 typename Domain,
1423 typename Payload,
1424 typename OperatorExemplar,
1425 typename Matrix>
1427linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
1428{
1430 // Initialize the payload based on the input exemplar matrix
1432 Payload(operator_exemplar, matrix)};
1433
1434 // Always store a reference to matrix and operator_exemplar in the lambda
1435 // functions. This ensures that a modification of the matrix after the
1436 // creation of a LinearOperator wrapper is respected - further a matrix
1437 // or an operator_exemplar cannot usually be copied...
1438
1439 return_op.reinit_range_vector =
1440 [&operator_exemplar](Range &v, bool omit_zeroing_entries) {
1442 Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1443 };
1444
1445 return_op.reinit_domain_vector = [&operator_exemplar](
1446 Domain &v, bool omit_zeroing_entries) {
1448 Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1449 };
1450
1451 std::conditional_t<
1455 .
1456 operator()(return_op, matrix);
1457
1458 return return_op;
1459}
1460
1461
1462
1479template <typename Range, typename Domain, typename Payload, typename Matrix>
1482 const Matrix &matrix)
1483{
1485 // Initialize the payload based on the LinearOperator exemplar
1486 auto return_op = operator_exemplar;
1487
1488 std::conditional_t<
1492 .
1493 operator()(return_op, matrix);
1494
1495 return return_op;
1496}
1497
1498
1501#ifndef DOXYGEN
1502
1503//
1504// Ensure that we never capture a reference to a temporary by accident.
1505// to avoid "stack use after free".
1506//
1507
1508template <
1509 typename Range = Vector<double>,
1510 typename Domain = Range,
1512 typename OperatorExemplar,
1513 typename Matrix,
1514 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1516linear_operator(const OperatorExemplar &, Matrix &&) = delete;
1517
1518template <
1519 typename Range = Vector<double>,
1520 typename Domain = Range,
1522 typename OperatorExemplar,
1523 typename Matrix,
1524 typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
1525 typename = std::enable_if_t<
1526 !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
1528linear_operator(OperatorExemplar &&, const Matrix &) = delete;
1529
1530template <
1531 typename Range = Vector<double>,
1532 typename Domain = Range,
1534 typename OperatorExemplar,
1535 typename Matrix,
1536 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>,
1537 typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
1538 typename = std::enable_if_t<
1539 !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
1541linear_operator(OperatorExemplar &&, Matrix &&) = delete;
1542
1543template <
1544 typename Range = Vector<double>,
1545 typename Domain = Range,
1547 typename Matrix,
1548 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1551 Matrix &&) = delete;
1552
1553template <
1554 typename Range = Vector<double>,
1555 typename Domain = Range,
1557 typename Matrix,
1558 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1560linear_operator(Matrix &&) = delete;
1561
1562template <
1563 typename Payload,
1564 typename Solver,
1565 typename Preconditioner,
1566 typename Range = typename Solver::vector_type,
1567 typename Domain = Range,
1568 typename = std::enable_if_t<!std::is_lvalue_reference_v<Preconditioner>>,
1569 typename =
1570 std::enable_if_t<!std::is_same_v<Preconditioner, PreconditionIdentity>>,
1571 typename = std::enable_if_t<
1572 !std::is_same_v<Preconditioner, LinearOperator<Range, Domain, Payload>>>>
1575 Solver &,
1576 Preconditioner &&) = delete;
1577
1578#endif // DOXYGEN
1579
1581
1582#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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#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)