Reference documentation for deal.II version 9.1.1
\(\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 
30 namespace internal
31 {
32  namespace BlockLinearOperatorImplementation
33  {
34  template <typename PayloadBlockType =
37  }
38 } // namespace internal
39 
40 template <typename Number>
42 
43 template <typename Range = BlockVector<double>,
44  typename Domain = Range,
45  typename BlockPayload =
46  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
48 
49 template <typename Range = BlockVector<double>,
50  typename Domain = Range,
51  typename BlockPayload =
52  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
53  typename BlockMatrixType>
55 block_operator(const BlockMatrixType &matrix);
56 
57 template <std::size_t m,
58  std::size_t n,
59  typename Range = BlockVector<double>,
60  typename Domain = Range,
61  typename BlockPayload =
65  const std::array<std::array<LinearOperator<typename Range::BlockType,
66  typename Domain::BlockType,
67  typename BlockPayload::BlockType>,
68  n>,
69  m> &);
70 
71 template <std::size_t m,
72  typename Range = BlockVector<double>,
73  typename Domain = Range,
74  typename BlockPayload =
78  const std::array<LinearOperator<typename Range::BlockType,
79  typename Domain::BlockType,
80  typename BlockPayload::BlockType>,
81  m> &);
82 
83 template <std::size_t m,
84  typename Range = BlockVector<double>,
85  typename Domain = Range,
86  typename BlockPayload =
90  const LinearOperator<typename Range::BlockType,
91  typename Domain::BlockType,
92  typename BlockPayload::BlockType> &op);
93 
94 // This is a workaround for a bug in <=gcc-4.7 that does not like partial
95 // template default values in combination with local lambda expressions [1]
96 //
97 // [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624
98 //
99 // Forward declare functions with partial template defaults:
100 
101 template <typename Range = BlockVector<double>,
102  typename Domain = Range,
103  typename BlockPayload =
104  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
105  typename BlockMatrixType>
107 block_diagonal_operator(const BlockMatrixType &block_matrix);
108 
109 template <typename Range = BlockVector<double>,
110  typename Domain = Range,
111  typename BlockPayload =
112  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
117 
118 template <typename Range = BlockVector<double>,
119  typename Domain = Range,
120  typename BlockPayload =
121  internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
126 
127 // end of workaround
128 
129 
130 
200 template <typename Range, typename Domain, typename BlockPayload>
202  : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
203 {
204 public:
205  using BlockType = LinearOperator<typename Range::BlockType,
206  typename Domain::BlockType,
207  typename BlockPayload::BlockType>;
208 
216  BlockLinearOperator(const BlockPayload &payload)
217  : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(
218  typename BlockPayload::BlockType(payload, payload))
219  {
220  n_block_rows = []() -> unsigned int {
221  Assert(
222  false,
223  ExcMessage(
224  "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
225  return 0;
226  };
227 
228  n_block_cols = []() -> unsigned int {
229  Assert(
230  false,
231  ExcMessage(
232  "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
233  return 0;
234  };
235 
236  block = [](unsigned int, unsigned int) -> BlockType {
237  Assert(
238  false,
239  ExcMessage(
240  "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
241  return BlockType();
242  };
243  }
244 
250 
256  template <typename Op>
257  BlockLinearOperator(const Op &op)
258  {
259  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
260  }
261 
267  template <std::size_t m, std::size_t n>
268  BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
269  {
270  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
271  }
272 
278  template <std::size_t m>
279  BlockLinearOperator(const std::array<BlockType, m> &ops)
280  {
281  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
282  }
283 
289 
294  template <typename Op>
296  operator=(const Op &op)
297  {
298  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
299  return *this;
300  }
301 
307  template <std::size_t m, std::size_t n>
309  operator=(const std::array<std::array<BlockType, n>, m> &ops)
310  {
311  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
312  return *this;
313  }
314 
320  template <std::size_t m>
322  operator=(const std::array<BlockType, m> &ops)
323  {
324  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
325  return *this;
326  }
327 
332  std::function<unsigned int()> n_block_rows;
333 
338  std::function<unsigned int()> n_block_cols;
339 
345  std::function<BlockType(unsigned int, unsigned int)> block;
346 };
347 
348 
349 namespace internal
350 {
351  namespace BlockLinearOperatorImplementation
352  {
353  // A helper function to apply a given vmult, or Tvmult to a vector with
354  // intermediate storage, similar to the corresponding helper
355  // function for LinearOperator. Here, two operators are used.
356  // The first one takes care of the first "column" and typically doesn't add.
357  // On the other hand, the second operator is normally an adding one.
358  template <typename Function1,
359  typename Function2,
360  typename Range,
361  typename Domain>
362  void
363  apply_with_intermediate_storage(const Function1 &first_op,
364  const Function2 &loop_op,
365  Range & v,
366  const Domain & u,
367  bool add)
368  {
369  GrowingVectorMemory<Range> vector_memory;
370 
371  typename VectorMemory<Range>::Pointer tmp(vector_memory);
372  tmp->reinit(v, /*bool omit_zeroing_entries =*/true);
373 
374  const unsigned int n = u.n_blocks();
375  const unsigned int m = v.n_blocks();
376 
377  for (unsigned int i = 0; i < m; ++i)
378  {
379  first_op(*tmp, u, i, 0);
380  for (unsigned int j = 1; j < n; ++j)
381  loop_op(*tmp, u, i, j);
382  }
383 
384  if (add)
385  v += *tmp;
386  else
387  v = *tmp;
388  }
389 
390  // Populate the LinearOperator interfaces with the help of the
391  // BlockLinearOperator functions
392  template <typename Range, typename Domain, typename BlockPayload>
393  inline void
394  populate_linear_operator_functions(
396  {
397  op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries) {
398  const unsigned int m = op.n_block_rows();
399 
400  // Reinitialize the block vector to m blocks:
401  v.reinit(m);
402 
403  // And reinitialize every individual block with reinit_range_vectors:
404  for (unsigned int i = 0; i < m; ++i)
405  op.block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
406 
407  v.collect_sizes();
408  };
409 
410  op.reinit_domain_vector = [=](Domain &v, bool omit_zeroing_entries) {
411  const unsigned int n = op.n_block_cols();
412 
413  // Reinitialize the block vector to n blocks:
414  v.reinit(n);
415 
416  // And reinitialize every individual block with reinit_domain_vectors:
417  for (unsigned int i = 0; i < n; ++i)
418  op.block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
419 
420  v.collect_sizes();
421  };
422 
423  op.vmult = [&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, false);
446  }
447  else
448  {
449  for (unsigned int i = 0; i < m; ++i)
450  {
451  op.block(i, 0).vmult(v.block(i), u.block(0));
452  for (unsigned int j = 1; j < n; ++j)
453  op.block(i, j).vmult_add(v.block(i), u.block(j));
454  }
455  }
456  };
457 
458  op.vmult_add = [&op](Range &v, const Domain &u) {
459  const unsigned int m = op.n_block_rows();
460  const unsigned int n = op.n_block_cols();
461  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
462  Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
463 
464  if (PointerComparison::equal(&v, &u))
465  {
466  const auto first_op = [&op](Range & v,
467  const Domain & u,
468  const unsigned int i,
469  const unsigned int j) {
470  op.block(i, j).vmult(v.block(i), u.block(j));
471  };
472 
473  const auto loop_op = [&op](Range & v,
474  const Domain & u,
475  const unsigned int i,
476  const unsigned int j) {
477  op.block(i, j).vmult_add(v.block(i), u.block(j));
478  };
479 
480  apply_with_intermediate_storage(first_op, loop_op, v, u, true);
481  }
482  else
483  {
484  for (unsigned int i = 0; i < m; ++i)
485  for (unsigned int j = 0; j < n; ++j)
486  op.block(i, j).vmult_add(v.block(i), u.block(j));
487  }
488  };
489 
490  op.Tvmult = [&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, false);
513  }
514  else
515  {
516  for (unsigned int i = 0; i < n; ++i)
517  {
518  op.block(0, i).Tvmult(v.block(i), u.block(0));
519  for (unsigned int j = 1; j < m; ++j)
520  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
521  }
522  }
523  };
524 
525  op.Tvmult_add = [&op](Domain &v, const Range &u) {
526  const unsigned int n = op.n_block_cols();
527  const unsigned int m = op.n_block_rows();
528  Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
529  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
530 
531  if (PointerComparison::equal(&v, &u))
532  {
533  const auto first_op = [&op](Range & v,
534  const Domain & u,
535  const unsigned int i,
536  const unsigned int j) {
537  op.block(j, i).Tvmult(v.block(i), u.block(j));
538  };
539 
540  const auto loop_op = [&op](Range & v,
541  const Domain & u,
542  const unsigned int i,
543  const unsigned int j) {
544  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
545  };
546 
547  apply_with_intermediate_storage(first_op, loop_op, v, u, true);
548  }
549  else
550  {
551  for (unsigned int i = 0; i < n; ++i)
552  for (unsigned int j = 0; j < m; ++j)
553  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
554  }
555  };
556  }
557 
558 
559 
574  template <typename PayloadBlockType>
575  class EmptyBlockPayload
576  {
577  public:
581  using BlockType = PayloadBlockType;
582 
590  template <typename... Args>
591  EmptyBlockPayload(const Args &...)
592  {}
593  };
594 
595  } // namespace BlockLinearOperatorImplementation
596 } // namespace internal
597 
598 
599 
604 
616 template <typename Range,
617  typename Domain,
618  typename BlockPayload,
619  typename BlockMatrixType>
621 block_operator(const BlockMatrixType &block_matrix)
622 {
623  using BlockType =
625 
627  BlockPayload(block_matrix, block_matrix)};
628 
629  return_op.n_block_rows = [&block_matrix]() -> unsigned int {
630  return block_matrix.n_block_rows();
631  };
632 
633  return_op.n_block_cols = [&block_matrix]() -> unsigned int {
634  return block_matrix.n_block_cols();
635  };
636 
637  return_op.block = [&block_matrix](unsigned int i,
638  unsigned int j) -> BlockType {
639 #ifdef DEBUG
640  const unsigned int m = block_matrix.n_block_rows();
641  const unsigned int n = block_matrix.n_block_cols();
642  Assert(i < m, ExcIndexRange(i, 0, m));
643  Assert(j < n, ExcIndexRange(j, 0, n));
644 #endif
645 
646  return BlockType(block_matrix.block(i, j));
647  };
648 
649  populate_linear_operator_functions(return_op);
650  return return_op;
651 }
652 
653 
654 
682 template <std::size_t m,
683  std::size_t n,
684  typename Range,
685  typename Domain,
686  typename BlockPayload>
689  const std::array<std::array<LinearOperator<typename Range::BlockType,
690  typename Domain::BlockType,
691  typename BlockPayload::BlockType>,
692  n>,
693  m> &ops)
694 {
695  static_assert(m > 0 && n > 0,
696  "a blocked LinearOperator must consist of at least one block");
697 
698  using BlockType =
700 
701  // TODO: Create block payload so that this can be initialized correctly
702  BlockLinearOperator<Range, Domain, BlockPayload> return_op{BlockPayload()};
703 
704  return_op.n_block_rows = []() -> unsigned int { return m; };
705 
706  return_op.n_block_cols = []() -> unsigned int { return n; };
707 
708  return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType {
709  Assert(i < m, ExcIndexRange(i, 0, m));
710  Assert(j < n, ExcIndexRange(j, 0, n));
711 
712  return ops[i][j];
713  };
714 
715  populate_linear_operator_functions(return_op);
716  return return_op;
717 }
718 
719 
720 
736 template <typename Range,
737  typename Domain,
738  typename BlockPayload,
739  typename BlockMatrixType>
741 block_diagonal_operator(const BlockMatrixType &block_matrix)
742 {
743  using BlockType =
745 
747  BlockPayload(block_matrix, block_matrix)};
748 
749  return_op.n_block_rows = [&block_matrix]() -> unsigned int {
750  return block_matrix.n_block_rows();
751  };
752 
753  return_op.n_block_cols = [&block_matrix]() -> unsigned int {
754  return block_matrix.n_block_cols();
755  };
756 
757  return_op.block = [&block_matrix](unsigned int i,
758  unsigned int j) -> BlockType {
759 #ifdef DEBUG
760  const unsigned int m = block_matrix.n_block_rows();
761  const unsigned int n = block_matrix.n_block_cols();
762  Assert(m == n, ExcDimensionMismatch(m, n));
763  Assert(i < m, ExcIndexRange(i, 0, m));
764  Assert(j < n, ExcIndexRange(j, 0, n));
765 #endif
766  if (i == j)
767  return BlockType(block_matrix.block(i, j));
768  else
769  return null_operator(BlockType(block_matrix.block(i, j)));
770  };
771 
772  populate_linear_operator_functions(return_op);
773  return return_op;
774 }
775 
776 
777 
795 template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
798  const std::array<LinearOperator<typename Range::BlockType,
799  typename Domain::BlockType,
800  typename BlockPayload::BlockType>,
801  m> &ops)
802 {
803  static_assert(
804  m > 0, "a blockdiagonal LinearOperator must consist of at least one block");
805 
806  using BlockType =
808 
809  std::array<std::array<BlockType, m>, m> new_ops;
810 
811  // This is a bit tricky. We have to make sure that the off-diagonal
812  // elements of return_op.ops are populated correctly. They must be
813  // null_operators, but with correct reinit_domain_vector and
814  // reinit_range_vector functions.
815  for (unsigned int i = 0; i < m; ++i)
816  for (unsigned int j = 0; j < m; ++j)
817  if (i == j)
818  {
819  // diagonal elements are easy:
820  new_ops[i][j] = ops[i];
821  }
822  else
823  {
824  // create a null-operator...
825  new_ops[i][j] = null_operator(ops[i]);
826  // ... and fix up reinit_domain_vector:
827  new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
828  }
829 
830  return block_operator<m, m, Range, Domain>(new_ops);
831 }
832 
833 
834 
844 template <std::size_t m, typename Range, typename Domain, typename BlockPayload>
847  const LinearOperator<typename Range::BlockType,
848  typename Domain::BlockType,
849  typename BlockPayload::BlockType> &op)
850 {
851  static_assert(m > 0,
852  "a blockdiagonal LinearOperator must consist of at least "
853  "one block");
854 
855  using BlockType =
857  std::array<BlockType, m> new_ops;
858  new_ops.fill(op);
859 
860  return block_diagonal_operator(new_ops);
861 }
862 
863 
864 
866 
870 
906 template <typename Range, typename Domain, typename BlockPayload>
910  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
911 {
913  typename BlockPayload::BlockType(diagonal_inverse)};
914 
915  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
916  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
917 
918  return_op.vmult = [block_operator, diagonal_inverse](Range & v,
919  const Range &u) {
920  const unsigned int m = block_operator.n_block_rows();
921  Assert(block_operator.n_block_cols() == m,
922  ExcDimensionMismatch(block_operator.n_block_cols(), m));
923  Assert(diagonal_inverse.n_block_rows() == m,
924  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
925  Assert(diagonal_inverse.n_block_cols() == m,
926  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
927  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
928  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
929 
930  if (m == 0)
931  return;
932 
933  diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
934  for (unsigned int i = 1; i < m; ++i)
935  {
936  auto &dst = v.block(i);
937  dst = u.block(i);
938  dst *= -1.;
939  for (unsigned int j = 0; j < i; ++j)
940  block_operator.block(i, j).vmult_add(dst, v.block(j));
941  dst *= -1.;
942  diagonal_inverse.block(i, i).vmult(dst,
943  dst); // uses intermediate storage
944  }
945  };
946 
947  return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
948  const Range &u) {
949  const unsigned int m = block_operator.n_block_rows();
950  Assert(block_operator.n_block_cols() == m,
951  ExcDimensionMismatch(block_operator.n_block_cols(), m));
952  Assert(diagonal_inverse.n_block_rows() == m,
953  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
954  Assert(diagonal_inverse.n_block_cols() == m,
955  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
956  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
957  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
958 
959  if (m == 0)
960  return;
961 
964  vector_memory);
965 
966  diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
967 
968  for (unsigned int i = 1; i < m; ++i)
969  {
970  diagonal_inverse.block(i, i).reinit_range_vector(
971  *tmp, /*bool omit_zeroing_entries=*/true);
972  *tmp = u.block(i);
973  *tmp *= -1.;
974  for (unsigned int j = 0; j < i; ++j)
975  block_operator.block(i, j).vmult_add(*tmp, v.block(j));
976  *tmp *= -1.;
977  diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
978  }
979  };
980 
981  return return_op;
982 }
983 
984 
985 
1021 template <typename Range, typename Domain, typename BlockPayload>
1025  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
1026 {
1028  typename BlockPayload::BlockType(diagonal_inverse)};
1029 
1030  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
1031  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
1032 
1033  return_op.vmult = [block_operator, diagonal_inverse](Range & v,
1034  const Range &u) {
1035  const unsigned int m = block_operator.n_block_rows();
1036  Assert(block_operator.n_block_cols() == m,
1037  ExcDimensionMismatch(block_operator.n_block_cols(), m));
1038  Assert(diagonal_inverse.n_block_rows() == m,
1039  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1040  Assert(diagonal_inverse.n_block_cols() == m,
1041  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1042  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1043  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1044 
1045  if (m == 0)
1046  return;
1047 
1048  diagonal_inverse.block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1049 
1050  for (int i = m - 2; i >= 0; --i)
1051  {
1052  auto &dst = v.block(i);
1053  dst = u.block(i);
1054  dst *= -1.;
1055  for (unsigned int j = i + 1; j < m; ++j)
1056  block_operator.block(i, j).vmult_add(dst, v.block(j));
1057  dst *= -1.;
1058  diagonal_inverse.block(i, i).vmult(dst,
1059  dst); // uses intermediate storage
1060  }
1061  };
1062 
1063  return_op.vmult_add = [block_operator, diagonal_inverse](Range & v,
1064  const Range &u) {
1065  const unsigned int m = block_operator.n_block_rows();
1066  Assert(block_operator.n_block_cols() == m,
1067  ExcDimensionMismatch(block_operator.n_block_cols(), m));
1068  Assert(diagonal_inverse.n_block_rows() == m,
1069  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
1070  Assert(diagonal_inverse.n_block_cols() == m,
1071  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
1072  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
1073  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
1076  vector_memory);
1077 
1078  if (m == 0)
1079  return;
1080 
1081  diagonal_inverse.block(m - 1, m - 1)
1082  .vmult_add(v.block(m - 1), u.block(m - 1));
1083 
1084  for (int i = m - 2; i >= 0; --i)
1085  {
1086  diagonal_inverse.block(i, i).reinit_range_vector(
1087  *tmp, /*bool omit_zeroing_entries=*/true);
1088  *tmp = u.block(i);
1089  *tmp *= -1.;
1090  for (unsigned int j = i + 1; j < m; ++j)
1091  block_operator.block(i, j).vmult_add(*tmp, v.block(j));
1092  *tmp *= -1.;
1093  diagonal_inverse.block(i, i).vmult_add(v.block(i), *tmp);
1094  }
1095  };
1096 
1097  return return_op;
1098 }
1099 
1101 
1102 DEAL_II_NAMESPACE_CLOSE
1103 
1104 #endif
BlockLinearOperator(const std::array< BlockType, m > &ops)
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
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< BlockType, m > &ops)
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 > &op)
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< std::array< BlockType, n >, m > &ops)
static ::ExceptionBase & ExcMessage(std::string arg1)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &block_matrix)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const BlockLinearOperator< Range, Domain, BlockPayload > &)=default
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
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
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)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const BlockMatrixType &block_matrix)
std::function< unsigned int()> n_block_rows
static bool equal(const T *p1, const T *p2)