Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
block_linear_operator.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2019 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 
26 DEAL_II_NAMESPACE_OPEN
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 // This is a workaround for a bug in <=gcc-4.7 that does not like partial
96 // template default values in combination with local lambda expressions [1]
97 //
98 // [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624
99 //
100 // Forward declare functions with partial template defaults:
101 
102 template <typename Range = BlockVector<double>,
103  typename Domain = Range,
104  typename BlockPayload =
105  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
106  typename BlockMatrixType>
108 block_diagonal_operator(const BlockMatrixType &block_matrix);
109 
110 template <typename Range = BlockVector<double>,
111  typename Domain = Range,
112  typename BlockPayload =
113  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
118 
119 template <typename Range = BlockVector<double>,
120  typename Domain = Range,
121  typename BlockPayload =
122  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
127 
128 // end of workaround
129 
130 
131 
201 template <typename Range, typename Domain, typename BlockPayload>
203  : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
204 {
205 public:
206  using BlockType = LinearOperator<typename Range::BlockType,
207  typename Domain::BlockType,
208  typename BlockPayload::BlockType>;
209 
217  BlockLinearOperator(const BlockPayload &payload)
218  : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(
219  typename BlockPayload::BlockType(payload, payload))
220  {
221  n_block_rows = []() -> unsigned int {
222  Assert(
223  false,
224  ExcMessage(
225  "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
226  return 0;
227  };
228 
229  n_block_cols = []() -> unsigned int {
230  Assert(
231  false,
232  ExcMessage(
233  "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
234  return 0;
235  };
236 
237  block = [](unsigned int, unsigned int) -> BlockType {
238  Assert(
239  false,
240  ExcMessage(
241  "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
242  return BlockType();
243  };
244  }
245 
251 
257  template <typename Op>
258  BlockLinearOperator(const Op &op)
259  {
260  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
261  }
262 
268  template <std::size_t m, std::size_t n>
269  BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
270  {
271  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
272  }
273 
279  template <std::size_t m>
280  BlockLinearOperator(const std::array<BlockType, m> &ops)
281  {
282  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
283  }
284 
289  operator=(const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;
290 
295  template <typename Op>
297  operator=(const Op &op)
298  {
299  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
300  return *this;
301  }
302 
308  template <std::size_t m, std::size_t n>
310  operator=(const std::array<std::array<BlockType, n>, m> &ops)
311  {
312  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
313  return *this;
314  }
315 
321  template <std::size_t m>
323  operator=(const std::array<BlockType, m> &ops)
324  {
325  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
326  return *this;
327  }
328 
333  std::function<unsigned int()> n_block_rows;
334 
339  std::function<unsigned int()> n_block_cols;
340 
346  std::function<BlockType(unsigned int, unsigned int)> block;
347 };
348 
349 
350 namespace internal
351 {
352  namespace BlockLinearOperatorImplementation
353  {
354  // A helper function to apply a given vmult, or Tvmult to a vector with
355  // intermediate storage, similar to the corresponding helper
356  // function for LinearOperator. Here, two operators are used.
357  // The first one takes care of the first "column" and typically doesn't add.
358  // On the other hand, the second operator is normally an adding one.
359  template <typename Function1,
360  typename Function2,
361  typename Range,
362  typename Domain>
363  void
364  apply_with_intermediate_storage(const Function1 &first_op,
365  const Function2 &loop_op,
366  Range & v,
367  const Domain & u,
368  bool add)
369  {
370  GrowingVectorMemory<Range> vector_memory;
371 
372  typename VectorMemory<Range>::Pointer tmp(vector_memory);
373  tmp->reinit(v, /*bool omit_zeroing_entries =*/true);
374 
375  const unsigned int n = u.n_blocks();
376  const unsigned int m = v.n_blocks();
377 
378  for (unsigned int i = 0; i < m; ++i)
379  {
380  first_op(*tmp, u, i, 0);
381  for (unsigned int j = 1; j < n; ++j)
382  loop_op(*tmp, u, i, j);
383  }
384 
385  if (add)
386  v += *tmp;
387  else
388  v = *tmp;
389  }
390 
391  // Populate the LinearOperator interfaces with the help of the
392  // BlockLinearOperator functions
393  template <typename Range, typename Domain, typename BlockPayload>
394  inline void
395  populate_linear_operator_functions(
397  {
398  op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries) {
399  const unsigned int m = op.n_block_rows();
400 
401  // Reinitialize the block vector to m blocks:
402  v.reinit(m);
403 
404  // And reinitialize every individual block with reinit_range_vectors:
405  for (unsigned int i = 0; i < m; ++i)
406  op.block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
407 
408  v.collect_sizes();
409  };
410 
411  op.reinit_domain_vector = [=](Domain &v, bool omit_zeroing_entries) {
412  const unsigned int n = op.n_block_cols();
413 
414  // Reinitialize the block vector to n blocks:
415  v.reinit(n);
416 
417  // And reinitialize every individual block with reinit_domain_vectors:
418  for (unsigned int i = 0; i < n; ++i)
419  op.block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
420 
421  v.collect_sizes();
422  };
423 
424  op.vmult = [&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, false);
447  }
448  else
449  {
450  for (unsigned int i = 0; i < m; ++i)
451  {
452  op.block(i, 0).vmult(v.block(i), u.block(0));
453  for (unsigned int j = 1; j < n; ++j)
454  op.block(i, j).vmult_add(v.block(i), u.block(j));
455  }
456  }
457  };
458 
459  op.vmult_add = [&op](Range &v, const Domain &u) {
460  const unsigned int m = op.n_block_rows();
461  const unsigned int n = op.n_block_cols();
462  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
463  Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
464 
465  if (PointerComparison::equal(&v, &u))
466  {
467  const auto first_op = [&op](Range & v,
468  const Domain & u,
469  const unsigned int i,
470  const unsigned int j) {
471  op.block(i, j).vmult(v.block(i), u.block(j));
472  };
473 
474  const auto loop_op = [&op](Range & v,
475  const Domain & u,
476  const unsigned int i,
477  const unsigned int j) {
478  op.block(i, j).vmult_add(v.block(i), u.block(j));
479  };
480 
481  apply_with_intermediate_storage(first_op, loop_op, v, u, true);
482  }
483  else
484  {
485  for (unsigned int i = 0; i < m; ++i)
486  for (unsigned int j = 0; j < n; ++j)
487  op.block(i, j).vmult_add(v.block(i), u.block(j));
488  }
489  };
490 
491  op.Tvmult = [&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, false);
514  }
515  else
516  {
517  for (unsigned int i = 0; i < n; ++i)
518  {
519  op.block(0, i).Tvmult(v.block(i), u.block(0));
520  for (unsigned int j = 1; j < m; ++j)
521  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
522  }
523  }
524  };
525 
526  op.Tvmult_add = [&op](Domain &v, const Range &u) {
527  const unsigned int n = op.n_block_cols();
528  const unsigned int m = op.n_block_rows();
529  Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
530  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
531 
532  if (PointerComparison::equal(&v, &u))
533  {
534  const auto first_op = [&op](Range & v,
535  const Domain & u,
536  const unsigned int i,
537  const unsigned int j) {
538  op.block(j, i).Tvmult(v.block(i), u.block(j));
539  };
540 
541  const auto loop_op = [&op](Range & v,
542  const Domain & u,
543  const unsigned int i,
544  const unsigned int j) {
545  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
546  };
547 
548  apply_with_intermediate_storage(first_op, loop_op, v, u, true);
549  }
550  else
551  {
552  for (unsigned int i = 0; i < n; ++i)
553  for (unsigned int j = 0; j < m; ++j)
554  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
555  }
556  };
557  }
558 
559 
560 
575  template <typename PayloadBlockType>
577  {
578  public:
582  using BlockType = PayloadBlockType;
583 
591  template <typename... Args>
592  EmptyBlockPayload(const Args &...)
593  {}
594  };
595 
596  } // namespace BlockLinearOperatorImplementation
597 } // namespace internal
598 
599 
600 
605 
617 template <typename Range,
618  typename Domain,
619  typename BlockPayload,
620  typename BlockMatrixType>
622 block_operator(const BlockMatrixType &block_matrix)
623 {
624  using BlockType =
626 
628  BlockPayload(block_matrix, block_matrix)};
629 
630  return_op.n_block_rows = [&block_matrix]() -> unsigned int {
631  return block_matrix.n_block_rows();
632  };
633 
634  return_op.n_block_cols = [&block_matrix]() -> unsigned int {
635  return block_matrix.n_block_cols();
636  };
637 
638  return_op.block = [&block_matrix](unsigned int i,
639  unsigned int j) -> BlockType {
640 #ifdef DEBUG
641  const unsigned int m = block_matrix.n_block_rows();
642  const unsigned int n = block_matrix.n_block_cols();
643  Assert(i < m, ExcIndexRange(i, 0, m));
644  Assert(j < n, ExcIndexRange(j, 0, n));
645 #endif
646 
647  return BlockType(block_matrix.block(i, j));
648  };
649 
650  populate_linear_operator_functions(return_op);
651  return return_op;
652 }
653 
654 
655 
683 template <std::size_t m,
684  std::size_t n,
685  typename Range,
686  typename Domain,
687  typename BlockPayload>
690  const std::array<std::array<LinearOperator<typename Range::BlockType,
691  typename Domain::BlockType,
692  typename BlockPayload::BlockType>,
693  n>,
694  m> &ops)
695 {
696  static_assert(m > 0 && n > 0,
697  "a blocked LinearOperator must consist of at least one block");
698 
699  using BlockType =
701 
702  // TODO: Create block payload so that this can be initialized correctly
703  BlockLinearOperator<Range, Domain, BlockPayload> return_op{BlockPayload()};
704 
705  return_op.n_block_rows = []() -> unsigned int { return m; };
706 
707  return_op.n_block_cols = []() -> unsigned int { return n; };
708 
709  return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType {
710  Assert(i < m, ExcIndexRange(i, 0, m));
711  Assert(j < n, ExcIndexRange(j, 0, n));
712 
713  return ops[i][j];
714  };
715 
716  populate_linear_operator_functions(return_op);
717  return return_op;
718 }
719 
720 
721 
737 template <typename Range,
738  typename Domain,
739  typename BlockPayload,
740  typename BlockMatrixType>
742 block_diagonal_operator(const BlockMatrixType &block_matrix)
743 {
744  using BlockType =
746 
748  BlockPayload(block_matrix, block_matrix)};
749 
750  return_op.n_block_rows = [&block_matrix]() -> unsigned int {
751  return block_matrix.n_block_rows();
752  };
753 
754  return_op.n_block_cols = [&block_matrix]() -> unsigned int {
755  return block_matrix.n_block_cols();
756  };
757 
758  return_op.block = [&block_matrix](unsigned int i,
759  unsigned int j) -> BlockType {
760 #ifdef DEBUG
761  const unsigned int m = block_matrix.n_block_rows();
762  const unsigned int n = block_matrix.n_block_cols();
763  Assert(m == n, ExcDimensionMismatch(m, n));
764  Assert(i < m, ExcIndexRange(i, 0, m));
765  Assert(j < n, ExcIndexRange(j, 0, n));
766 #endif
767  if (i == j)
768  return BlockType(block_matrix.block(i, j));
769  else
770  return null_operator(BlockType(block_matrix.block(i, j)));
771  };
772 
773  populate_linear_operator_functions(return_op);
774  return return_op;
775 }
776 
777 
778 
796 template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
799  const std::array<LinearOperator<typename Range::BlockType,
800  typename Domain::BlockType,
801  typename BlockPayload::BlockType>,
802  m> &ops)
803 {
804  static_assert(
805  m > 0, "a blockdiagonal LinearOperator must consist of at least one block");
806 
807  using BlockType =
809 
810  std::array<std::array<BlockType, m>, m> new_ops;
811 
812  // This is a bit tricky. We have to make sure that the off-diagonal
813  // elements of return_op.ops are populated correctly. They must be
814  // null_operators, but with correct reinit_domain_vector and
815  // reinit_range_vector functions.
816  for (unsigned int i = 0; i < m; ++i)
817  for (unsigned int j = 0; j < m; ++j)
818  if (i == j)
819  {
820  // diagonal elements are easy:
821  new_ops[i][j] = ops[i];
822  }
823  else
824  {
825  // create a null-operator...
826  new_ops[i][j] = null_operator(ops[i]);
827  // ... and fix up reinit_domain_vector:
828  new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
829  }
830 
831  return block_operator<m, m, Range, Domain>(new_ops);
832 }
833 
834 
835 
845 template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
848  const LinearOperator<typename Range::BlockType,
849  typename Domain::BlockType,
850  typename BlockPayload::BlockType> &op)
851 {
852  static_assert(m > 0,
853  "a blockdiagonal LinearOperator must consist of at least "
854  "one block");
855 
856  using BlockType =
858  std::array<BlockType, m> new_ops;
859  new_ops.fill(op);
860 
861  return block_diagonal_operator(new_ops);
862 }
863 
864 
865 
867 
871 
907 template <typename Range, typename Domain, typename BlockPayload>
911  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
912 {
914  typename BlockPayload::BlockType(diagonal_inverse)};
915 
916  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
917  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
918 
919  return_op.vmult = [block_operator, diagonal_inverse](Range & v,
920  const Range &u) {
921  const unsigned int m = block_operator.n_block_rows();
922  Assert(block_operator.n_block_cols() == m,
923  ExcDimensionMismatch(block_operator.n_block_cols(), m));
924  Assert(diagonal_inverse.n_block_rows() == m,
925  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
926  Assert(diagonal_inverse.n_block_cols() == m,
927  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
928  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
929  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
930 
931  if (m == 0)
932  return;
933 
934  diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
935  for (unsigned int i = 1; i < m; ++i)
936  {
937  auto &dst = v.block(i);
938  dst = u.block(i);
939  dst *= -1.;
940  for (unsigned int j = 0; j < i; ++j)
941  block_operator.block(i, j).vmult_add(dst, v.block(j));
942  dst *= -1.;
943  diagonal_inverse.block(i, i).vmult(dst,
944  dst); // uses intermediate storage
945  }
946  };
947 
948  return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
949  const Range &u) {
950  const unsigned int m = block_operator.n_block_rows();
951  Assert(block_operator.n_block_cols() == m,
952  ExcDimensionMismatch(block_operator.n_block_cols(), m));
953  Assert(diagonal_inverse.n_block_rows() == m,
954  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
955  Assert(diagonal_inverse.n_block_cols() == m,
956  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
957  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
958  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
959 
960  if (m == 0)
961  return;
962 
965  vector_memory);
966 
967  diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
968 
969  for (unsigned int i = 1; i < m; ++i)
970  {
971  diagonal_inverse.block(i, i).reinit_range_vector(
972  *tmp, /*bool omit_zeroing_entries=*/true);
973  *tmp = u.block(i);
974  *tmp *= -1.;
975  for (unsigned int j = 0; j < i; ++j)
976  block_operator.block(i, j).vmult_add(*tmp, v.block(j));
977  *tmp *= -1.;
978  diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
979  }
980  };
981 
982  return return_op;
983 }
984 
985 
986 
1022 template <typename Range, typename Domain, typename BlockPayload>
1026  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
1027 {
1029  typename BlockPayload::BlockType(diagonal_inverse)};
1030 
1031  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
1032  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
1033 
1034  return_op.vmult = [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 
1046  if (m == 0)
1047  return;
1048 
1049  diagonal_inverse.block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1050 
1051  for (int i = m - 2; i >= 0; --i)
1052  {
1053  auto &dst = v.block(i);
1054  dst = u.block(i);
1055  dst *= -1.;
1056  for (unsigned int j = i + 1; j < m; ++j)
1057  block_operator.block(i, j).vmult_add(dst, v.block(j));
1058  dst *= -1.;
1059  diagonal_inverse.block(i, i).vmult(dst,
1060  dst); // uses intermediate storage
1061  }
1062  };
1063 
1064  return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
1065  const Range &u) {
1066  const unsigned int m = block_operator.n_block_rows();
1067  Assert(block_operator.n_block_cols() == m,
1068  ExcDimensionMismatch(block_operator.n_block_cols(), m));
1069  Assert(diagonal_inverse.n_block_rows() == m,
1070  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1071  Assert(diagonal_inverse.n_block_cols() == m,
1072  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1073  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1074  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1077  vector_memory);
1078 
1079  if (m == 0)
1080  return;
1081 
1082  diagonal_inverse.block(m - 1, m - 1)
1083  .vmult_add(v.block(m - 1), u.block(m - 1));
1084 
1085  for (int i = m - 2; i >= 0; --i)
1086  {
1087  diagonal_inverse.block(i, i).reinit_range_vector(
1088  *tmp, /*bool omit_zeroing_entries=*/true);
1089  *tmp = u.block(i);
1090  *tmp *= -1.;
1091  for (unsigned int j = i + 1; j < m; ++j)
1092  block_operator.block(i, j).vmult_add(*tmp, v.block(j));
1093  *tmp *= -1.;
1094  diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
1095  }
1096  };
1097 
1098  return return_op;
1099 }
1100 
1102 
1103 DEAL_II_NAMESPACE_CLOSE
1104 
1105 #endif
BlockLinearOperator(const std::array< BlockType, m > &ops)
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
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< BlockType, m > &ops)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &matrix)
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const Op &op)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::function< void(Range &v, const Domain &u)> vmult
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
BlockLinearOperator(const BlockPayload &payload)
std::function< unsigned int()> n_block_rows
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &, const BlockLinearOperator< Domain, Range, BlockPayload > &)
static bool equal(const T *p1, const T *p2)
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &, const BlockLinearOperator< Domain, Range, BlockPayload > &)