Reference documentation for deal.II version Git 92c1895bf8 2020-06-04 17:48:28 -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 
166 template <typename Range, typename Domain, typename BlockPayload>
168  : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
169 {
170 public:
171  using BlockType = LinearOperator<typename Range::BlockType,
172  typename Domain::BlockType,
173  typename BlockPayload::BlockType>;
174 
182  BlockLinearOperator(const BlockPayload &payload)
183  : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(
184  typename BlockPayload::BlockType(payload, payload))
185  {
186  n_block_rows = []() -> unsigned int {
187  Assert(
188  false,
189  ExcMessage(
190  "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
191  return 0;
192  };
193 
194  n_block_cols = []() -> unsigned int {
195  Assert(
196  false,
197  ExcMessage(
198  "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
199  return 0;
200  };
201 
202  block = [](unsigned int, unsigned int) -> BlockType {
203  Assert(
204  false,
205  ExcMessage(
206  "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
207  return BlockType();
208  };
209  }
210 
216 
222  template <typename Op>
223  BlockLinearOperator(const Op &op)
224  {
225  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
226  }
227 
233  template <std::size_t m, std::size_t n>
234  BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
235  {
236  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
237  }
238 
244  template <std::size_t m>
245  BlockLinearOperator(const std::array<BlockType, m> &ops)
246  {
247  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
248  }
249 
254  operator=(const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;
255 
260  template <typename Op>
262  operator=(const Op &op)
263  {
264  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
265  return *this;
266  }
267 
273  template <std::size_t m, std::size_t n>
275  operator=(const std::array<std::array<BlockType, n>, m> &ops)
276  {
277  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
278  return *this;
279  }
280 
286  template <std::size_t m>
288  operator=(const std::array<BlockType, m> &ops)
289  {
290  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
291  return *this;
292  }
293 
298  std::function<unsigned int()> n_block_rows;
299 
304  std::function<unsigned int()> n_block_cols;
305 
311  std::function<BlockType(unsigned int, unsigned int)> block;
312 };
313 
314 
315 namespace internal
316 {
317  namespace BlockLinearOperatorImplementation
318  {
319  // A helper function to apply a given vmult, or Tvmult to a vector with
320  // intermediate storage, similar to the corresponding helper
321  // function for LinearOperator. Here, two operators are used.
322  // The first one takes care of the first "column" and typically doesn't add.
323  // On the other hand, the second operator is normally an adding one.
324  template <typename Function1,
325  typename Function2,
326  typename Range,
327  typename Domain>
328  void
329  apply_with_intermediate_storage(const Function1 &first_op,
330  const Function2 &loop_op,
331  Range & v,
332  const Domain & u,
333  bool add)
334  {
335  GrowingVectorMemory<Range> vector_memory;
336 
337  typename VectorMemory<Range>::Pointer tmp(vector_memory);
338  tmp->reinit(v, /*bool omit_zeroing_entries =*/true);
339 
340  const unsigned int n = u.n_blocks();
341  const unsigned int m = v.n_blocks();
342 
343  for (unsigned int i = 0; i < m; ++i)
344  {
345  first_op(*tmp, u, i, 0);
346  for (unsigned int j = 1; j < n; ++j)
347  loop_op(*tmp, u, i, j);
348  }
349 
350  if (add)
351  v += *tmp;
352  else
353  v = *tmp;
354  }
355 
356  // Populate the LinearOperator interfaces with the help of the
357  // BlockLinearOperator functions
358  template <typename Range, typename Domain, typename BlockPayload>
359  inline void
362  {
363  op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries) {
364  const unsigned int m = op.n_block_rows();
365 
366  // Reinitialize the block vector to m blocks:
367  v.reinit(m);
368 
369  // And reinitialize every individual block with reinit_range_vectors:
370  for (unsigned int i = 0; i < m; ++i)
371  op.block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
372 
373  v.collect_sizes();
374  };
375 
376  op.reinit_domain_vector = [=](Domain &v, bool omit_zeroing_entries) {
377  const unsigned int n = op.n_block_cols();
378 
379  // Reinitialize the block vector to n blocks:
380  v.reinit(n);
381 
382  // And reinitialize every individual block with reinit_domain_vectors:
383  for (unsigned int i = 0; i < n; ++i)
384  op.block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
385 
386  v.collect_sizes();
387  };
388 
389  op.vmult = [&op](Range &v, const Domain &u) {
390  const unsigned int m = op.n_block_rows();
391  const unsigned int n = op.n_block_cols();
392  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
393  Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
394 
395  if (PointerComparison::equal(&v, &u))
396  {
397  const auto first_op = [&op](Range & v,
398  const Domain & u,
399  const unsigned int i,
400  const unsigned int j) {
401  op.block(i, j).vmult(v.block(i), u.block(j));
402  };
403 
404  const auto loop_op = [&op](Range & v,
405  const Domain & u,
406  const unsigned int i,
407  const unsigned int j) {
408  op.block(i, j).vmult_add(v.block(i), u.block(j));
409  };
410 
411  apply_with_intermediate_storage(first_op, loop_op, v, u, false);
412  }
413  else
414  {
415  for (unsigned int i = 0; i < m; ++i)
416  {
417  op.block(i, 0).vmult(v.block(i), u.block(0));
418  for (unsigned int j = 1; j < n; ++j)
419  op.block(i, j).vmult_add(v.block(i), u.block(j));
420  }
421  }
422  };
423 
424  op.vmult_add = [&op](Range &v, const Domain &u) {
425  const unsigned int m = op.n_block_rows();
426  const unsigned int n = op.n_block_cols();
427  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
428  Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
429 
430  if (PointerComparison::equal(&v, &u))
431  {
432  const auto first_op = [&op](Range & v,
433  const Domain & u,
434  const unsigned int i,
435  const unsigned int j) {
436  op.block(i, j).vmult(v.block(i), u.block(j));
437  };
438 
439  const auto loop_op = [&op](Range & v,
440  const Domain & u,
441  const unsigned int i,
442  const unsigned int j) {
443  op.block(i, j).vmult_add(v.block(i), u.block(j));
444  };
445 
446  apply_with_intermediate_storage(first_op, loop_op, v, u, true);
447  }
448  else
449  {
450  for (unsigned int i = 0; i < m; ++i)
451  for (unsigned int j = 0; j < n; ++j)
452  op.block(i, j).vmult_add(v.block(i), u.block(j));
453  }
454  };
455 
456  op.Tvmult = [&op](Domain &v, const Range &u) {
457  const unsigned int n = op.n_block_cols();
458  const unsigned int m = op.n_block_rows();
459  Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
460  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
461 
462  if (PointerComparison::equal(&v, &u))
463  {
464  const auto first_op = [&op](Range & v,
465  const Domain & u,
466  const unsigned int i,
467  const unsigned int j) {
468  op.block(j, i).Tvmult(v.block(i), u.block(j));
469  };
470 
471  const auto loop_op = [&op](Range & v,
472  const Domain & u,
473  const unsigned int i,
474  const unsigned int j) {
475  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
476  };
477 
478  apply_with_intermediate_storage(first_op, loop_op, v, u, false);
479  }
480  else
481  {
482  for (unsigned int i = 0; i < n; ++i)
483  {
484  op.block(0, i).Tvmult(v.block(i), u.block(0));
485  for (unsigned int j = 1; j < m; ++j)
486  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
487  }
488  }
489  };
490 
491  op.Tvmult_add = [&op](Domain &v, const Range &u) {
492  const unsigned int n = op.n_block_cols();
493  const unsigned int m = op.n_block_rows();
494  Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
495  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
496 
497  if (PointerComparison::equal(&v, &u))
498  {
499  const auto first_op = [&op](Range & v,
500  const Domain & u,
501  const unsigned int i,
502  const unsigned int j) {
503  op.block(j, i).Tvmult(v.block(i), u.block(j));
504  };
505 
506  const auto loop_op = [&op](Range & v,
507  const Domain & u,
508  const unsigned int i,
509  const unsigned int j) {
510  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
511  };
512 
513  apply_with_intermediate_storage(first_op, loop_op, v, u, true);
514  }
515  else
516  {
517  for (unsigned int i = 0; i < n; ++i)
518  for (unsigned int j = 0; j < m; ++j)
519  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
520  }
521  };
522  }
523 
524 
525 
540  template <typename PayloadBlockType>
542  {
543  public:
547  using BlockType = PayloadBlockType;
548 
556  template <typename... Args>
557  EmptyBlockPayload(const Args &...)
558  {}
559  };
560 
561  } // namespace BlockLinearOperatorImplementation
562 } // namespace internal
563 
564 
565 
570 
582 template <typename Range,
583  typename Domain,
584  typename BlockPayload,
585  typename BlockMatrixType>
587 block_operator(const BlockMatrixType &block_matrix)
588 {
589  using BlockType =
591 
593  BlockPayload(block_matrix, block_matrix)};
594 
595  return_op.n_block_rows = [&block_matrix]() -> unsigned int {
596  return block_matrix.n_block_rows();
597  };
598 
599  return_op.n_block_cols = [&block_matrix]() -> unsigned int {
600  return block_matrix.n_block_cols();
601  };
602 
603  return_op.block = [&block_matrix](unsigned int i,
604  unsigned int j) -> BlockType {
605 #ifdef DEBUG
606  const unsigned int m = block_matrix.n_block_rows();
607  const unsigned int n = block_matrix.n_block_cols();
608  AssertIndexRange(i, m);
609  AssertIndexRange(j, n);
610 #endif
611 
612  return BlockType(block_matrix.block(i, j));
613  };
614 
616  return return_op;
617 }
618 
619 
620 
648 template <std::size_t m,
649  std::size_t n,
650  typename Range,
651  typename Domain,
652  typename BlockPayload>
655  const std::array<std::array<LinearOperator<typename Range::BlockType,
656  typename Domain::BlockType,
657  typename BlockPayload::BlockType>,
658  n>,
659  m> &ops)
660 {
661  static_assert(m > 0 && n > 0,
662  "a blocked LinearOperator must consist of at least one block");
663 
664  using BlockType =
666 
667  // TODO: Create block payload so that this can be initialized correctly
668  BlockLinearOperator<Range, Domain, BlockPayload> return_op{BlockPayload()};
669 
670  return_op.n_block_rows = []() -> unsigned int { return m; };
671 
672  return_op.n_block_cols = []() -> unsigned int { return n; };
673 
674  return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType {
675  AssertIndexRange(i, m);
676  AssertIndexRange(j, n);
677 
678  return ops[i][j];
679  };
680 
682  return return_op;
683 }
684 
685 
686 
702 template <typename Range = BlockVector<double>,
703  typename Domain = Range,
704  typename BlockPayload =
705  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
706  typename BlockMatrixType>
708 block_diagonal_operator(const BlockMatrixType &block_matrix)
709 {
710  using BlockType =
712 
714  BlockPayload(block_matrix, block_matrix)};
715 
716  return_op.n_block_rows = [&block_matrix]() -> unsigned int {
717  return block_matrix.n_block_rows();
718  };
719 
720  return_op.n_block_cols = [&block_matrix]() -> unsigned int {
721  return block_matrix.n_block_cols();
722  };
723 
724  return_op.block = [&block_matrix](unsigned int i,
725  unsigned int j) -> BlockType {
726 #ifdef DEBUG
727  const unsigned int m = block_matrix.n_block_rows();
728  const unsigned int n = block_matrix.n_block_cols();
729  Assert(m == n, ExcDimensionMismatch(m, n));
730  AssertIndexRange(i, m);
731  AssertIndexRange(j, n);
732 #endif
733  if (i == j)
734  return BlockType(block_matrix.block(i, j));
735  else
736  return null_operator(BlockType(block_matrix.block(i, j)));
737  };
738 
740  return return_op;
741 }
742 
743 
744 
762 template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
765  const std::array<LinearOperator<typename Range::BlockType,
766  typename Domain::BlockType,
767  typename BlockPayload::BlockType>,
768  m> &ops)
769 {
770  static_assert(
771  m > 0, "a blockdiagonal LinearOperator must consist of at least one block");
772 
773  using BlockType =
775 
776  std::array<std::array<BlockType, m>, m> new_ops;
777 
778  // This is a bit tricky. We have to make sure that the off-diagonal
779  // elements of return_op.ops are populated correctly. They must be
780  // null_operators, but with correct reinit_domain_vector and
781  // reinit_range_vector functions.
782  for (unsigned int i = 0; i < m; ++i)
783  for (unsigned int j = 0; j < m; ++j)
784  if (i == j)
785  {
786  // diagonal elements are easy:
787  new_ops[i][j] = ops[i];
788  }
789  else
790  {
791  // create a null-operator...
792  new_ops[i][j] = null_operator(ops[i]);
793  // ... and fix up reinit_domain_vector:
794  new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
795  }
796 
797  return block_operator<m, m, Range, Domain>(new_ops);
798 }
799 
800 
801 
811 template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
814  const LinearOperator<typename Range::BlockType,
815  typename Domain::BlockType,
816  typename BlockPayload::BlockType> &op)
817 {
818  static_assert(m > 0,
819  "a blockdiagonal LinearOperator must consist of at least "
820  "one block");
821 
822  using BlockType =
824  std::array<BlockType, m> new_ops;
825  new_ops.fill(op);
826 
827  return block_diagonal_operator(new_ops);
828 }
829 
830 
831 
833 
837 
873 template <typename Range = BlockVector<double>,
874  typename Domain = Range,
875  typename BlockPayload =
876  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
880  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
881 {
883  typename BlockPayload::BlockType(diagonal_inverse)};
884 
885  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
886  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
887 
888  return_op.vmult = [block_operator, diagonal_inverse](Range & v,
889  const Range &u) {
890  const unsigned int m = block_operator.n_block_rows();
891  Assert(block_operator.n_block_cols() == m,
892  ExcDimensionMismatch(block_operator.n_block_cols(), m));
893  Assert(diagonal_inverse.n_block_rows() == m,
894  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
895  Assert(diagonal_inverse.n_block_cols() == m,
896  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
897  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
898  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
899 
900  if (m == 0)
901  return;
902 
903  diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
904  for (unsigned int i = 1; i < m; ++i)
905  {
906  auto &dst = v.block(i);
907  dst = u.block(i);
908  dst *= -1.;
909  for (unsigned int j = 0; j < i; ++j)
910  block_operator.block(i, j).vmult_add(dst, v.block(j));
911  dst *= -1.;
912  diagonal_inverse.block(i, i).vmult(dst,
913  dst); // uses intermediate storage
914  }
915  };
916 
917  return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
918  const Range &u) {
919  const unsigned int m = block_operator.n_block_rows();
920  Assert(block_operator.n_block_cols() == m,
921  ExcDimensionMismatch(block_operator.n_block_cols(), m));
922  Assert(diagonal_inverse.n_block_rows() == m,
923  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
924  Assert(diagonal_inverse.n_block_cols() == m,
925  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
926  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
927  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
928 
929  if (m == 0)
930  return;
931 
934  vector_memory);
935 
936  diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
937 
938  for (unsigned int i = 1; i < m; ++i)
939  {
940  diagonal_inverse.block(i, i).reinit_range_vector(
941  *tmp, /*bool omit_zeroing_entries=*/true);
942  *tmp = u.block(i);
943  *tmp *= -1.;
944  for (unsigned int j = 0; j < i; ++j)
945  block_operator.block(i, j).vmult_add(*tmp, v.block(j));
946  *tmp *= -1.;
947  diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
948  }
949  };
950 
951  return return_op;
952 }
953 
954 
955 
991 template <typename Range = BlockVector<double>,
992  typename Domain = Range,
993  typename BlockPayload =
994  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
998  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
999 {
1001  typename BlockPayload::BlockType(diagonal_inverse)};
1002 
1003  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
1004  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
1005 
1006  return_op.vmult = [block_operator, diagonal_inverse](Range & v,
1007  const Range &u) {
1008  const unsigned int m = block_operator.n_block_rows();
1009  Assert(block_operator.n_block_cols() == m,
1010  ExcDimensionMismatch(block_operator.n_block_cols(), m));
1011  Assert(diagonal_inverse.n_block_rows() == m,
1012  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1013  Assert(diagonal_inverse.n_block_cols() == m,
1014  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1015  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1016  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1017 
1018  if (m == 0)
1019  return;
1020 
1021  diagonal_inverse.block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1022 
1023  for (int i = m - 2; i >= 0; --i)
1024  {
1025  auto &dst = v.block(i);
1026  dst = u.block(i);
1027  dst *= -1.;
1028  for (unsigned int j = i + 1; j < m; ++j)
1029  block_operator.block(i, j).vmult_add(dst, v.block(j));
1030  dst *= -1.;
1031  diagonal_inverse.block(i, i).vmult(dst,
1032  dst); // uses intermediate storage
1033  }
1034  };
1035 
1036  return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
1037  const Range &u) {
1038  const unsigned int m = block_operator.n_block_rows();
1039  Assert(block_operator.n_block_cols() == m,
1040  ExcDimensionMismatch(block_operator.n_block_cols(), m));
1041  Assert(diagonal_inverse.n_block_rows() == m,
1042  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1043  Assert(diagonal_inverse.n_block_cols() == m,
1044  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1045  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1046  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1049  vector_memory);
1050 
1051  if (m == 0)
1052  return;
1053 
1054  diagonal_inverse.block(m - 1, m - 1)
1055  .vmult_add(v.block(m - 1), u.block(m - 1));
1056 
1057  for (int i = m - 2; i >= 0; --i)
1058  {
1059  diagonal_inverse.block(i, i).reinit_range_vector(
1060  *tmp, /*bool omit_zeroing_entries=*/true);
1061  *tmp = u.block(i);
1062  *tmp *= -1.;
1063  for (unsigned int j = i + 1; j < m; ++j)
1064  block_operator.block(i, j).vmult_add(*tmp, v.block(j));
1065  *tmp *= -1.;
1066  diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
1067  }
1068  };
1069 
1070  return return_op;
1071 }
1072 
1074 
1076 
1077 #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:1638
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:1408
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:361
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:360
BlockLinearOperator(const BlockPayload &payload)
std::function< unsigned int()> n_block_rows
static bool equal(const T *p1, const T *p2)