Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3067-g24f3489fdc 2025-04-14 21:10:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
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 if constexpr (running_in_debug_mode())
603 {
604 const unsigned int m = block_matrix.n_block_rows();
605 const unsigned int n = block_matrix.n_block_cols();
606 AssertIndexRange(i, m);
607 AssertIndexRange(j, n);
608 }
609
610 return BlockType(block_matrix.block(i, j));
611 };
612
613 populate_linear_operator_functions(return_op);
614 return return_op;
615}
616
617
618
646template <std::size_t m,
647 std::size_t n,
648 typename Range,
649 typename Domain,
650 typename BlockPayload>
653 const std::array<std::array<LinearOperator<typename Range::BlockType,
654 typename Domain::BlockType,
655 typename BlockPayload::BlockType>,
656 n>,
657 m> &ops)
658{
659 static_assert(m > 0 && n > 0,
660 "a blocked LinearOperator must consist of at least one block");
661
662 using BlockType =
664
665 // TODO: Create block payload so that this can be initialized correctly
666 BlockLinearOperator<Range, Domain, BlockPayload> return_op{BlockPayload()};
667
668 return_op.n_block_rows = []() -> unsigned int { return m; };
669
670 return_op.n_block_cols = []() -> unsigned int { return n; };
671
672 return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType {
673 AssertIndexRange(i, m);
674 AssertIndexRange(j, n);
675
676 return ops[i][j];
677 };
678
679 populate_linear_operator_functions(return_op);
680 return return_op;
681}
682
683
684
700template <typename Range = BlockVector<double>,
701 typename Domain = Range,
702 typename BlockPayload =
703 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
704 typename BlockMatrixType>
706block_diagonal_operator(const BlockMatrixType &block_matrix)
707{
708 using BlockType =
710
712 BlockPayload(block_matrix, block_matrix)};
713
714 return_op.n_block_rows = [&block_matrix]() -> unsigned int {
715 return block_matrix.n_block_rows();
716 };
717
718 return_op.n_block_cols = [&block_matrix]() -> unsigned int {
719 return block_matrix.n_block_cols();
720 };
721
722 return_op.block = [&block_matrix](unsigned int i,
723 unsigned int j) -> BlockType {
724 if constexpr (running_in_debug_mode())
725 {
726 const unsigned int m = block_matrix.n_block_rows();
727 const unsigned int n = block_matrix.n_block_cols();
728 Assert(m == n, ExcDimensionMismatch(m, n));
729 AssertIndexRange(i, m);
730 AssertIndexRange(j, n);
731 }
732 if (i == j)
733 return BlockType(block_matrix.block(i, j));
734 else
735 return null_operator(BlockType(block_matrix.block(i, j)));
736 };
737
738 populate_linear_operator_functions(return_op);
739 return return_op;
740}
741
742
743
761template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
764 const std::array<LinearOperator<typename Range::BlockType,
765 typename Domain::BlockType,
766 typename BlockPayload::BlockType>,
767 m> &ops)
768{
769 static_assert(
770 m > 0, "a blockdiagonal LinearOperator must consist of at least one block");
771
772 using BlockType =
774
775 std::array<std::array<BlockType, m>, m> new_ops;
776
777 // This is a bit tricky. We have to make sure that the off-diagonal
778 // elements of return_op.ops are populated correctly. They must be
779 // null_operators, but with correct reinit_domain_vector and
780 // reinit_range_vector functions.
781 for (unsigned int i = 0; i < m; ++i)
782 for (unsigned int j = 0; j < m; ++j)
783 if (i == j)
784 {
785 // diagonal elements are easy:
786 new_ops[i][j] = ops[i];
787 }
788 else
789 {
790 // create a null-operator...
791 new_ops[i][j] = null_operator(ops[i]);
792 // ... and fix up reinit_domain_vector:
793 new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
794 }
795
796 return block_operator<m, m, Range, Domain>(new_ops);
797}
798
799
800
810template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
813 const LinearOperator<typename Range::BlockType,
814 typename Domain::BlockType,
815 typename BlockPayload::BlockType> &op)
816{
817 static_assert(m > 0,
818 "a blockdiagonal LinearOperator must consist of at least "
819 "one block");
820
821 using BlockType =
823 std::array<BlockType, m> new_ops;
824 new_ops.fill(op);
825
826 return block_diagonal_operator(new_ops);
827}
828
829
830
871template <typename Range = BlockVector<double>,
872 typename Domain = Range,
873 typename BlockPayload =
874 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
879{
881 typename BlockPayload::BlockType(diagonal_inverse)};
882
883 return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
884 return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
885
886 return_op.vmult = [block_operator, diagonal_inverse](Range &v,
887 const Range &u) {
888 const unsigned int m = block_operator.n_block_rows();
889 Assert(block_operator.n_block_cols() == m,
890 ExcDimensionMismatch(block_operator.n_block_cols(), m));
891 Assert(diagonal_inverse.n_block_rows() == m,
892 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
893 Assert(diagonal_inverse.n_block_cols() == m,
894 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
895 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
896 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
897
898 if (m == 0)
899 return;
900
901 diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
902 for (unsigned int i = 1; i < m; ++i)
903 {
904 auto &dst = v.block(i);
905 dst = u.block(i);
906 dst *= -1.;
907 for (unsigned int j = 0; j < i; ++j)
908 block_operator.block(i, j).vmult_add(dst, v.block(j));
909 dst *= -1.;
910 diagonal_inverse.block(i, i).vmult(dst,
911 dst); // uses intermediate storage
912 }
913 };
914
915 return_op.vmult_add = [block_operator, diagonal_inverse](Range &v,
916 const Range &u) {
917 const unsigned int m = block_operator.n_block_rows();
918 Assert(block_operator.n_block_cols() == m,
919 ExcDimensionMismatch(block_operator.n_block_cols(), m));
920 Assert(diagonal_inverse.n_block_rows() == m,
921 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
922 Assert(diagonal_inverse.n_block_cols() == m,
923 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
924 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
925 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
926
927 if (m == 0)
928 return;
929
932 vector_memory);
933
934 diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
935
936 for (unsigned int i = 1; i < m; ++i)
937 {
938 diagonal_inverse.block(i, i).reinit_range_vector(
939 *tmp, /*bool omit_zeroing_entries=*/true);
940 *tmp = u.block(i);
941 *tmp *= -1.;
942 for (unsigned int j = 0; j < i; ++j)
943 block_operator.block(i, j).vmult_add(*tmp, v.block(j));
944 *tmp *= -1.;
945 diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
946 }
947 };
948
949 return return_op;
950}
951
952
953
988template <typename Range = BlockVector<double>,
989 typename Domain = Range,
990 typename BlockPayload =
991 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
996{
998 typename BlockPayload::BlockType(diagonal_inverse)};
999
1000 return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
1001 return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
1002
1003 return_op.vmult = [block_operator, diagonal_inverse](Range &v,
1004 const Range &u) {
1005 const unsigned int m = block_operator.n_block_rows();
1006 Assert(block_operator.n_block_cols() == m,
1007 ExcDimensionMismatch(block_operator.n_block_cols(), m));
1008 Assert(diagonal_inverse.n_block_rows() == m,
1009 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1010 Assert(diagonal_inverse.n_block_cols() == m,
1011 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1012 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1013 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1014
1015 if (m == 0)
1016 return;
1017
1018 diagonal_inverse.block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1019
1020 for (int i = m - 2; i >= 0; --i)
1021 {
1022 auto &dst = v.block(i);
1023 dst = u.block(i);
1024 dst *= -1.;
1025 for (unsigned int j = i + 1; j < m; ++j)
1026 block_operator.block(i, j).vmult_add(dst, v.block(j));
1027 dst *= -1.;
1028 diagonal_inverse.block(i, i).vmult(dst,
1029 dst); // uses intermediate storage
1030 }
1031 };
1032
1033 return_op.vmult_add = [block_operator, diagonal_inverse](Range &v,
1034 const Range &u) {
1035 const unsigned int m = block_operator.n_block_rows();
1036 Assert(block_operator.n_block_cols() == m,
1037 ExcDimensionMismatch(block_operator.n_block_cols(), m));
1038 Assert(diagonal_inverse.n_block_rows() == m,
1039 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1040 Assert(diagonal_inverse.n_block_cols() == m,
1041 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1042 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1043 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1046 vector_memory);
1047
1048 if (m == 0)
1049 return;
1050
1051 diagonal_inverse.block(m - 1, m - 1)
1052 .vmult_add(v.block(m - 1), u.block(m - 1));
1053
1054 for (int i = m - 2; i >= 0; --i)
1055 {
1056 diagonal_inverse.block(i, i).reinit_range_vector(
1057 *tmp, /*bool omit_zeroing_entries=*/true);
1058 *tmp = u.block(i);
1059 *tmp *= -1.;
1060 for (unsigned int j = i + 1; j < m; ++j)
1061 block_operator.block(i, j).vmult_add(*tmp, v.block(j));
1062 *tmp *= -1.;
1063 diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
1064 }
1065 };
1066
1067 return return_op;
1068}
1069
1073
1074#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:35
constexpr bool running_in_debug_mode()
Definition config.h:73
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#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)