Reference documentation for deal.II version Git 2618e0f 2017-11-23 17:25:26 +0100
block_linear_operator.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2010 - 2017 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 at
12 // the top level of the deal.II distribution.
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 #include <deal.II/base/exceptions.h>
21 
22 #include <deal.II/lac/linear_operator.h>
23 
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 // Forward declarations:
28 
29 namespace internal
30 {
31  namespace BlockLinearOperator
32  {
33  template <typename PayloadBlockType = internal::LinearOperator::EmptyPayload>
35  }
36 }
37 
38 template <typename Number> class BlockVector;
39 
40 template <typename Range = BlockVector<double>,
41  typename Domain = Range,
42  typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
44 
45 template <typename Range = BlockVector<double>,
46  typename Domain = Range,
47  typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<>,
48  typename BlockMatrixType>
50 block_operator(const BlockMatrixType &matrix);
51 
52 template <size_t m, size_t n,
53  typename Range = BlockVector<double>,
54  typename Domain = Range,
58 
59 template <size_t m,
60  typename Range = BlockVector<double>,
61  typename Domain = Range,
65 
66 template <size_t m,
67  typename Range = BlockVector<double>,
68  typename Domain = Range,
72 
73 // This is a workaround for a bug in <=gcc-4.7 that does not like partial
74 // template default values in combination with local lambda expressions [1]
75 //
76 // [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624
77 //
78 // Forward declare functions with partial template defaults:
79 
80 template <typename Range = BlockVector<double>,
81  typename Domain = Range,
82  typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<>,
83  typename BlockMatrixType>
85 block_diagonal_operator(const BlockMatrixType &block_matrix);
86 
87 template <typename Range = BlockVector<double>,
88  typename Domain = Range,
89  typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
93 
94 template <typename Range = BlockVector<double>,
95  typename Domain = Range,
96  typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
100 
101 // end of workaround
102 
103 
104 
173 template <typename Range, typename Domain, typename BlockPayload>
174 class BlockLinearOperator : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
175 {
176 public:
177 
179 
187  BlockLinearOperator(const BlockPayload &payload)
188  : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(typename BlockPayload::BlockType(payload,payload))
189  {
190  n_block_rows = []() -> unsigned int
191  {
192  Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
193  return 0;
194  };
195 
196  n_block_cols = []() -> unsigned int
197  {
198  Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
199  return 0;
200  };
201 
202  block = [](unsigned int, unsigned int) -> BlockType
203  {
204  Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::block called"));
205  return BlockType();
206  };
207  }
208 
213  default;
214 
220  template <typename Op>
221  BlockLinearOperator(const Op &op)
222  {
223  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
224  }
225 
231  template <size_t m, size_t n>
232  BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
233  {
234  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
235  }
236 
242  template <size_t m>
243  BlockLinearOperator(const std::array<BlockType, m> &ops)
244  {
245  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
246  }
247 
253 
258  template <typename Op>
260  {
261  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
262  return *this;
263  }
264 
270  template <size_t m, size_t n>
272  operator=(const std::array<std::array<BlockType, n>, m> &ops)
273  {
274  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
275  return *this;
276  }
277 
283  template <size_t m>
285  operator=(const std::array<BlockType, m> &ops)
286  {
287  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
288  return *this;
289  }
290 
295  std::function<unsigned int()> n_block_rows;
296 
301  std::function<unsigned int()> n_block_cols;
302 
308  std::function<BlockType(unsigned int, unsigned int)> block;
309 };
310 
311 
312 
313 namespace internal
314 {
315  namespace BlockLinearOperator
316  {
317  // Populate the LinearOperator interfaces with the help of the
318  // BlockLinearOperator functions
319  template <typename Range, typename Domain, typename BlockPayload>
320  inline void
321  populate_linear_operator_functions(
323  {
324  op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries)
325  {
326  const unsigned int m = op.n_block_rows();
327 
328  // Reinitialize the block vector to m blocks:
329  v.reinit(m);
330 
331  // And reinitialize every individual block with reinit_range_vectors:
332  for (unsigned int i = 0; i < m; ++i)
333  op.block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
334 
335  v.collect_sizes();
336  };
337 
338  op.reinit_domain_vector = [=](Domain &v, bool omit_zeroing_entries)
339  {
340  const unsigned int n = op.n_block_cols();
341 
342  // Reinitialize the block vector to n blocks:
343  v.reinit(n);
344 
345  // And reinitialize every individual block with reinit_domain_vectors:
346  for (unsigned int i = 0; i < n; ++i)
347  op.block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
348 
349  v.collect_sizes();
350  };
351 
352  op.vmult = [=](Range &v, const Domain &u)
353  {
354  const unsigned int m = op.n_block_rows();
355  const unsigned int n = op.n_block_cols();
356  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
357  Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
358 
359  for (unsigned int i = 0; i < m; ++i)
360  {
361  op.block(i, 0).vmult(v.block(i), u.block(0));
362  for (unsigned int j = 1; j < n; ++j)
363  op.block(i, j).vmult_add(v.block(i), u.block(j));
364  }
365  };
366 
367  op.vmult_add = [=](Range &v, const Domain &u)
368  {
369  const unsigned int m = op.n_block_rows();
370  const unsigned int n = op.n_block_cols();
371  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
372  Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
373 
374  for (unsigned int i = 0; i < m; ++i)
375  for (unsigned int j = 0; j < n; ++j)
376  op.block(i, j).vmult_add(v.block(i), u.block(j));
377  };
378 
379  op.Tvmult = [=](Domain &v, const Range &u)
380  {
381  const unsigned int n = op.n_block_cols();
382  const unsigned int m = op.n_block_rows();
383  Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
384  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
385 
386  for (unsigned int i = 0; i < n; ++i)
387  {
388  op.block(0, i).Tvmult(v.block(i), u.block(0));
389  for (unsigned int j = 1; j < m; ++j)
390  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
391  }
392  };
393 
394  op.Tvmult_add = [=](Domain &v, const Range &u)
395  {
396  const unsigned int n = op.n_block_cols();
397  const unsigned int m = op.n_block_rows();
398  Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
399  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
400 
401  for (unsigned int i = 0; i < n; ++i)
402  for (unsigned int j = 0; j < m; ++j)
403  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
404  };
405  }
406 
407 
422  template <typename PayloadBlockType>
423  class EmptyBlockPayload
424  {
425  public:
429  typedef PayloadBlockType BlockType;
430 
438  template <typename... Args>
439  EmptyBlockPayload (const Args &...)
440  { }
441  };
442 
443  } /*namespace BlockLinearOperator*/
444 } /*namespace internal*/
445 
446 
447 
452 
464 template <typename Range,
465  typename Domain,
466  typename BlockPayload,
467  typename BlockMatrixType>
469 block_operator(const BlockMatrixType &block_matrix)
470 {
472 
473  BlockLinearOperator<Range, Domain, BlockPayload> return_op (BlockPayload(block_matrix,block_matrix));
474 
475  return_op.n_block_rows = [&block_matrix]() -> unsigned int
476  {
477  return block_matrix.n_block_rows();
478  };
479 
480  return_op.n_block_cols = [&block_matrix]() -> unsigned int
481  {
482  return block_matrix.n_block_cols();
483  };
484 
485  return_op.block = [&block_matrix](unsigned int i, unsigned int j) -> BlockType
486  {
487 #ifdef DEBUG
488  const unsigned int m = block_matrix.n_block_rows();
489  const unsigned int n = block_matrix.n_block_cols();
490  Assert(i < m, ExcIndexRange (i, 0, m));
491  Assert(j < n, ExcIndexRange (j, 0, n));
492 #endif
493 
494  return BlockType(block_matrix.block(i, j));
495  };
496 
497  internal::BlockLinearOperator::populate_linear_operator_functions(return_op);
498  return return_op;
499 }
500 
501 
502 
530 template <size_t m, size_t n, typename Range, typename Domain, typename BlockPayload>
533 {
534  static_assert(m > 0 && n > 0,
535  "a blocked LinearOperator must consist of at least one block");
536 
538 
539  BlockLinearOperator<Range, Domain, BlockPayload> return_op ((BlockPayload())); // TODO: Create block payload so that this can be initialized correctly
540 
541  return_op.n_block_rows = []() -> unsigned int
542  {
543  return m;
544  };
545 
546  return_op.n_block_cols = []() -> unsigned int
547  {
548  return n;
549  };
550 
551  return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType
552  {
553  Assert(i < m, ExcIndexRange (i, 0, m));
554  Assert(j < n, ExcIndexRange (j, 0, n));
555 
556  return ops[i][j];
557  };
558 
559  internal::BlockLinearOperator::populate_linear_operator_functions(return_op);
560  return return_op;
561 }
562 
563 
564 
580 template <typename Range,
581  typename Domain,
582  typename BlockPayload,
583  typename BlockMatrixType>
585 block_diagonal_operator(const BlockMatrixType &block_matrix)
586 {
588 
589  BlockLinearOperator<Range, Domain, BlockPayload> return_op (BlockPayload(block_matrix,block_matrix));
590 
591  return_op.n_block_rows = [&block_matrix]() -> unsigned int
592  {
593  return block_matrix.n_block_rows();
594  };
595 
596  return_op.n_block_cols = [&block_matrix]() -> unsigned int
597  {
598  return block_matrix.n_block_cols();
599  };
600 
601  return_op.block = [&block_matrix](unsigned int i, unsigned int j) -> BlockType
602  {
603 #ifdef DEBUG
604  const unsigned int m = block_matrix.n_block_rows();
605  const unsigned int n = block_matrix.n_block_cols();
606  Assert(m == n, ExcDimensionMismatch(m, n));
607  Assert(i < m, ExcIndexRange (i, 0, m));
608  Assert(j < n, ExcIndexRange (j, 0, n));
609 #endif
610  if (i == j)
611  return BlockType(block_matrix.block(i, j));
612  else
613  return null_operator(BlockType(block_matrix.block(i, j)));
614  };
615 
616  internal::BlockLinearOperator::populate_linear_operator_functions(return_op);
617  return return_op;
618 }
619 
620 
621 
639 template <size_t m, typename Range, typename Domain, typename BlockPayload>
642 {
643  static_assert(m > 0,
644  "a blockdiagonal LinearOperator must consist of at least one block");
645 
647 
648  std::array<std::array<BlockType, m>, m> new_ops;
649 
650  // This is a bit tricky. We have to make sure that the off-diagonal
651  // elements of return_op.ops are populated correctly. They must be
652  // null_operators, but with correct reinit_domain_vector and
653  // reinit_range_vector functions.
654  for (unsigned int i = 0; i < m; ++i)
655  for (unsigned int j = 0; j < m; ++j)
656  if (i == j)
657  {
658  // diagonal elements are easy:
659  new_ops[i][j] = ops[i];
660  }
661  else
662  {
663  // create a null-operator...
664  new_ops[i][j] = null_operator(ops[i]);
665  // ... and fix up reinit_domain_vector:
666  new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
667  }
668 
669  return block_operator<m,m,Range,Domain>(new_ops);
670 }
671 
672 
673 
683 template <size_t m, typename Range, typename Domain, typename BlockPayload>
686 {
687  static_assert(m > 0,
688  "a blockdiagonal LinearOperator must consist of at least "
689  "one block");
690 
692  std::array<BlockType, m> new_ops;
693  new_ops.fill(op);
694 
695  return block_diagonal_operator(new_ops);
696 }
697 
698 
699 
701 
705 
741 template <typename Range, typename Domain, typename BlockPayload>
744  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
745 {
746  LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op ((typename BlockPayload::BlockType(diagonal_inverse)));
747 
748  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
749  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
750 
751  return_op.vmult = [block_operator, diagonal_inverse](Range &v, const Range &u)
752  {
753  const unsigned int m = block_operator.n_block_rows();
754  Assert(block_operator.n_block_cols() == m,
755  ExcDimensionMismatch(block_operator.n_block_cols(), m));
756  Assert(diagonal_inverse.n_block_rows() == m,
757  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
758  Assert(diagonal_inverse.n_block_cols() == m,
759  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
760  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
761  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
762 
763  if (m == 0)
764  return;
765 
766  diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
767  for (unsigned int i = 1; i < m; ++i)
768  {
769  auto &dst = v.block(i);
770  dst = u.block(i);
771  dst *= -1.;
772  for (unsigned int j = 0; j < i; ++j)
773  block_operator.block(i, j).vmult_add(dst, v.block(j));
774  dst *= -1.;
775  diagonal_inverse.block(i, i).vmult(dst, dst); // uses intermediate storage
776  }
777  };
778 
779  return_op.vmult_add = [block_operator, diagonal_inverse](Range &v, const Range &u)
780  {
781  const unsigned int m = block_operator.n_block_rows();
782  Assert(block_operator.n_block_cols() == m,
783  ExcDimensionMismatch(block_operator.n_block_cols(), m));
784  Assert(diagonal_inverse.n_block_rows() == m,
785  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
786  Assert(diagonal_inverse.n_block_cols() == m,
787  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
788  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
789  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
790 
791  if (m == 0)
792  return;
793 
795  typename VectorMemory<typename Range::BlockType>::Pointer tmp (vector_memory);
796 
797  diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
798 
799  for (unsigned int i = 1; i < m; ++i)
800  {
801  diagonal_inverse.block(i, i).reinit_range_vector(*tmp, /*bool omit_zeroing_entries=*/ true);
802  *tmp = u.block(i);
803  *tmp *= -1.;
804  for (unsigned int j = 0; j < i; ++j)
805  block_operator.block(i, j).vmult_add(*tmp, v.block(j));
806  *tmp *= -1.;
807  diagonal_inverse.block(i, i).vmult_add(v.block(i),*tmp);
808  }
809  };
810 
811  return return_op;
812 }
813 
814 
815 
851 template <typename Range, typename Domain, typename BlockPayload>
854  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
855 {
856  LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op ((typename BlockPayload::BlockType(diagonal_inverse)));
857 
858  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
859  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
860 
861  return_op.vmult = [block_operator, diagonal_inverse](Range &v, const Range &u)
862  {
863  const unsigned int m = block_operator.n_block_rows();
864  Assert(block_operator.n_block_cols() == m,
865  ExcDimensionMismatch(block_operator.n_block_cols(), m));
866  Assert(diagonal_inverse.n_block_rows() == m,
867  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
868  Assert(diagonal_inverse.n_block_cols() == m,
869  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
870  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
871  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
872 
873  if (m == 0)
874  return;
875 
876  diagonal_inverse.block(m-1, m-1).vmult(v.block(m-1),u.block(m-1));
877 
878  for (int i = m - 2; i >= 0; --i)
879  {
880  auto &dst = v.block(i);
881  dst = u.block(i);
882  dst *= -1.;
883  for (unsigned int j = i + 1; j < m; ++j)
884  block_operator.block(i, j).vmult_add(dst, v.block(j));
885  dst *= -1.;
886  diagonal_inverse.block(i, i).vmult(dst, dst); // uses intermediate storage
887  }
888  };
889 
890  return_op.vmult_add = [block_operator, diagonal_inverse](Range &v, const Range &u)
891  {
892  const unsigned int m = block_operator.n_block_rows();
893  Assert(block_operator.n_block_cols() == m,
894  ExcDimensionMismatch(block_operator.n_block_cols(), m));
895  Assert(diagonal_inverse.n_block_rows() == m,
896  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
897  Assert(diagonal_inverse.n_block_cols() == m,
898  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
899  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
900  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
902  typename VectorMemory<typename Range::BlockType>::Pointer tmp (vector_memory);
903 
904  if (m == 0)
905  return;
906 
907  diagonal_inverse.block(m-1, m-1).vmult_add(v.block(m-1),u.block(m-1));
908 
909  for (int i = m - 2; i >= 0; --i)
910  {
911  diagonal_inverse.block(i, i).reinit_range_vector(*tmp, /*bool omit_zeroing_entries=*/ true);
912  *tmp = u.block(i);
913  *tmp *= -1.;
914  for (unsigned int j = i + 1; j < m; ++j)
915  block_operator.block(i, j).vmult_add(*tmp,v.block(j));
916  *tmp *= -1.;
917  diagonal_inverse.block(i, i).vmult_add(v.block(i),*tmp);
918  }
919  };
920 
921  return return_op;
922 }
923 
925 
926 DEAL_II_NAMESPACE_CLOSE
927 
928 #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
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, n >, m > &ops)
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:337
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
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > &op)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &ops)