Reference documentation for deal.II version 9.5.1
\(\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// Copyright (C) 2010 - 2022 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
16#ifndef dealii_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>
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
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
870template <typename Range = BlockVector<double>,
871 typename Domain = Range,
872 typename BlockPayload =
873 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
874LinearOperator<Domain, Range, typename BlockPayload::BlockType>
875block_forward_substitution(
876 const BlockLinearOperator<Range, Domain, BlockPayload> &block_operator,
877 const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
878{
879 LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op{
880 typename BlockPayload::BlockType(diagonal_inverse)};
881
882 return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
883 return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
884
885 return_op.vmult = [block_operator, diagonal_inverse](Range & v,
886 const Range &u) {
887 const unsigned int m = block_operator.n_block_rows();
888 Assert(block_operator.n_block_cols() == m,
889 ExcDimensionMismatch(block_operator.n_block_cols(), m));
890 Assert(diagonal_inverse.n_block_rows() == m,
891 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
892 Assert(diagonal_inverse.n_block_cols() == m,
893 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
894 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
895 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
896
897 if (m == 0)
898 return;
899
900 diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
901 for (unsigned int i = 1; i < m; ++i)
902 {
903 auto &dst = v.block(i);
904 dst = u.block(i);
905 dst *= -1.;
906 for (unsigned int j = 0; j < i; ++j)
907 block_operator.block(i, j).vmult_add(dst, v.block(j));
908 dst *= -1.;
909 diagonal_inverse.block(i, i).vmult(dst,
910 dst); // uses intermediate storage
911 }
912 };
913
914 return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
915 const Range &u) {
916 const unsigned int m = block_operator.n_block_rows();
917 Assert(block_operator.n_block_cols() == m,
918 ExcDimensionMismatch(block_operator.n_block_cols(), m));
919 Assert(diagonal_inverse.n_block_rows() == m,
920 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
921 Assert(diagonal_inverse.n_block_cols() == m,
922 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
923 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
924 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
925
926 if (m == 0)
927 return;
928
929 GrowingVectorMemory<typename Range::BlockType> vector_memory;
930 typename VectorMemory<typename Range::BlockType>::Pointer tmp(
931 vector_memory);
932
933 diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
934
935 for (unsigned int i = 1; i < m; ++i)
936 {
937 diagonal_inverse.block(i, i).reinit_range_vector(
938 *tmp, /*bool omit_zeroing_entries=*/true);
939 *tmp = u.block(i);
940 *tmp *= -1.;
941 for (unsigned int j = 0; j < i; ++j)
942 block_operator.block(i, j).vmult_add(*tmp, v.block(j));
943 *tmp *= -1.;
944 diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
945 }
946 };
947
948 return return_op;
949}
950
951
952
987template <typename Range = BlockVector<double>,
988 typename Domain = Range,
989 typename BlockPayload =
990 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
991LinearOperator<Domain, Range, typename BlockPayload::BlockType>
992block_back_substitution(
993 const BlockLinearOperator<Range, Domain, BlockPayload> &block_operator,
994 const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
995{
996 LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op{
997 typename BlockPayload::BlockType(diagonal_inverse)};
998
999 return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
1000 return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
1001
1002 return_op.vmult = [block_operator, diagonal_inverse](Range & v,
1003 const Range &u) {
1004 const unsigned int m = block_operator.n_block_rows();
1005 Assert(block_operator.n_block_cols() == m,
1006 ExcDimensionMismatch(block_operator.n_block_cols(), m));
1007 Assert(diagonal_inverse.n_block_rows() == m,
1008 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1009 Assert(diagonal_inverse.n_block_cols() == m,
1010 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1011 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1012 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1013
1014 if (m == 0)
1015 return;
1016
1017 diagonal_inverse.block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1018
1019 for (int i = m - 2; i >= 0; --i)
1020 {
1021 auto &dst = v.block(i);
1022 dst = u.block(i);
1023 dst *= -1.;
1024 for (unsigned int j = i + 1; j < m; ++j)
1025 block_operator.block(i, j).vmult_add(dst, v.block(j));
1026 dst *= -1.;
1027 diagonal_inverse.block(i, i).vmult(dst,
1028 dst); // uses intermediate storage
1029 }
1030 };
1031
1032 return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
1033 const Range &u) {
1034 const unsigned int m = block_operator.n_block_rows();
1035 Assert(block_operator.n_block_cols() == m,
1036 ExcDimensionMismatch(block_operator.n_block_cols(), m));
1037 Assert(diagonal_inverse.n_block_rows() == m,
1038 ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1039 Assert(diagonal_inverse.n_block_cols() == m,
1040 ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1041 Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1042 Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1043 GrowingVectorMemory<typename Range::BlockType> vector_memory;
1044 typename VectorMemory<typename Range::BlockType>::Pointer tmp(
1045 vector_memory);
1046
1047 if (m == 0)
1048 return;
1049
1050 diagonal_inverse.block(m - 1, m - 1)
1051 .vmult_add(v.block(m - 1), u.block(m - 1));
1052
1053 for (int i = m - 2; i >= 0; --i)
1054 {
1055 diagonal_inverse.block(i, i).reinit_range_vector(
1056 *tmp, /*bool omit_zeroing_entries=*/true);
1057 *tmp = u.block(i);
1058 *tmp *= -1.;
1059 for (unsigned int j = i + 1; j < m; ++j)
1060 block_operator.block(i, j).vmult_add(*tmp, v.block(j));
1061 *tmp *= -1.;
1062 diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
1063 }
1064 };
1065
1066 return return_op;
1067}
1068
1071DEAL_II_NAMESPACE_CLOSE
1072
1073#endif
ArrayView::reinit
void reinit(value_type *starting_element, const std::size_t n_elements)
Definition array_view.h:413
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 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
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exceptions.h:1855
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
block_diagonal_operator
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &)
Definition block_linear_operator.h:762
LinearOperator::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:875
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
null_operator
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &)
Definition linear_operator.h:919
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:992
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
block_operator
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &matrix)
Definition block_linear_operator.h:585
linear_operator.h
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:390
Generated by doxygen 1.9.8