Reference documentation for deal.II version GIT relicensing-1054-gedd46abbd6 2024-07-12 01:20:01+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
block_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_block_linear_operator_h
16#define dealii_block_linear_operator_h
17
18#include <deal.II/base/config.h>
19
21
23
24
26
27// Forward declarations:
28#ifndef DOXYGEN
29namespace internal
30{
31 namespace BlockLinearOperatorImplementation
32 {
33 template <typename PayloadBlockType =
35 class EmptyBlockPayload;
36 }
37} // namespace internal
38
39template <typename Number>
40class BlockVector;
41
42template <typename Range = BlockVector<double>,
43 typename Domain = Range,
44 typename BlockPayload =
45 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
47#endif
48
49template <typename Range = BlockVector<double>,
50 typename Domain = Range,
51 typename BlockPayload =
52 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
53 typename BlockMatrixType>
55block_operator(const BlockMatrixType &matrix);
56
57template <std::size_t m,
58 std::size_t n,
59 typename Range = BlockVector<double>,
60 typename Domain = Range,
61 typename BlockPayload =
65 const std::array<std::array<LinearOperator<typename Range::BlockType,
66 typename Domain::BlockType,
67 typename BlockPayload::BlockType>,
68 n>,
69 m> &);
70
71template <std::size_t m,
72 typename Range = BlockVector<double>,
73 typename Domain = Range,
74 typename BlockPayload =
78 const std::array<LinearOperator<typename Range::BlockType,
79 typename Domain::BlockType,
80 typename BlockPayload::BlockType>,
81 m> &);
82
83template <std::size_t m,
84 typename Range = BlockVector<double>,
85 typename Domain = Range,
86 typename BlockPayload =
90 const LinearOperator<typename Range::BlockType,
91 typename Domain::BlockType,
92 typename BlockPayload::BlockType> &op);
93
94
95
164template <typename Range, typename Domain, typename BlockPayload>
166 : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
167{
168public:
169 using BlockType = LinearOperator<typename Range::BlockType,
170 typename Domain::BlockType,
171 typename BlockPayload::BlockType>;
172
180 BlockLinearOperator(const BlockPayload &payload)
181 : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(
182 typename BlockPayload::BlockType(payload, payload))
183 {
184 n_block_rows = []() -> unsigned int {
185 Assert(
186 false,
188 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
189 return 0;
190 };
191
192 n_block_cols = []() -> unsigned int {
193 Assert(
194 false,
196 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
197 return 0;
198 };
199
200 block = [](unsigned int, unsigned int) -> BlockType {
201 Assert(
202 false,
204 "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
205 return BlockType();
206 };
207 }
208
214
220 template <typename Op>
222 {
223 *this = block_operator<Range, Domain, BlockPayload, Op>(op);
224 }
225
231 template <std::size_t m, std::size_t n>
232 BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
233 {
234 *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
235 }
236
242 template <std::size_t m>
243 BlockLinearOperator(const std::array<BlockType, m> &ops)
244 {
245 *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
246 }
247
253
258 template <typename Op>
260 operator=(const Op &op)
261 {
262 *this = block_operator<Range, Domain, BlockPayload, Op>(op);
263 return *this;
264 }
265
271 template <std::size_t m, std::size_t n>
273 operator=(const std::array<std::array<BlockType, n>, m> &ops)
274 {
275 *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
276 return *this;
277 }
278
284 template <std::size_t m>
286 operator=(const std::array<BlockType, m> &ops)
287 {
288 *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
289 return *this;
290 }
291
296 std::function<unsigned int()> n_block_rows;
297
302 std::function<unsigned int()> n_block_cols;
303
309 std::function<BlockType(unsigned int, unsigned int)> block;
310};
311
312
313namespace internal
314{
315 namespace BlockLinearOperatorImplementation
316 {
317 // A helper function to apply a given vmult, or Tvmult to a vector with
318 // intermediate storage, similar to the corresponding helper
319 // function for LinearOperator. Here, two operators are used.
320 // The first one takes care of the first "column" and typically doesn't add.
321 // On the other hand, the second operator is normally an adding one.
322 template <typename Function1,
323 typename Function2,
324 typename Range,
325 typename Domain>
326 void
327 apply_with_intermediate_storage(const Function1 &first_op,
328 const Function2 &loop_op,
329 Range &v,
330 const Domain &u,
331 bool add)
332 {
333 GrowingVectorMemory<Range> vector_memory;
334
335 typename VectorMemory<Range>::Pointer tmp(vector_memory);
336 tmp->reinit(v, /*bool omit_zeroing_entries =*/true);
337
338 const unsigned int n = u.n_blocks();
339 const unsigned int m = v.n_blocks();
340
341 for (unsigned int i = 0; i < m; ++i)
342 {
343 first_op(*tmp, u, i, 0);
344 for (unsigned int j = 1; j < n; ++j)
345 loop_op(*tmp, u, i, j);
346 }
347
348 if (add)
349 v += *tmp;
350 else
351 v = *tmp;
352 }
353
354 // Populate the LinearOperator interfaces with the help of the
355 // BlockLinearOperator functions
356 template <typename Range, typename Domain, typename BlockPayload>
357 inline void
360 {
361 op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries) {
362 const unsigned int m = op.n_block_rows();
363
364 // Reinitialize the block vector to m blocks:
365 v.reinit(m);
366
367 // And reinitialize every individual block with reinit_range_vectors:
368 for (unsigned int i = 0; i < m; ++i)
369 op.block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
370
371 v.collect_sizes();
372 };
373
374 op.reinit_domain_vector = [=](Domain &v, bool omit_zeroing_entries) {
375 const unsigned int n = op.n_block_cols();
376
377 // Reinitialize the block vector to n blocks:
378 v.reinit(n);
379
380 // And reinitialize every individual block with reinit_domain_vectors:
381 for (unsigned int i = 0; i < n; ++i)
382 op.block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
383
384 v.collect_sizes();
385 };
386
387 op.vmult = [&op](Range &v, const Domain &u) {
388 const unsigned int m = op.n_block_rows();
389 const unsigned int n = op.n_block_cols();
390 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
391 Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
392
393 if (PointerComparison::equal(&v, &u))
394 {
395 const auto first_op = [&op](Range &v,
396 const Domain &u,
397 const unsigned int i,
398 const unsigned int j) {
399 op.block(i, j).vmult(v.block(i), u.block(j));
400 };
401
402 const auto loop_op = [&op](Range &v,
403 const Domain &u,
404 const unsigned int i,
405 const unsigned int j) {
406 op.block(i, j).vmult_add(v.block(i), u.block(j));
407 };
408
409 apply_with_intermediate_storage(first_op, loop_op, v, u, false);
410 }
411 else
412 {
413 for (unsigned int i = 0; i < m; ++i)
414 {
415 op.block(i, 0).vmult(v.block(i), u.block(0));
416 for (unsigned int j = 1; j < n; ++j)
417 op.block(i, j).vmult_add(v.block(i), u.block(j));
418 }
419 }
420 };
421
422 op.vmult_add = [&op](Range &v, const Domain &u) {
423 const unsigned int m = op.n_block_rows();
424 const unsigned int n = op.n_block_cols();
425 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
426 Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
427
428 if (PointerComparison::equal(&v, &u))
429 {
430 const auto first_op = [&op](Range &v,
431 const Domain &u,
432 const unsigned int i,
433 const unsigned int j) {
434 op.block(i, j).vmult(v.block(i), u.block(j));
435 };
436
437 const auto loop_op = [&op](Range &v,
438 const Domain &u,
439 const unsigned int i,
440 const unsigned int j) {
441 op.block(i, j).vmult_add(v.block(i), u.block(j));
442 };
443
444 apply_with_intermediate_storage(first_op, loop_op, v, u, true);
445 }
446 else
447 {
448 for (unsigned int i = 0; i < m; ++i)
449 for (unsigned int j = 0; j < n; ++j)
450 op.block(i, j).vmult_add(v.block(i), u.block(j));
451 }
452 };
453
454 op.Tvmult = [&op](Domain &v, const Range &u) {
455 const unsigned int n = op.n_block_cols();
456 const unsigned int m = op.n_block_rows();
457 Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
458 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
459
460 if (PointerComparison::equal(&v, &u))
461 {
462 const auto first_op = [&op](Range &v,
463 const Domain &u,
464 const unsigned int i,
465 const unsigned int j) {
466 op.block(j, i).Tvmult(v.block(i), u.block(j));
467 };
468
469 const auto loop_op = [&op](Range &v,
470 const Domain &u,
471 const unsigned int i,
472 const unsigned int j) {
473 op.block(j, i).Tvmult_add(v.block(i), u.block(j));
474 };
475
476 apply_with_intermediate_storage(first_op, loop_op, v, u, false);
477 }
478 else
479 {
480 for (unsigned int i = 0; i < n; ++i)
481 {
482 op.block(0, i).Tvmult(v.block(i), u.block(0));
483 for (unsigned int j = 1; j < m; ++j)
484 op.block(j, i).Tvmult_add(v.block(i), u.block(j));
485 }
486 }
487 };
488
489 op.Tvmult_add = [&op](Domain &v, const Range &u) {
490 const unsigned int n = op.n_block_cols();
491 const unsigned int m = op.n_block_rows();
492 Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
493 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
494
495 if (PointerComparison::equal(&v, &u))
496 {
497 const auto first_op = [&op](Range &v,
498 const Domain &u,
499 const unsigned int i,
500 const unsigned int j) {
501 op.block(j, i).Tvmult(v.block(i), u.block(j));
502 };
503
504 const auto loop_op = [&op](Range &v,
505 const Domain &u,
506 const unsigned int i,
507 const unsigned int j) {
508 op.block(j, i).Tvmult_add(v.block(i), u.block(j));
509 };
510
511 apply_with_intermediate_storage(first_op, loop_op, v, u, true);
512 }
513 else
514 {
515 for (unsigned int i = 0; i < n; ++i)
516 for (unsigned int j = 0; j < m; ++j)
517 op.block(j, i).Tvmult_add(v.block(i), u.block(j));
518 }
519 };
520 }
521
522
523
537 template <typename PayloadBlockType>
539 {
540 public:
544 using BlockType = PayloadBlockType;
545
553 template <typename... Args>
554 EmptyBlockPayload(const Args &...)
555 {}
556 };
557
558 } // namespace BlockLinearOperatorImplementation
559} // namespace internal
560
561
562
579template <typename Range,
580 typename Domain,
581 typename BlockPayload,
582 typename BlockMatrixType>
584block_operator(const BlockMatrixType &block_matrix)
585{
586 using BlockType =
588
590 BlockPayload(block_matrix, block_matrix)};
591
592 return_op.n_block_rows = [&block_matrix]() -> unsigned int {
593 return block_matrix.n_block_rows();
594 };
595
596 return_op.n_block_cols = [&block_matrix]() -> unsigned int {
597 return block_matrix.n_block_cols();
598 };
599
600 return_op.block = [&block_matrix](unsigned int i,
601 unsigned int j) -> BlockType {
602#ifdef DEBUG
603 const unsigned int m = block_matrix.n_block_rows();
604 const unsigned int n = block_matrix.n_block_cols();
605 AssertIndexRange(i, m);
606 AssertIndexRange(j, n);
607#endif
608
609 return BlockType(block_matrix.block(i, j));
610 };
611
612 populate_linear_operator_functions(return_op);
613 return return_op;
614}
615
616
617
645template <std::size_t m,
646 std::size_t n,
647 typename Range,
648 typename Domain,
649 typename BlockPayload>
652 const std::array<std::array<LinearOperator<typename Range::BlockType,
653 typename Domain::BlockType,
654 typename BlockPayload::BlockType>,
655 n>,
656 m> &ops)
657{
658 static_assert(m > 0 && n > 0,
659 "a blocked LinearOperator must consist of at least one block");
660
661 using BlockType =
663
664 // TODO: Create block payload so that this can be initialized correctly
665 BlockLinearOperator<Range, Domain, BlockPayload> return_op{BlockPayload()};
666
667 return_op.n_block_rows = []() -> unsigned int { return m; };
668
669 return_op.n_block_cols = []() -> unsigned int { return n; };
670
671 return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType {
672 AssertIndexRange(i, m);
673 AssertIndexRange(j, n);
674
675 return ops[i][j];
676 };
677
678 populate_linear_operator_functions(return_op);
679 return return_op;
680}
681
682
683
699template <typename Range = BlockVector<double>,
700 typename Domain = Range,
701 typename BlockPayload =
702 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
703 typename BlockMatrixType>
705block_diagonal_operator(const BlockMatrixType &block_matrix)
706{
707 using BlockType =
709
711 BlockPayload(block_matrix, block_matrix)};
712
713 return_op.n_block_rows = [&block_matrix]() -> unsigned int {
714 return block_matrix.n_block_rows();
715 };
716
717 return_op.n_block_cols = [&block_matrix]() -> unsigned int {
718 return block_matrix.n_block_cols();
719 };
720
721 return_op.block = [&block_matrix](unsigned int i,
722 unsigned int j) -> BlockType {
723#ifdef DEBUG
724 const unsigned int m = block_matrix.n_block_rows();
725 const unsigned int n = block_matrix.n_block_cols();
726 Assert(m == n, ExcDimensionMismatch(m, n));
727 AssertIndexRange(i, m);
728 AssertIndexRange(j, n);
729#endif
730 if (i == j)
731 return BlockType(block_matrix.block(i, j));
732 else
733 return null_operator(BlockType(block_matrix.block(i, j)));
734 };
735
736 populate_linear_operator_functions(return_op);
737 return return_op;
738}
739
740
741
759template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
762 const std::array<LinearOperator<typename Range::BlockType,
763 typename Domain::BlockType,
764 typename BlockPayload::BlockType>,
765 m> &ops)
766{
767 static_assert(
768 m > 0, "a blockdiagonal LinearOperator must consist of at least one block");
769
770 using BlockType =
772
773 std::array<std::array<BlockType, m>, m> new_ops;
774
775 // This is a bit tricky. We have to make sure that the off-diagonal
776 // elements of return_op.ops are populated correctly. They must be
777 // null_operators, but with correct reinit_domain_vector and
778 // reinit_range_vector functions.
779 for (unsigned int i = 0; i < m; ++i)
780 for (unsigned int j = 0; j < m; ++j)
781 if (i == j)
782 {
783 // diagonal elements are easy:
784 new_ops[i][j] = ops[i];
785 }
786 else
787 {
788 // create a null-operator...
789 new_ops[i][j] = null_operator(ops[i]);
790 // ... and fix up reinit_domain_vector:
791 new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
792 }
793
794 return block_operator<m, m, Range, Domain>(new_ops);
795}
796
797
798
808template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
811 const LinearOperator<typename Range::BlockType,
812 typename Domain::BlockType,
813 typename BlockPayload::BlockType> &op)
814{
815 static_assert(m > 0,
816 "a blockdiagonal LinearOperator must consist of at least "
817 "one block");
818
819 using BlockType =
821 std::array<BlockType, m> new_ops;
822 new_ops.fill(op);
823
824 return block_diagonal_operator(new_ops);
825}
826
827
828
869template <typename Range = BlockVector<double>,
870 typename Domain = Range,
871 typename BlockPayload =
872 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
877{
879 typename BlockPayload::BlockType(diagonal_inverse)};
880
881 return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
882 return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
883
884 return_op.vmult = [block_operator, diagonal_inverse](Range &v,
885 const Range &u) {
886 const unsigned int m = block_operator.n_block_rows();
887 Assert(block_operator.n_block_cols() == m,
888 ExcDimensionMismatch(block_operator.n_block_cols(), m));
889 Assert(diagonal_inverse.n_block_rows() == m,
890 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
891 Assert(diagonal_inverse.n_block_cols() == m,
892 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
893 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
894 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
895
896 if (m == 0)
897 return;
898
899 diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
900 for (unsigned int i = 1; i < m; ++i)
901 {
902 auto &dst = v.block(i);
903 dst = u.block(i);
904 dst *= -1.;
905 for (unsigned int j = 0; j < i; ++j)
906 block_operator.block(i, j).vmult_add(dst, v.block(j));
907 dst *= -1.;
908 diagonal_inverse.block(i, i).vmult(dst,
909 dst); // uses intermediate storage
910 }
911 };
912
913 return_op.vmult_add = [block_operator, diagonal_inverse](Range &v,
914 const Range &u) {
915 const unsigned int m = block_operator.n_block_rows();
916 Assert(block_operator.n_block_cols() == m,
917 ExcDimensionMismatch(block_operator.n_block_cols(), m));
918 Assert(diagonal_inverse.n_block_rows() == m,
919 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
920 Assert(diagonal_inverse.n_block_cols() == m,
921 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
922 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
923 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
924
925 if (m == 0)
926 return;
927
930 vector_memory);
931
932 diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
933
934 for (unsigned int i = 1; i < m; ++i)
935 {
936 diagonal_inverse.block(i, i).reinit_range_vector(
937 *tmp, /*bool omit_zeroing_entries=*/true);
938 *tmp = u.block(i);
939 *tmp *= -1.;
940 for (unsigned int j = 0; j < i; ++j)
941 block_operator.block(i, j).vmult_add(*tmp, v.block(j));
942 *tmp *= -1.;
943 diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
944 }
945 };
946
947 return return_op;
948}
949
950
951
986template <typename Range = BlockVector<double>,
987 typename Domain = Range,
988 typename BlockPayload =
989 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
994{
996 typename BlockPayload::BlockType(diagonal_inverse)};
997
998 return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
999 return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
1000
1001 return_op.vmult = [block_operator, diagonal_inverse](Range &v,
1002 const Range &u) {
1003 const unsigned int m = block_operator.n_block_rows();
1004 Assert(block_operator.n_block_cols() == m,
1005 ExcDimensionMismatch(block_operator.n_block_cols(), m));
1006 Assert(diagonal_inverse.n_block_rows() == m,
1007 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1008 Assert(diagonal_inverse.n_block_cols() == m,
1009 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1010 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1011 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1012
1013 if (m == 0)
1014 return;
1015
1016 diagonal_inverse.block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1017
1018 for (int i = m - 2; i >= 0; --i)
1019 {
1020 auto &dst = v.block(i);
1021 dst = u.block(i);
1022 dst *= -1.;
1023 for (unsigned int j = i + 1; j < m; ++j)
1024 block_operator.block(i, j).vmult_add(dst, v.block(j));
1025 dst *= -1.;
1026 diagonal_inverse.block(i, i).vmult(dst,
1027 dst); // uses intermediate storage
1028 }
1029 };
1030
1031 return_op.vmult_add = [block_operator, diagonal_inverse](Range &v,
1032 const Range &u) {
1033 const unsigned int m = block_operator.n_block_rows();
1034 Assert(block_operator.n_block_cols() == m,
1035 ExcDimensionMismatch(block_operator.n_block_cols(), m));
1036 Assert(diagonal_inverse.n_block_rows() == m,
1037 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1038 Assert(diagonal_inverse.n_block_cols() == m,
1039 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1040 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1041 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1044 vector_memory);
1045
1046 if (m == 0)
1047 return;
1048
1049 diagonal_inverse.block(m - 1, m - 1)
1050 .vmult_add(v.block(m - 1), u.block(m - 1));
1051
1052 for (int i = m - 2; i >= 0; --i)
1053 {
1054 diagonal_inverse.block(i, i).reinit_range_vector(
1055 *tmp, /*bool omit_zeroing_entries=*/true);
1056 *tmp = u.block(i);
1057 *tmp *= -1.;
1058 for (unsigned int j = i + 1; j < m; ++j)
1059 block_operator.block(i, j).vmult_add(*tmp, v.block(j));
1060 *tmp *= -1.;
1061 diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
1062 }
1063 };
1064
1065 return return_op;
1066}
1067
1071
1072#endif
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< std::array< BlockType, n >, m > &ops)
BlockLinearOperator(const std::array< BlockType, m > &ops)
LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > BlockType
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< BlockType, m > &ops)
std::function< BlockType(unsigned int, unsigned int)> block
BlockLinearOperator(const BlockLinearOperator< Range, Domain, BlockPayload > &)=default
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const BlockLinearOperator< Range, Domain, BlockPayload > &)=default
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const Op &op)
std::function< unsigned int()> n_block_cols
std::function< unsigned int()> n_block_rows
BlockLinearOperator(const std::array< std::array< BlockType, n >, m > &ops)
BlockLinearOperator(const BlockPayload &payload)
std::function< void(Range &v, const Domain &u)> vmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
std::function< void(Range &v, const Domain &u)> vmult
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &)
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, n >, m > &ops)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &ops)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &)
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const BlockMatrixType &block_matrix)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > &op)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &block_matrix)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &matrix)
void populate_linear_operator_functions(::BlockLinearOperator< Range, Domain, BlockPayload > &op)
void apply_with_intermediate_storage(const Function1 &first_op, const Function2 &loop_op, Range &v, const Domain &u, bool add)
static bool equal(const T *p1, const T *p2)