Reference documentation for deal.II version 9.3.3
\(\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\}}\)
block_linear_operator.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2010 - 2020 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_block_linear_operator_h
17#define dealii_block_linear_operator_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/exceptions.h>
22
23#include <deal.II/lac/linear_operator.h>
24
25
26DEAL_II_NAMESPACE_OPEN
27
28// Forward declarations:
29#ifndef DOXYGEN
30namespace internal
31{
32 namespace BlockLinearOperatorImplementation
33 {
34 template <typename PayloadBlockType =
35 internal::LinearOperatorImplementation::EmptyPayload>
36 class EmptyBlockPayload;
37 }
38} // namespace internal
39
40template <typename Number>
41class BlockVector;
42
43template <typename Range = BlockVector<double>,
44 typename Domain = Range,
45 typename BlockPayload =
46 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
47class BlockLinearOperator;
48#endif
49
50template <typename Range = BlockVector<double>,
51 typename Domain = Range,
52 typename BlockPayload =
53 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
54 typename BlockMatrixType>
55BlockLinearOperator<Range, Domain, BlockPayload>
56block_operator(const BlockMatrixType &matrix);
57
58template <std::size_t m,
59 std::size_t n,
60 typename Range = BlockVector<double>,
61 typename Domain = Range,
62 typename BlockPayload =
63 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
64BlockLinearOperator<Range, Domain, BlockPayload>
65block_operator(
66 const std::array<std::array<LinearOperator<typename Range::BlockType,
67 typename Domain::BlockType,
68 typename BlockPayload::BlockType>,
69 n>,
70 m> &);
71
72template <std::size_t m,
73 typename Range = BlockVector<double>,
74 typename Domain = Range,
75 typename BlockPayload =
76 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
77BlockLinearOperator<Range, Domain, BlockPayload>
78block_diagonal_operator(
79 const std::array<LinearOperator<typename Range::BlockType,
80 typename Domain::BlockType,
81 typename BlockPayload::BlockType>,
82 m> &);
83
84template <std::size_t m,
85 typename Range = BlockVector<double>,
86 typename Domain = Range,
87 typename BlockPayload =
88 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
89BlockLinearOperator<Range, Domain, BlockPayload>
90block_diagonal_operator(
91 const LinearOperator<typename Range::BlockType,
92 typename Domain::BlockType,
93 typename BlockPayload::BlockType> &op);
94
95
96
165template <typename Range, typename Domain, typename BlockPayload>
166class BlockLinearOperator
167 : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
168{
169public:
170 using BlockType = LinearOperator<typename Range::BlockType,
171 typename Domain::BlockType,
172 typename BlockPayload::BlockType>;
173
181 BlockLinearOperator(const BlockPayload &payload)
182 : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(
183 typename BlockPayload::BlockType(payload, payload))
184 {
185 n_block_rows = []() -> unsigned int {
186 Assert(
187 false,
188 ExcMessage(
189 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
190 return 0;
191 };
192
193 n_block_cols = []() -> unsigned int {
194 Assert(
195 false,
196 ExcMessage(
197 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
198 return 0;
199 };
200
201 block = [](unsigned int, unsigned int) -> BlockType {
202 Assert(
203 false,
204 ExcMessage(
205 "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
206 return BlockType();
207 };
208 }
209
213 BlockLinearOperator(
214 const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;
215
221 template <typename Op>
222 BlockLinearOperator(const Op &op)
223 {
224 *this = block_operator<Range, Domain, BlockPayload, Op>(op);
225 }
226
232 template <std::size_t m, std::size_t n>
233 BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
234 {
235 *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
236 }
237
243 template <std::size_t m>
244 BlockLinearOperator(const std::array<BlockType, m> &ops)
245 {
246 *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
247 }
248
252 BlockLinearOperator<Range, Domain, BlockPayload> &
253 operator=(const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;
254
259 template <typename Op>
260 BlockLinearOperator<Range, Domain, BlockPayload> &
261 operator=(const Op &op)
262 {
263 *this = block_operator<Range, Domain, BlockPayload, Op>(op);
264 return *this;
265 }
266
272 template <std::size_t m, std::size_t n>
273 BlockLinearOperator<Range, Domain, BlockPayload> &
274 operator=(const std::array<std::array<BlockType, n>, m> &ops)
275 {
276 *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
277 return *this;
278 }
279
285 template <std::size_t m>
286 BlockLinearOperator<Range, Domain, BlockPayload> &
287 operator=(const std::array<BlockType, m> &ops)
288 {
289 *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
290 return *this;
291 }
292
297 std::function<unsigned int()> n_block_rows;
298
303 std::function<unsigned int()> n_block_cols;
304
310 std::function<BlockType(unsigned int, unsigned int)> block;
311};
312
313
314namespace internal
315{
316 namespace BlockLinearOperatorImplementation
317 {
318 // A helper function to apply a given vmult, or Tvmult to a vector with
319 // intermediate storage, similar to the corresponding helper
320 // function for LinearOperator. Here, two operators are used.
321 // The first one takes care of the first "column" and typically doesn't add.
322 // On the other hand, the second operator is normally an adding one.
323 template <typename Function1,
324 typename Function2,
325 typename Range,
326 typename Domain>
327 void
328 apply_with_intermediate_storage(const Function1 &first_op,
329 const Function2 &loop_op,
330 Range & v,
331 const Domain & u,
332 bool add)
333 {
334 GrowingVectorMemory<Range> vector_memory;
335
336 typename VectorMemory<Range>::Pointer tmp(vector_memory);
337 tmp->reinit(v, /*bool omit_zeroing_entries =*/true);
338
339 const unsigned int n = u.n_blocks();
340 const unsigned int m = v.n_blocks();
341
342 for (unsigned int i = 0; i < m; ++i)
343 {
344 first_op(*tmp, u, i, 0);
345 for (unsigned int j = 1; j < n; ++j)
346 loop_op(*tmp, u, i, j);
347 }
348
349 if (add)
350 v += *tmp;
351 else
352 v = *tmp;
353 }
354
355 // Populate the LinearOperator interfaces with the help of the
356 // BlockLinearOperator functions
357 template <typename Range, typename Domain, typename BlockPayload>
358 inline void
359 populate_linear_operator_functions(
360 ::BlockLinearOperator<Range, Domain, BlockPayload> &op)
361 {
362 op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries) {
363 const unsigned int m = op.n_block_rows();
364
365 // Reinitialize the block vector to m blocks:
366 v.reinit(m);
367
368 // And reinitialize every individual block with reinit_range_vectors:
369 for (unsigned int i = 0; i < m; ++i)
370 op.block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
371
372 v.collect_sizes();
373 };
374
375 op.reinit_domain_vector = [=](Domain &v, bool omit_zeroing_entries) {
376 const unsigned int n = op.n_block_cols();
377
378 // Reinitialize the block vector to n blocks:
379 v.reinit(n);
380
381 // And reinitialize every individual block with reinit_domain_vectors:
382 for (unsigned int i = 0; i < n; ++i)
383 op.block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
384
385 v.collect_sizes();
386 };
387
388 op.vmult = [&op](Range &v, const Domain &u) {
389 const unsigned int m = op.n_block_rows();
390 const unsigned int n = op.n_block_cols();
391 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
392 Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
393
394 if (PointerComparison::equal(&v, &u))
395 {
396 const auto first_op = [&op](Range & v,
397 const Domain & u,
398 const unsigned int i,
399 const unsigned int j) {
400 op.block(i, j).vmult(v.block(i), u.block(j));
401 };
402
403 const auto loop_op = [&op](Range & v,
404 const Domain & u,
405 const unsigned int i,
406 const unsigned int j) {
407 op.block(i, j).vmult_add(v.block(i), u.block(j));
408 };
409
410 apply_with_intermediate_storage(first_op, loop_op, v, u, false);
411 }
412 else
413 {
414 for (unsigned int i = 0; i < m; ++i)
415 {
416 op.block(i, 0).vmult(v.block(i), u.block(0));
417 for (unsigned int j = 1; j < n; ++j)
418 op.block(i, j).vmult_add(v.block(i), u.block(j));
419 }
420 }
421 };
422
423 op.vmult_add = [&op](Range &v, const Domain &u) {
424 const unsigned int m = op.n_block_rows();
425 const unsigned int n = op.n_block_cols();
426 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
427 Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
428
429 if (PointerComparison::equal(&v, &u))
430 {
431 const auto first_op = [&op](Range & v,
432 const Domain & u,
433 const unsigned int i,
434 const unsigned int j) {
435 op.block(i, j).vmult(v.block(i), u.block(j));
436 };
437
438 const auto loop_op = [&op](Range & v,
439 const Domain & u,
440 const unsigned int i,
441 const unsigned int j) {
442 op.block(i, j).vmult_add(v.block(i), u.block(j));
443 };
444
445 apply_with_intermediate_storage(first_op, loop_op, v, u, true);
446 }
447 else
448 {
449 for (unsigned int i = 0; i < m; ++i)
450 for (unsigned int j = 0; j < n; ++j)
451 op.block(i, j).vmult_add(v.block(i), u.block(j));
452 }
453 };
454
455 op.Tvmult = [&op](Domain &v, const Range &u) {
456 const unsigned int n = op.n_block_cols();
457 const unsigned int m = op.n_block_rows();
458 Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
459 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
460
461 if (PointerComparison::equal(&v, &u))
462 {
463 const auto first_op = [&op](Range & v,
464 const Domain & u,
465 const unsigned int i,
466 const unsigned int j) {
467 op.block(j, i).Tvmult(v.block(i), u.block(j));
468 };
469
470 const auto loop_op = [&op](Range & v,
471 const Domain & u,
472 const unsigned int i,
473 const unsigned int j) {
474 op.block(j, i).Tvmult_add(v.block(i), u.block(j));
475 };
476
477 apply_with_intermediate_storage(first_op, loop_op, v, u, false);
478 }
479 else
480 {
481 for (unsigned int i = 0; i < n; ++i)
482 {
483 op.block(0, i).Tvmult(v.block(i), u.block(0));
484 for (unsigned int j = 1; j < m; ++j)
485 op.block(j, i).Tvmult_add(v.block(i), u.block(j));
486 }
487 }
488 };
489
490 op.Tvmult_add = [&op](Domain &v, const Range &u) {
491 const unsigned int n = op.n_block_cols();
492 const unsigned int m = op.n_block_rows();
493 Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
494 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
495
496 if (PointerComparison::equal(&v, &u))
497 {
498 const auto first_op = [&op](Range & v,
499 const Domain & u,
500 const unsigned int i,
501 const unsigned int j) {
502 op.block(j, i).Tvmult(v.block(i), u.block(j));
503 };
504
505 const auto loop_op = [&op](Range & v,
506 const Domain & u,
507 const unsigned int i,
508 const unsigned int j) {
509 op.block(j, i).Tvmult_add(v.block(i), u.block(j));
510 };
511
512 apply_with_intermediate_storage(first_op, loop_op, v, u, true);
513 }
514 else
515 {
516 for (unsigned int i = 0; i < n; ++i)
517 for (unsigned int j = 0; j < m; ++j)
518 op.block(j, i).Tvmult_add(v.block(i), u.block(j));
519 }
520 };
521 }
522
523
524
538 template <typename PayloadBlockType>
539 class EmptyBlockPayload
540 {
541 public:
545 using BlockType = PayloadBlockType;
546
554 template <typename... Args>
555 EmptyBlockPayload(const Args &...)
556 {}
557 };
558
559 } // namespace BlockLinearOperatorImplementation
560} // namespace internal
561
562
563
568
580template <typename Range,
581 typename Domain,
582 typename BlockPayload,
583 typename BlockMatrixType>
584BlockLinearOperator<Range, Domain, BlockPayload>
585block_operator(const BlockMatrixType &block_matrix)
586{
587 using BlockType =
588 typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;
589
590 BlockLinearOperator<Range, Domain, BlockPayload> return_op{
591 BlockPayload(block_matrix, block_matrix)};
592
593 return_op.n_block_rows = [&block_matrix]() -> unsigned int {
594 return block_matrix.n_block_rows();
595 };
596
597 return_op.n_block_cols = [&block_matrix]() -> unsigned int {
598 return block_matrix.n_block_cols();
599 };
600
601 return_op.block = [&block_matrix](unsigned int i,
602 unsigned int j) -> BlockType {
603#ifdef DEBUG
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#endif
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>
651BlockLinearOperator<Range, Domain, BlockPayload>
652block_operator(
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 =
663 typename BlockLinearOperator<Range, Domain, BlockPayload>::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>
705BlockLinearOperator<Range, Domain, BlockPayload>
706block_diagonal_operator(const BlockMatrixType &block_matrix)
707{
708 using BlockType =
709 typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;
710
711 BlockLinearOperator<Range, Domain, BlockPayload> return_op{
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#ifdef DEBUG
725 const unsigned int m = block_matrix.n_block_rows();
726 const unsigned int n = block_matrix.n_block_cols();
727 Assert(m == n, ExcDimensionMismatch(m, n));
728 AssertIndexRange(i, m);
729 AssertIndexRange(j, n);
730#endif
731 if (i == j)
732 return BlockType(block_matrix.block(i, j));
733 else
734 return null_operator(BlockType(block_matrix.block(i, j)));
735 };
736
737 populate_linear_operator_functions(return_op);
738 return return_op;
739}
740
741
742
760template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
761BlockLinearOperator<Range, Domain, BlockPayload>
762block_diagonal_operator(
763 const std::array<LinearOperator<typename Range::BlockType,
764 typename Domain::BlockType,
765 typename BlockPayload::BlockType>,
766 m> &ops)
767{
768 static_assert(
769 m > 0, "a blockdiagonal LinearOperator must consist of at least one block");
770
771 using BlockType =
772 typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;
773
774 std::array<std::array<BlockType, m>, m> new_ops;
775
776 // This is a bit tricky. We have to make sure that the off-diagonal
777 // elements of return_op.ops are populated correctly. They must be
778 // null_operators, but with correct reinit_domain_vector and
779 // reinit_range_vector functions.
780 for (unsigned int i = 0; i < m; ++i)
781 for (unsigned int j = 0; j < m; ++j)
782 if (i == j)
783 {
784 // diagonal elements are easy:
785 new_ops[i][j] = ops[i];
786 }
787 else
788 {
789 // create a null-operator...
790 new_ops[i][j] = null_operator(ops[i]);
791 // ... and fix up reinit_domain_vector:
792 new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
793 }
794
795 return block_operator<m, m, Range, Domain>(new_ops);
796}
797
798
799
809template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
810BlockLinearOperator<Range, Domain, BlockPayload>
811block_diagonal_operator(
812 const LinearOperator<typename Range::BlockType,
813 typename Domain::BlockType,
814 typename BlockPayload::BlockType> &op)
815{
816 static_assert(m > 0,
817 "a blockdiagonal LinearOperator must consist of at least "
818 "one block");
819
820 using BlockType =
821 typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType;
822 std::array<BlockType, m> new_ops;
823 new_ops.fill(op);
824
825 return block_diagonal_operator(new_ops);
826}
827
828
829
831
835
871template <typename Range = BlockVector<double>,
872 typename Domain = Range,
873 typename BlockPayload =
874 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
875LinearOperator<Domain, Range, typename BlockPayload::BlockType>
876block_forward_substitution(
877 const BlockLinearOperator<Range, Domain, BlockPayload> &block_operator,
878 const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
879{
880 LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op{
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
930 GrowingVectorMemory<typename Range::BlockType> vector_memory;
931 typename VectorMemory<typename Range::BlockType>::Pointer tmp(
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
989template <typename Range = BlockVector<double>,
990 typename Domain = Range,
991 typename BlockPayload =
992 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
993LinearOperator<Domain, Range, typename BlockPayload::BlockType>
994block_back_substitution(
995 const BlockLinearOperator<Range, Domain, BlockPayload> &block_operator,
996 const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
997{
998 LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op{
999 typename BlockPayload::BlockType(diagonal_inverse)};
1000
1001 return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
1002 return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
1003
1004 return_op.vmult = [block_operator, diagonal_inverse](Range & v,
1005 const Range &u) {
1006 const unsigned int m = block_operator.n_block_rows();
1007 Assert(block_operator.n_block_cols() == m,
1008 ExcDimensionMismatch(block_operator.n_block_cols(), m));
1009 Assert(diagonal_inverse.n_block_rows() == m,
1010 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1011 Assert(diagonal_inverse.n_block_cols() == m,
1012 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1013 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1014 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1015
1016 if (m == 0)
1017 return;
1018
1019 diagonal_inverse.block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1020
1021 for (int i = m - 2; i >= 0; --i)
1022 {
1023 auto &dst = v.block(i);
1024 dst = u.block(i);
1025 dst *= -1.;
1026 for (unsigned int j = i + 1; j < m; ++j)
1027 block_operator.block(i, j).vmult_add(dst, v.block(j));
1028 dst *= -1.;
1029 diagonal_inverse.block(i, i).vmult(dst,
1030 dst); // uses intermediate storage
1031 }
1032 };
1033
1034 return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
1035 const Range &u) {
1036 const unsigned int m = block_operator.n_block_rows();
1037 Assert(block_operator.n_block_cols() == m,
1038 ExcDimensionMismatch(block_operator.n_block_cols(), m));
1039 Assert(diagonal_inverse.n_block_rows() == m,
1040 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1041 Assert(diagonal_inverse.n_block_cols() == m,
1042 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1043 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1044 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1045 GrowingVectorMemory<typename Range::BlockType> vector_memory;
1046 typename VectorMemory<typename Range::BlockType>::Pointer tmp(
1047 vector_memory);
1048
1049 if (m == 0)
1050 return;
1051
1052 diagonal_inverse.block(m - 1, m - 1)
1053 .vmult_add(v.block(m - 1), u.block(m - 1));
1054
1055 for (int i = m - 2; i >= 0; --i)
1056 {
1057 diagonal_inverse.block(i, i).reinit_range_vector(
1058 *tmp, /*bool omit_zeroing_entries=*/true);
1059 *tmp = u.block(i);
1060 *tmp *= -1.;
1061 for (unsigned int j = i + 1; j < m; ++j)
1062 block_operator.block(i, j).vmult_add(*tmp, v.block(j));
1063 *tmp *= -1.;
1064 diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
1065 }
1066 };
1067
1068 return return_op;
1069}
1070
1072
1073DEAL_II_NAMESPACE_CLOSE
1074
1075#endif
BlockLinearOperator::operator=
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< std::array< BlockType, n >, m > &ops)
Definition: block_linear_operator.h:274
BlockLinearOperator::BlockLinearOperator
BlockLinearOperator(const std::array< BlockType, m > &ops)
Definition: block_linear_operator.h:244
BlockLinearOperator::BlockType
LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > BlockType
Definition: block_linear_operator.h:172
BlockLinearOperator::operator=
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< BlockType, m > &ops)
Definition: block_linear_operator.h:287
BlockLinearOperator::block
std::function< BlockType(unsigned int, unsigned int)> block
Definition: block_linear_operator.h:310
BlockLinearOperator::BlockLinearOperator
BlockLinearOperator(const BlockLinearOperator< Range, Domain, BlockPayload > &)=default
BlockLinearOperator::operator=
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const BlockLinearOperator< Range, Domain, BlockPayload > &)=default
BlockLinearOperator::operator=
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const Op &op)
Definition: block_linear_operator.h:261
BlockLinearOperator::n_block_cols
std::function< unsigned int()> n_block_cols
Definition: block_linear_operator.h:303
BlockLinearOperator::n_block_rows
std::function< unsigned int()> n_block_rows
Definition: block_linear_operator.h:297
BlockLinearOperator::BlockLinearOperator
BlockLinearOperator(const Op &op)
Definition: block_linear_operator.h:222
BlockLinearOperator::BlockLinearOperator
BlockLinearOperator(const std::array< std::array< BlockType, n >, m > &ops)
Definition: block_linear_operator.h:233
BlockLinearOperator::BlockLinearOperator
BlockLinearOperator(const BlockPayload &payload)
Definition: block_linear_operator.h:181
LinearOperator::vmult_add
std::function< void(Range &v, const Domain &u)> vmult_add
Definition: linear_operator.h:281
LinearOperator::Tvmult
std::function< void(Domain &v, const Range &u)> Tvmult
Definition: linear_operator.h:287
LinearOperator::reinit_domain_vector
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
Definition: linear_operator.h:312
LinearOperator::vmult
std::function< void(Range &v, const Domain &u)> vmult
Definition: linear_operator.h:275
LinearOperator::reinit_range_vector
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
Definition: linear_operator.h:302
LinearOperator::Tvmult_add
std::function< void(Domain &v, const Range &u)> Tvmult_add
Definition: linear_operator.h:293
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
LinearOperator::null_operator
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &op)
Definition: linear_operator.h:919
BlockLinearOperator::block_forward_substitution
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
Definition: block_linear_operator.h:876
BlockLinearOperator::block_operator
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, n >, m > &ops)
Definition: block_linear_operator.h:652
BlockLinearOperator::block_diagonal_operator
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &ops)
Definition: block_linear_operator.h:762
BlockLinearOperator::block_back_substitution
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
Definition: block_linear_operator.h:994
BlockLinearOperator::block_diagonal_operator
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const BlockMatrixType &block_matrix)
Definition: block_linear_operator.h:706
BlockLinearOperator::block_diagonal_operator
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > &op)
Definition: block_linear_operator.h:811
BlockLinearOperator::block_operator
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &block_matrix)
Definition: block_linear_operator.h:585
linear_operator.h
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:58
internal::BlockLinearOperatorImplementation::populate_linear_operator_functions
void populate_linear_operator_functions(::BlockLinearOperator< Range, Domain, BlockPayload > &op)
Definition: block_linear_operator.h:359
internal::BlockLinearOperatorImplementation::apply_with_intermediate_storage
void apply_with_intermediate_storage(const Function1 &first_op, const Function2 &loop_op, Range &v, const Domain &u, bool add)
Definition: block_linear_operator.h:328
PointerComparison::equal
static bool equal(const T *p1, const T *p2)
Definition: template_constraints.h:228
Generated by doxygen 1.9.4