Reference documentation for deal.II version Git 7f58f27 2017-09-20 15:35:49 -0600
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 
140 template <typename Range, typename Domain, typename BlockPayload>
141 class BlockLinearOperator : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
142 {
143 public:
144 
146 
154  BlockLinearOperator(const BlockPayload &payload)
155  : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(typename BlockPayload::BlockType(payload,payload))
156  {
157  n_block_rows = []() -> unsigned int
158  {
159  Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
160  return 0;
161  };
162 
163  n_block_cols = []() -> unsigned int
164  {
165  Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
166  return 0;
167  };
168 
169  block = [](unsigned int, unsigned int) -> BlockType
170  {
171  Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::block called"));
172  return BlockType();
173  };
174  }
175 
180  default;
181 
187  template <typename Op>
188  BlockLinearOperator(const Op &op)
189  {
190  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
191  }
192 
198  template <size_t m, size_t n>
199  BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
200  {
201  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
202  }
203 
209  template <size_t m>
210  BlockLinearOperator(const std::array<BlockType, m> &ops)
211  {
212  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
213  }
214 
220 
225  template <typename Op>
227  {
228  *this = block_operator<Range, Domain, BlockPayload, Op>(op);
229  return *this;
230  }
231 
237  template <size_t m, size_t n>
239  operator=(const std::array<std::array<BlockType, n>, m> &ops)
240  {
241  *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
242  return *this;
243  }
244 
250  template <size_t m>
252  operator=(const std::array<BlockType, m> &ops)
253  {
254  *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
255  return *this;
256  }
257 
262  std::function<unsigned int()> n_block_rows;
263 
268  std::function<unsigned int()> n_block_cols;
269 
275  std::function<BlockType(unsigned int, unsigned int)> block;
276 };
277 
278 
279 
280 namespace internal
281 {
282  namespace BlockLinearOperator
283  {
284  // Populate the LinearOperator interfaces with the help of the
285  // BlockLinearOperator functions
286  template <typename Range, typename Domain, typename BlockPayload>
287  inline void
288  populate_linear_operator_functions(
290  {
291  op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries)
292  {
293  const unsigned int m = op.n_block_rows();
294 
295  // Reinitialize the block vector to m blocks:
296  v.reinit(m);
297 
298  // And reinitialize every individual block with reinit_range_vectors:
299  for (unsigned int i = 0; i < m; ++i)
300  op.block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
301 
302  v.collect_sizes();
303  };
304 
305  op.reinit_domain_vector = [=](Domain &v, bool omit_zeroing_entries)
306  {
307  const unsigned int n = op.n_block_cols();
308 
309  // Reinitialize the block vector to n blocks:
310  v.reinit(n);
311 
312  // And reinitialize every individual block with reinit_domain_vectors:
313  for (unsigned int i = 0; i < n; ++i)
314  op.block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
315 
316  v.collect_sizes();
317  };
318 
319  op.vmult = [=](Range &v, const Domain &u)
320  {
321  const unsigned int m = op.n_block_rows();
322  const unsigned int n = op.n_block_cols();
323  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
324  Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
325 
326  for (unsigned int i = 0; i < m; ++i)
327  {
328  op.block(i, 0).vmult(v.block(i), u.block(0));
329  for (unsigned int j = 1; j < n; ++j)
330  op.block(i, j).vmult_add(v.block(i), u.block(j));
331  }
332  };
333 
334  op.vmult_add = [=](Range &v, const Domain &u)
335  {
336  const unsigned int m = op.n_block_rows();
337  const unsigned int n = op.n_block_cols();
338  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
339  Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n));
340 
341  for (unsigned int i = 0; i < m; ++i)
342  for (unsigned int j = 0; j < n; ++j)
343  op.block(i, j).vmult_add(v.block(i), u.block(j));
344  };
345 
346  op.Tvmult = [=](Domain &v, const Range &u)
347  {
348  const unsigned int n = op.n_block_cols();
349  const unsigned int m = op.n_block_rows();
350  Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
351  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
352 
353  for (unsigned int i = 0; i < n; ++i)
354  {
355  op.block(0, i).Tvmult(v.block(i), u.block(0));
356  for (unsigned int j = 1; j < m; ++j)
357  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
358  }
359  };
360 
361  op.Tvmult_add = [=](Domain &v, const Range &u)
362  {
363  const unsigned int n = op.n_block_cols();
364  const unsigned int m = op.n_block_rows();
365  Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n));
366  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
367 
368  for (unsigned int i = 0; i < n; ++i)
369  for (unsigned int j = 0; j < m; ++j)
370  op.block(j, i).Tvmult_add(v.block(i), u.block(j));
371  };
372  }
373 
374 
389  template <typename PayloadBlockType>
390  class EmptyBlockPayload
391  {
392  public:
396  typedef PayloadBlockType BlockType;
397 
405  template <typename... Args>
406  EmptyBlockPayload (const Args &...)
407  { }
408  };
409 
410  } /*namespace BlockLinearOperator*/
411 } /*namespace internal*/
412 
413 
414 
419 
431 template <typename Range,
432  typename Domain,
433  typename BlockPayload,
434  typename BlockMatrixType>
436 block_operator(const BlockMatrixType &block_matrix)
437 {
439 
440  BlockLinearOperator<Range, Domain, BlockPayload> return_op (BlockPayload(block_matrix,block_matrix));
441 
442  return_op.n_block_rows = [&block_matrix]() -> unsigned int
443  {
444  return block_matrix.n_block_rows();
445  };
446 
447  return_op.n_block_cols = [&block_matrix]() -> unsigned int
448  {
449  return block_matrix.n_block_cols();
450  };
451 
452  return_op.block = [&block_matrix](unsigned int i, unsigned int j) -> BlockType
453  {
454 #ifdef DEBUG
455  const unsigned int m = block_matrix.n_block_rows();
456  const unsigned int n = block_matrix.n_block_cols();
457  Assert(i < m, ExcIndexRange (i, 0, m));
458  Assert(j < n, ExcIndexRange (j, 0, n));
459 #endif
460 
461  return BlockType(block_matrix.block(i, j));
462  };
463 
464  internal::BlockLinearOperator::populate_linear_operator_functions(return_op);
465  return return_op;
466 }
467 
468 
469 
497 template <size_t m, size_t n, typename Range, typename Domain, typename BlockPayload>
500 {
501  static_assert(m > 0 && n > 0,
502  "a blocked LinearOperator must consist of at least one block");
503 
505 
506  BlockLinearOperator<Range, Domain, BlockPayload> return_op ((BlockPayload())); // TODO: Create block payload so that this can be initialized correctly
507 
508  return_op.n_block_rows = []() -> unsigned int
509  {
510  return m;
511  };
512 
513  return_op.n_block_cols = []() -> unsigned int
514  {
515  return n;
516  };
517 
518  return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType
519  {
520  Assert(i < m, ExcIndexRange (i, 0, m));
521  Assert(j < n, ExcIndexRange (j, 0, n));
522 
523  return ops[i][j];
524  };
525 
526  internal::BlockLinearOperator::populate_linear_operator_functions(return_op);
527  return return_op;
528 }
529 
530 
531 
547 template <typename Range,
548  typename Domain,
549  typename BlockPayload,
550  typename BlockMatrixType>
552 block_diagonal_operator(const BlockMatrixType &block_matrix)
553 {
555 
556  BlockLinearOperator<Range, Domain, BlockPayload> return_op (BlockPayload(block_matrix,block_matrix));
557 
558  return_op.n_block_rows = [&block_matrix]() -> unsigned int
559  {
560  return block_matrix.n_block_rows();
561  };
562 
563  return_op.n_block_cols = [&block_matrix]() -> unsigned int
564  {
565  return block_matrix.n_block_cols();
566  };
567 
568  return_op.block = [&block_matrix](unsigned int i, unsigned int j) -> BlockType
569  {
570 #ifdef DEBUG
571  const unsigned int m = block_matrix.n_block_rows();
572  const unsigned int n = block_matrix.n_block_cols();
573  Assert(m == n, ExcDimensionMismatch(m, n));
574  Assert(i < m, ExcIndexRange (i, 0, m));
575  Assert(j < n, ExcIndexRange (j, 0, n));
576 #endif
577  if (i == j)
578  return BlockType(block_matrix.block(i, j));
579  else
580  return null_operator(BlockType(block_matrix.block(i, j)));
581  };
582 
583  internal::BlockLinearOperator::populate_linear_operator_functions(return_op);
584  return return_op;
585 }
586 
587 
588 
606 template <size_t m, typename Range, typename Domain, typename BlockPayload>
609 {
610  static_assert(m > 0,
611  "a blockdiagonal LinearOperator must consist of at least one block");
612 
614 
615  std::array<std::array<BlockType, m>, m> new_ops;
616 
617  // This is a bit tricky. We have to make sure that the off-diagonal
618  // elements of return_op.ops are populated correctly. They must be
619  // null_operators, but with correct reinit_domain_vector and
620  // reinit_range_vector functions.
621  for (unsigned int i = 0; i < m; ++i)
622  for (unsigned int j = 0; j < m; ++j)
623  if (i == j)
624  {
625  // diagonal elements are easy:
626  new_ops[i][j] = ops[i];
627  }
628  else
629  {
630  // create a null-operator...
631  new_ops[i][j] = null_operator(ops[i]);
632  // ... and fix up reinit_domain_vector:
633  new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
634  }
635 
636  return block_operator<m,m,Range,Domain>(new_ops);
637 }
638 
639 
640 
650 template <size_t m, typename Range, typename Domain, typename BlockPayload>
653 {
654  static_assert(m > 0,
655  "a blockdiagonal LinearOperator must consist of at least "
656  "one block");
657 
659  std::array<BlockType, m> new_ops;
660  new_ops.fill(op);
661 
662  return block_diagonal_operator(new_ops);
663 }
664 
665 
666 
668 
672 
708 template <typename Range, typename Domain, typename BlockPayload>
711  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
712 {
713  LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op ((typename BlockPayload::BlockType(diagonal_inverse)));
714 
715  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
716  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
717 
718  return_op.vmult = [block_operator, diagonal_inverse](Range &v, const Range &u)
719  {
720  const unsigned int m = block_operator.n_block_rows();
721  Assert(block_operator.n_block_cols() == m,
722  ExcDimensionMismatch(block_operator.n_block_cols(), m));
723  Assert(diagonal_inverse.n_block_rows() == m,
724  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
725  Assert(diagonal_inverse.n_block_cols() == m,
726  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
727  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
728  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
729 
730  if (m == 0)
731  return;
732 
733  diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
734  for (unsigned int i = 1; i < m; ++i)
735  {
736  auto &dst = v.block(i);
737  dst = u.block(i);
738  dst *= -1.;
739  for (unsigned int j = 0; j < i; ++j)
740  block_operator.block(i, j).vmult_add(dst, v.block(j));
741  dst *= -1.;
742  diagonal_inverse.block(i, i).vmult(dst, dst); // uses intermediate storage
743  }
744  };
745 
746  return_op.vmult_add = [block_operator, diagonal_inverse](Range &v, const Range &u)
747  {
748  const unsigned int m = block_operator.n_block_rows();
749  Assert(block_operator.n_block_cols() == m,
750  ExcDimensionMismatch(block_operator.n_block_cols(), m));
751  Assert(diagonal_inverse.n_block_rows() == m,
752  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
753  Assert(diagonal_inverse.n_block_cols() == m,
754  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
755  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
756  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
757 
758  if (m == 0)
759  return;
760 
762  typename VectorMemory<typename Range::BlockType>::Pointer tmp (vector_memory);
763 
764  diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0));
765 
766  for (unsigned int i = 1; i < m; ++i)
767  {
768  diagonal_inverse.block(i, i).reinit_range_vector(*tmp, /*bool omit_zeroing_entries=*/ true);
769  *tmp = u.block(i);
770  *tmp *= -1.;
771  for (unsigned int j = 0; j < i; ++j)
772  block_operator.block(i, j).vmult_add(*tmp, v.block(j));
773  *tmp *= -1.;
774  diagonal_inverse.block(i, i).vmult_add(v.block(i),*tmp);
775  }
776  };
777 
778  return return_op;
779 }
780 
781 
782 
818 template <typename Range, typename Domain, typename BlockPayload>
821  const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
822 {
823  LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op ((typename BlockPayload::BlockType(diagonal_inverse)));
824 
825  return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
826  return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
827 
828  return_op.vmult = [block_operator, diagonal_inverse](Range &v, const Range &u)
829  {
830  const unsigned int m = block_operator.n_block_rows();
831  Assert(block_operator.n_block_cols() == m,
832  ExcDimensionMismatch(block_operator.n_block_cols(), m));
833  Assert(diagonal_inverse.n_block_rows() == m,
834  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
835  Assert(diagonal_inverse.n_block_cols() == m,
836  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
837  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
838  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
839 
840  if (m == 0)
841  return;
842 
843  diagonal_inverse.block(m-1, m-1).vmult(v.block(m-1),u.block(m-1));
844 
845  for (int i = m - 2; i >= 0; --i)
846  {
847  auto &dst = v.block(i);
848  dst = u.block(i);
849  dst *= -1.;
850  for (unsigned int j = i + 1; j < m; ++j)
851  block_operator.block(i, j).vmult_add(dst, v.block(j));
852  dst *= -1.;
853  diagonal_inverse.block(i, i).vmult(dst, dst); // uses intermediate storage
854  }
855  };
856 
857  return_op.vmult_add = [block_operator, diagonal_inverse](Range &v, const Range &u)
858  {
859  const unsigned int m = block_operator.n_block_rows();
860  Assert(block_operator.n_block_cols() == m,
861  ExcDimensionMismatch(block_operator.n_block_cols(), m));
862  Assert(diagonal_inverse.n_block_rows() == m,
863  ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m));
864  Assert(diagonal_inverse.n_block_cols() == m,
865  ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m));
866  Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m));
867  Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m));
869  typename VectorMemory<typename Range::BlockType>::Pointer tmp (vector_memory);
870 
871  if (m == 0)
872  return;
873 
874  diagonal_inverse.block(m-1, m-1).vmult_add(v.block(m-1),u.block(m-1));
875 
876  for (int i = m - 2; i >= 0; --i)
877  {
878  diagonal_inverse.block(i, i).reinit_range_vector(*tmp, /*bool omit_zeroing_entries=*/ true);
879  *tmp = u.block(i);
880  *tmp *= -1.;
881  for (unsigned int j = i + 1; j < m; ++j)
882  block_operator.block(i, j).vmult_add(*tmp,v.block(j));
883  *tmp *= -1.;
884  diagonal_inverse.block(i, i).vmult_add(v.block(i),*tmp);
885  }
886  };
887 
888  return return_op;
889 }
890 
892 
893 DEAL_II_NAMESPACE_CLOSE
894 
895 #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)