Reference documentation for deal.II version Git bed997f895 2020-09-22 11:49:20 -0400
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
matrix_block.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2007 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_matrix_block_h
17 #define dealii_matrix_block_h
18 
19 #include <deal.II/base/config.h>
20 
22 
26 
31 
32 #include <memory>
33 
34 
36 
37 // Forward declarations
38 #ifndef DOXYGEN
39 template <typename MatrixType>
40 class MatrixBlock;
41 #endif
42 
43 namespace internal
44 {
45  template <typename MatrixType>
46  void
48 
49  template <typename number>
50  void
52  const BlockSparsityPattern & p);
53 } // namespace internal
54 
110 template <typename MatrixType>
111 class MatrixBlock : public Subscriptor
112 {
113 public:
118 
122  using value_type = typename MatrixType::value_type;
123 
127  MatrixBlock();
128 
132  MatrixBlock(const MatrixBlock<MatrixType> &M) = default;
133 
138  operator=(const MatrixBlock<MatrixType> &) = default;
139 
145 
153  void
154  reinit(const BlockSparsityPattern &sparsity);
155 
156  operator MatrixType &();
157  operator const MatrixType &() const;
158 
163  void
164  add(const size_type i,
165  const size_type j,
166  const typename MatrixType::value_type value);
167 
183  template <typename number>
184  void
185  add(const std::vector<size_type> &indices,
186  const FullMatrix<number> & full_matrix,
187  const bool elide_zero_values = true);
188 
203  template <typename number>
204  void
205  add(const std::vector<size_type> &row_indices,
206  const std::vector<size_type> &col_indices,
207  const FullMatrix<number> & full_matrix,
208  const bool elide_zero_values = true);
209 
226  template <typename number>
227  void
228  add(const size_type row_index,
229  const std::vector<size_type> &col_indices,
230  const std::vector<number> & values,
231  const bool elide_zero_values = true);
232 
242  template <typename number>
243  void
244  add(const size_type row,
245  const size_type n_cols,
246  const size_type *col_indices,
247  const number * values,
248  const bool elide_zero_values = true,
249  const bool col_indices_are_sorted = false);
250 
256  template <class VectorType>
257  void
258  vmult(VectorType &w, const VectorType &v) const;
259 
265  template <class VectorType>
266  void
267  vmult_add(VectorType &w, const VectorType &v) const;
268 
274  template <class VectorType>
275  void
276  Tvmult(VectorType &w, const VectorType &v) const;
277 
283  template <class VectorType>
284  void
285  Tvmult_add(VectorType &w, const VectorType &v) const;
286 
290  std::size_t
291  memory_consumption() const;
292 
297  DeclException2(ExcBlockIndexMismatch,
298  size_type,
299  size_type,
300  << "Block index " << arg1 << " does not match " << arg2);
301 
312 
316  MatrixType matrix;
317 
318 private:
330 
331  template <class OTHER_MatrixType>
332  friend void
334  const BlockSparsityPattern &);
335 
336  template <typename number>
337  friend void
339  const BlockSparsityPattern & p);
340 };
341 
342 
351 template <typename MatrixType>
352 class MatrixBlockVector : private AnyData
353 {
354 public:
359 
364 
369  using ptr_type = std::shared_ptr<value_type>;
370 
375  void
376  add(size_type row, size_type column, const std::string &name);
377 
382  void
383  reinit(const BlockSparsityPattern &sparsity);
384 
395  void
396  clear(bool really_clean = false);
397 
401  std::size_t
402  memory_consumption() const;
403 
407  const value_type &
408  block(size_type i) const;
409 
413  value_type &
414  block(size_type i);
415 
419  MatrixType &
420  matrix(size_type i);
421 
425  using AnyData::name;
426  using AnyData::size;
427  using AnyData::subscribe;
428  using AnyData::unsubscribe;
429 };
430 
431 
440 template <typename MatrixType>
442 {
443 public:
448 
462  MGMatrixBlockVector(const bool edge_matrices = false,
463  const bool edge_flux_matrices = false);
464 
468  unsigned int
469  size() const;
470 
476  void
477  add(size_type row, size_type column, const std::string &name);
478 
485  void
486  reinit_matrix(const MGLevelObject<BlockSparsityPattern> &sparsity);
494  void
495  reinit_edge(const MGLevelObject<BlockSparsityPattern> &sparsity);
502  void
503  reinit_edge_flux(const MGLevelObject<BlockSparsityPattern> &sparsity);
504 
515  void
516  clear(bool really_clean = false);
517 
521  const value_type &
522  block(size_type i) const;
523 
527  value_type &
528  block(size_type i);
529 
534  const value_type &
535  block_in(size_type i) const;
536 
540  value_type &
541  block_in(size_type i);
542 
547  const value_type &
548  block_out(size_type i) const;
549 
553  value_type &
554  block_out(size_type i);
555 
560  const value_type &
561  block_up(size_type i) const;
562 
566  value_type &
567  block_up(size_type i);
568 
573  const value_type &
574  block_down(size_type i) const;
575 
579  value_type &
580  block_down(size_type i);
581 
585  std::size_t
586  memory_consumption() const;
587 
588 private:
590  void
591  clear_object(AnyData &);
592 
594  const bool edge_matrices;
595 
597  const bool edge_flux_matrices;
598 
609 };
610 
611 
612 //----------------------------------------------------------------------//
613 
614 namespace internal
615 {
616  template <typename MatrixType>
617  void
619  {
622  }
623 
624 
625  template <typename number>
626  void
628  const BlockSparsityPattern & p)
629  {
630  v.row_indices = p.get_row_indices();
631  v.column_indices = p.get_column_indices();
632  v.matrix.reinit(p.block(v.row, v.column));
633  }
634 } // namespace internal
635 
636 
637 template <typename MatrixType>
639  : row(numbers::invalid_size_type)
640  , column(numbers::invalid_size_type)
641 {}
642 
643 
644 template <typename MatrixType>
646  : row(i)
647  , column(j)
648 {}
649 
650 
651 template <typename MatrixType>
652 inline void
654 {
655  internal::reinit(*this, sparsity);
656 }
657 
658 
659 template <typename MatrixType>
661 {
662  return matrix;
663 }
664 
665 
666 template <typename MatrixType>
667 inline MatrixBlock<MatrixType>::operator const MatrixType &() const
668 {
669  return matrix;
670 }
671 
672 
673 template <typename MatrixType>
674 inline void
676  const size_type gj,
677  const typename MatrixType::value_type value)
678 {
681 
682  const std::pair<unsigned int, size_type> bi = row_indices.global_to_local(gi);
683  const std::pair<unsigned int, size_type> bj =
685 
686  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
687  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
688 
689  matrix.add(bi.second, bj.second, value);
690 }
691 
692 
693 template <typename MatrixType>
694 template <typename number>
695 inline void
696 MatrixBlock<MatrixType>::add(const std::vector<size_type> &r_indices,
697  const std::vector<size_type> &c_indices,
698  const FullMatrix<number> & values,
699  const bool elide_zero_values)
700 {
703 
704  AssertDimension(r_indices.size(), values.m());
705  AssertDimension(c_indices.size(), values.n());
706 
707  for (size_type i = 0; i < row_indices.size(); ++i)
708  add(r_indices[i],
709  c_indices.size(),
710  c_indices.data(),
711  &values(i, 0),
712  elide_zero_values);
713 }
714 
715 
716 template <typename MatrixType>
717 template <typename number>
718 inline void
720  const size_type n_cols,
721  const size_type *col_indices,
722  const number * values,
723  const bool,
724  const bool)
725 {
728 
729  const std::pair<unsigned int, size_type> bi =
731 
732  // In debug mode, we check whether
733  // all indices are in the correct
734  // block.
735 
736  // Actually, for the time being, we
737  // leave it at this. While it may
738  // not be the most efficient way,
739  // it is at least thread safe.
740  //#ifdef DEBUG
741  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
742 
743  for (size_type j = 0; j < n_cols; ++j)
744  {
745  const std::pair<unsigned int, size_type> bj =
746  column_indices.global_to_local(col_indices[j]);
747  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
748 
749  matrix.add(bi.second, bj.second, values[j]);
750  }
751  //#endif
752 }
753 
754 
755 template <typename MatrixType>
756 template <typename number>
757 inline void
758 MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
759  const FullMatrix<number> & values,
760  const bool elide_zero_values)
761 {
764 
765  AssertDimension(indices.size(), values.m());
766  Assert(values.n() == values.m(), ExcNotQuadratic());
767 
768  for (size_type i = 0; i < indices.size(); ++i)
769  add(indices[i],
770  indices.size(),
771  indices.data(),
772  &values(i, 0),
773  elide_zero_values);
774 }
775 
776 
777 
778 template <typename MatrixType>
779 template <typename number>
780 inline void
782  const std::vector<size_type> &col_indices,
783  const std::vector<number> & values,
784  const bool elide_zero_values)
785 {
788 
789  AssertDimension(col_indices.size(), values.size());
790  add(row,
791  col_indices.size(),
792  col_indices.data(),
793  values.data(),
794  elide_zero_values);
795 }
796 
797 
798 template <typename MatrixType>
799 template <class VectorType>
800 inline void
802 {
803  matrix.vmult(w, v);
804 }
805 
806 
807 template <typename MatrixType>
808 template <class VectorType>
809 inline void
811 {
812  matrix.vmult_add(w, v);
813 }
814 
815 
816 template <typename MatrixType>
817 template <class VectorType>
818 inline void
820 {
821  matrix.Tvmult(w, v);
822 }
823 
824 
825 template <typename MatrixType>
826 template <class VectorType>
827 inline void
829 {
830  matrix.Tvmult_add(w, v);
831 }
832 
833 
834 template <typename MatrixType>
835 inline std::size_t
837 {
838  return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
839  sizeof(matrix));
840 }
841 
842 //----------------------------------------------------------------------//
843 
844 template <typename MatrixType>
845 inline void
848  const std::string &name)
849 {
850  ptr_type p(new value_type(row, column));
851  AnyData::add(p, name);
852 }
853 
854 
855 template <typename MatrixType>
856 inline void
858 {
859  for (size_type i = 0; i < this->size(); ++i)
860  {
861  block(i).reinit(sparsity);
862  }
863 }
864 
865 
866 template <typename MatrixType>
867 inline void
869 {
870  if (really_clean)
871  {
872  Assert(false, ExcNotImplemented());
873  }
874  else
875  {
876  for (size_type i = 0; i < this->size(); ++i)
877  matrix(i).clear();
878  }
879 }
880 
881 
882 
883 template <typename MatrixType>
884 inline const MatrixBlock<MatrixType> &
886 {
887  return *this->read<ptr_type>(i);
888 }
889 
890 
891 template <typename MatrixType>
894 {
895  return *this->entry<ptr_type>(i);
896 }
897 
898 
899 template <typename MatrixType>
900 inline MatrixType &
902 {
903  return this->entry<ptr_type>(i)->matrix;
904 }
905 
906 
907 
908 //----------------------------------------------------------------------//
909 
910 template <typename MatrixType>
912  const bool f)
913  : edge_matrices(e)
914  , edge_flux_matrices(f)
915 {}
916 
917 
918 template <typename MatrixType>
919 inline unsigned int
921 {
922  return matrices.size();
923 }
924 
925 
926 template <typename MatrixType>
927 inline void
929  size_type column,
930  const std::string &name)
931 {
933  p[0].row = row;
934  p[0].column = column;
935 
936  matrices.add(p, name);
937  if (edge_matrices)
938  {
939  matrices_in.add(p, name);
940  matrices_out.add(p, name);
941  }
942  if (edge_flux_matrices)
943  {
944  flux_matrices_up.add(p, name);
945  flux_matrices_down.add(p, name);
946  }
947 }
948 
949 
950 template <typename MatrixType>
953 {
954  return *matrices.read<const MGLevelObject<MatrixType> *>(i);
955 }
956 
957 
958 template <typename MatrixType>
961 {
963 }
964 
965 
966 template <typename MatrixType>
969 {
970  return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
971 }
972 
973 
974 template <typename MatrixType>
977 {
979 }
980 
981 
982 template <typename MatrixType>
985 {
986  return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
987 }
988 
989 
990 template <typename MatrixType>
993 {
995 }
996 
997 
998 template <typename MatrixType>
1001 {
1002  return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
1003 }
1004 
1005 
1006 template <typename MatrixType>
1009 {
1011 }
1012 
1013 
1014 template <typename MatrixType>
1017 {
1019 }
1020 
1021 
1022 template <typename MatrixType>
1025 {
1027 }
1028 
1029 
1030 template <typename MatrixType>
1031 inline void
1033  const MGLevelObject<BlockSparsityPattern> &sparsity)
1034 {
1035  for (size_type i = 0; i < this->size(); ++i)
1036  {
1038  const size_type row = o[o.min_level()].row;
1039  const size_type col = o[o.min_level()].column;
1040 
1041  o.resize(sparsity.min_level(), sparsity.max_level());
1042  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1043  {
1044  o[level].row = row;
1045  o[level].column = col;
1046  internal::reinit(o[level], sparsity[level]);
1047  }
1048  }
1049 }
1050 
1051 
1052 template <typename MatrixType>
1053 inline void
1055  const MGLevelObject<BlockSparsityPattern> &sparsity)
1056 {
1057  for (size_type i = 0; i < this->size(); ++i)
1058  {
1060  const size_type row = o[o.min_level()].row;
1061  const size_type col = o[o.min_level()].column;
1062 
1063  block_in(i).resize(sparsity.min_level(), sparsity.max_level());
1064  block_out(i).resize(sparsity.min_level(), sparsity.max_level());
1065  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1066  {
1067  block_in(i)[level].row = row;
1068  block_in(i)[level].column = col;
1069  internal::reinit(block_in(i)[level], sparsity[level]);
1070  block_out(i)[level].row = row;
1071  block_out(i)[level].column = col;
1072  internal::reinit(block_out(i)[level], sparsity[level]);
1073  }
1074  }
1075 }
1076 
1077 
1078 template <typename MatrixType>
1079 inline void
1081  const MGLevelObject<BlockSparsityPattern> &sparsity)
1082 {
1083  for (size_type i = 0; i < this->size(); ++i)
1084  {
1086  const size_type row = o[o.min_level()].row;
1087  const size_type col = o[o.min_level()].column;
1088 
1089  block_up(i).resize(sparsity.min_level(), sparsity.max_level());
1090  block_down(i).resize(sparsity.min_level(), sparsity.max_level());
1091  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1092  {
1093  block_up(i)[level].row = row;
1094  block_up(i)[level].column = col;
1095  internal::reinit(block_up(i)[level], sparsity[level]);
1096  block_down(i)[level].row = row;
1097  block_down(i)[level].column = col;
1098  internal::reinit(block_down(i)[level], sparsity[level]);
1099  }
1100  }
1101 }
1102 
1103 
1104 template <typename MatrixType>
1105 inline void
1107 {
1108  for (size_type i = 0; i < mo.size(); ++i)
1109  {
1112  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1113  o[level].matrix.clear();
1114  }
1115 }
1116 
1117 
1118 template <typename MatrixType>
1119 inline void
1121 {
1122  if (really_clean)
1123  {
1124  Assert(false, ExcNotImplemented());
1125  }
1126  else
1127  {
1133  }
1134 }
1135 
1136 
1137 
1139 
1140 #endif
void reinit_matrix(const MGLevelObject< BlockSparsityPattern > &sparsity)
unsigned int max_level() const
size_type m() const
const types::global_dof_index invalid_size_type
Definition: types.h:205
const bool edge_matrices
Flag for storing matrices_in and matrices_out.
Definition: matrix_block.h:594
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618
type entry(const std::string &name)
Access to stored data object by name.
Definition: any_data.h:351
AnyData matrices
The level matrices.
Definition: matrix_block.h:600
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:538
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1568
Contents is actually a matrix.
AnyData matrices_out
The matrix from the refinement edge to the interior of a level.
Definition: matrix_block.h:604
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
const type read(const std::string &name) const
Dedicated read only access by name.
Definition: any_data.h:374
MatrixType & matrix(size_type i)
Definition: matrix_block.h:901
const value_type & block(size_type i) const
Definition: matrix_block.h:885
void Tvmult(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:819
void vmult(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:801
const value_type & block(size_type i) const
Definition: matrix_block.h:952
void unsubscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
Definition: subscriptor.cc:156
AnyData matrices_in
The matrix from the interior of a level to the refinement edge.
Definition: matrix_block.h:602
static ::ExceptionBase & ExcNotInitialized()
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
Definition: matrix_block.h:675
size_type column
Definition: matrix_block.h:311
unsigned int min_level() const
SparsityPatternType & block(const size_type row, const size_type column)
const bool edge_flux_matrices
Flag for storing flux_matrices_up and flux_matrices_down.
Definition: matrix_block.h:597
void clear(bool really_clean=false)
Definition: matrix_block.h:868
const value_type & block_out(size_type i) const
Definition: matrix_block.h:984
void add(size_type row, size_type column, const std::string &name)
Definition: matrix_block.h:846
void reinit_edge(const MGLevelObject< BlockSparsityPattern > &sparsity)
MGMatrixBlockVector(const bool edge_matrices=false, const bool edge_flux_matrices=false)
Definition: matrix_block.h:911
size_type n() const
void subscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
Definition: subscriptor.cc:136
#define Assert(cond, exc)
Definition: exceptions.h:1411
AnyData flux_matrices_down
The DG flux from a level to the lower level.
Definition: matrix_block.h:606
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
void reinit(MatrixBlock<::SparseMatrix< number >> &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:627
const BlockIndices & get_row_indices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
unsigned int level
Definition: grid_out.cc:4338
std::size_t memory_consumption() const
Definition: matrix_block.h:836
BlockIndices column_indices
Definition: matrix_block.h:329
const value_type & block_in(size_type i) const
Definition: matrix_block.h:968
void vmult_add(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:810
std::shared_ptr< value_type > ptr_type
Definition: matrix_block.h:369
BlockIndices row_indices
Definition: matrix_block.h:323
typename FullMatrix< number > ::value_type value_type
Definition: matrix_block.h:122
MatrixType matrix
Definition: matrix_block.h:316
static ::ExceptionBase & ExcNotQuadratic()
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int global_dof_index
Definition: types.h:76
void add(size_type row, size_type column, const std::string &name)
Definition: matrix_block.h:928
static ::ExceptionBase & ExcBlockIndexMismatch(size_type arg1, size_type arg2)
void clear(bool really_clean=false)
unsigned int size() const
Definition: matrix_block.h:920
void reinit_edge_flux(const MGLevelObject< BlockSparsityPattern > &sparsity)
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:857
const value_type & block_up(size_type i) const
void add(type entry, const std::string &name)
Add a new data object.
Definition: any_data.h:432
unsigned int size() const
Number of stored data objects.
Definition: any_data.h:223
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:653
AnyData flux_matrices_up
The DG flux from the lower level to a level.
Definition: matrix_block.h:608
const value_type & block_down(size_type i) const
static ::ExceptionBase & ExcNotImplemented()
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel)
void clear_object(AnyData &)
Clear one of the matrix objects.
unsigned int size() const
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
const BlockIndices & get_column_indices() const
void Tvmult_add(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:828
const std::string & name(const unsigned int i) const
Name of object at index.
Definition: any_data.h:311
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
size_type row
Definition: matrix_block.h:306