Reference documentation for deal.II version Git e7bb9ce7b3 2020-09-18 12:07:32 -0400
\(\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 
22 
24 
25 
27 
28 // Forward declarations:
29 #ifndef DOXYGEN
30 namespace internal
31 {
32  namespace BlockLinearOperatorImplementation
33  {
34  template <typename PayloadBlockType =
36  class EmptyBlockPayload;
37  }
38 } // namespace internal
39 
40 template <typename Number>
41 class BlockVector;
42 
43 template <typename Range = BlockVector<double>,
44  typename Domain = Range,
45  typename BlockPayload =
46  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
48 #endif
49 
50 template <typename Range = BlockVector<double>,
51  typename Domain = Range,
52  typename BlockPayload =
53  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
54  typename BlockMatrixType>
56 block_operator(const BlockMatrixType &matrix);
57 
58 template <std::size_t m,
59  std::size_t n,
60  typename Range = BlockVector<double>,
61  typename Domain = Range,
62  typename BlockPayload =
66  const std::array<std::array<LinearOperator<typename Range::BlockType,
67  typename Domain::BlockType,
68  typename BlockPayload::BlockType>,
69  n>,
70  m> &);
71 
72 template <std::size_t m,
73  typename Range = BlockVector<double>,
74  typename Domain = Range,
75  typename BlockPayload =
79  const std::array<LinearOperator<typename Range::BlockType,
80  typename Domain::BlockType,
81  typename BlockPayload::BlockType>,
82  m> &);
83 
84 template <std::size_t m,
85  typename Range = BlockVector<double>,
86  typename Domain = Range,
87  typename BlockPayload =
91  const LinearOperator<typename Range::BlockType,
92  typename Domain::BlockType,
93  typename BlockPayload::BlockType> &op);
94 
95 
96 
165 template <typename Range, typename Domain, typename BlockPayload>
167  : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
168 {
169 public:
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 
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 
253  operator=(const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;
254 
259  template <typename Op>
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>
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>
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 
314 namespace 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
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>
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 
580 template <typename Range,
581  typename Domain,
582  typename BlockPayload,
583  typename BlockMatrixType>
585 block_operator(const BlockMatrixType &block_matrix)
586 {
587  using BlockType =
589 
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 
614  return return_op;
615 }
616 
617 
618 
646 template <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 
680  return return_op;
681 }
682 
683 
684 
700 template <typename Range = BlockVector<double>,
701  typename Domain = Range,
702  typename BlockPayload =
703  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
704  typename BlockMatrixType>
706 block_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 #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 
738  return return_op;
739 }
740 
741 
742 
760 template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
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 =
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 
809 template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
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 =
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 
871 template <typename Range = BlockVector<double>,
872  typename Domain = Range,
873  typename BlockPayload =
874  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
878  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
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 
989 template <typename Range = BlockVector<double>,
990  typename Domain = Range,
991  typename BlockPayload =
992  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
996  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
997 {
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));
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 
1074 
1075 #endif
BlockLinearOperator(const std::array< BlockType, m > &ops)
Contents is actually a matrix.
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
BlockLinearOperator(const Op &op)
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< BlockType, m > &ops)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &matrix)
void apply_with_intermediate_storage(const Function1 &first_op, const Function2 &loop_op, Range &v, const Domain &u, bool add)
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const Op &op)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &)
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< std::array< BlockType, n >, m > &ops)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1411
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
std::function< void(Range &v, const Domain &u)> vmult
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
void populate_linear_operator_functions(::BlockLinearOperator< Range, Domain, BlockPayload > &op)
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &)
BlockLinearOperator(const std::array< std::array< BlockType, n >, m > &ops)
std::function< void(Range &v, const Domain &u)> vmult_add
std::function< BlockType(unsigned int, unsigned int)> block
std::function< unsigned int()> n_block_cols
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
BlockLinearOperator(const BlockPayload &payload)
std::function< unsigned int()> n_block_rows
static bool equal(const T *p1, const T *p2)