Reference documentation for deal.II version 9.2.0
\(\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 
111 template <typename MatrixType>
112 class MatrixBlock : public Subscriptor
113 {
114 public:
119 
123  using value_type = typename MatrixType::value_type;
124 
128  MatrixBlock();
129 
133  MatrixBlock(const MatrixBlock<MatrixType> &M) = default;
134 
139  operator=(const MatrixBlock<MatrixType> &) = default;
140 
146 
154  void
155  reinit(const BlockSparsityPattern &sparsity);
156 
157  operator MatrixType &();
158  operator const MatrixType &() const;
159 
164  void
165  add(const size_type i,
166  const size_type j,
167  const typename MatrixType::value_type value);
168 
184  template <typename number>
185  void
186  add(const std::vector<size_type> &indices,
187  const FullMatrix<number> & full_matrix,
188  const bool elide_zero_values = true);
189 
204  template <typename number>
205  void
206  add(const std::vector<size_type> &row_indices,
207  const std::vector<size_type> &col_indices,
208  const FullMatrix<number> & full_matrix,
209  const bool elide_zero_values = true);
210 
227  template <typename number>
228  void
229  add(const size_type row_index,
230  const std::vector<size_type> &col_indices,
231  const std::vector<number> & values,
232  const bool elide_zero_values = true);
233 
243  template <typename number>
244  void
245  add(const size_type row,
246  const size_type n_cols,
247  const size_type *col_indices,
248  const number * values,
249  const bool elide_zero_values = true,
250  const bool col_indices_are_sorted = false);
251 
257  template <class VectorType>
258  void
259  vmult(VectorType &w, const VectorType &v) const;
260 
266  template <class VectorType>
267  void
268  vmult_add(VectorType &w, const VectorType &v) const;
269 
275  template <class VectorType>
276  void
277  Tvmult(VectorType &w, const VectorType &v) const;
278 
284  template <class VectorType>
285  void
286  Tvmult_add(VectorType &w, const VectorType &v) const;
287 
291  std::size_t
292  memory_consumption() const;
293 
298  DeclException2(ExcBlockIndexMismatch,
299  size_type,
300  size_type,
301  << "Block index " << arg1 << " does not match " << arg2);
302 
313 
317  MatrixType matrix;
318 
319 private:
331 
332  template <class OTHER_MatrixType>
333  friend void
335  const BlockSparsityPattern &);
336 
337  template <typename number>
338  friend void
340  const BlockSparsityPattern & p);
341 };
342 
343 
353 template <typename MatrixType>
354 class MatrixBlockVector : private AnyData
355 {
356 public:
361 
366 
371  using ptr_type = std::shared_ptr<value_type>;
372 
377  void
378  add(size_type row, size_type column, const std::string &name);
379 
384  void
385  reinit(const BlockSparsityPattern &sparsity);
386 
397  void
398  clear(bool really_clean = false);
399 
403  std::size_t
404  memory_consumption() const;
405 
409  const value_type &
410  block(size_type i) const;
411 
415  value_type &
416  block(size_type i);
417 
421  MatrixType &
422  matrix(size_type i);
423 
427  using AnyData::name;
428  using AnyData::size;
429  using AnyData::subscribe;
430  using AnyData::unsubscribe;
431 };
432 
433 
443 template <typename MatrixType>
445 {
446 public:
451 
465  MGMatrixBlockVector(const bool edge_matrices = false,
466  const bool edge_flux_matrices = false);
467 
471  unsigned int
472  size() const;
473 
479  void
480  add(size_type row, size_type column, const std::string &name);
481 
488  void
489  reinit_matrix(const MGLevelObject<BlockSparsityPattern> &sparsity);
497  void
498  reinit_edge(const MGLevelObject<BlockSparsityPattern> &sparsity);
505  void
506  reinit_edge_flux(const MGLevelObject<BlockSparsityPattern> &sparsity);
507 
518  void
519  clear(bool really_clean = false);
520 
524  const value_type &
525  block(size_type i) const;
526 
530  value_type &
531  block(size_type i);
532 
537  const value_type &
538  block_in(size_type i) const;
539 
543  value_type &
544  block_in(size_type i);
545 
550  const value_type &
551  block_out(size_type i) const;
552 
556  value_type &
557  block_out(size_type i);
558 
563  const value_type &
564  block_up(size_type i) const;
565 
569  value_type &
570  block_up(size_type i);
571 
576  const value_type &
577  block_down(size_type i) const;
578 
582  value_type &
583  block_down(size_type i);
584 
588  std::size_t
589  memory_consumption() const;
590 
591 private:
593  void
594  clear_object(AnyData &);
595 
597  const bool edge_matrices;
598 
600  const bool edge_flux_matrices;
601 
612 };
613 
614 
615 //----------------------------------------------------------------------//
616 
617 namespace internal
618 {
619  template <typename MatrixType>
620  void
622  {
625  }
626 
627 
628  template <typename number>
629  void
631  const BlockSparsityPattern & p)
632  {
633  v.row_indices = p.get_row_indices();
634  v.column_indices = p.get_column_indices();
635  v.matrix.reinit(p.block(v.row, v.column));
636  }
637 } // namespace internal
638 
639 
640 template <typename MatrixType>
642  : row(numbers::invalid_size_type)
643  , column(numbers::invalid_size_type)
644 {}
645 
646 
647 template <typename MatrixType>
649  : row(i)
650  , column(j)
651 {}
652 
653 
654 template <typename MatrixType>
655 inline void
657 {
658  internal::reinit(*this, sparsity);
659 }
660 
661 
662 template <typename MatrixType>
664 {
665  return matrix;
666 }
667 
668 
669 template <typename MatrixType>
670 inline MatrixBlock<MatrixType>::operator const MatrixType &() const
671 {
672  return matrix;
673 }
674 
675 
676 template <typename MatrixType>
677 inline void
679  const size_type gj,
680  const typename MatrixType::value_type value)
681 {
684 
685  const std::pair<unsigned int, size_type> bi = row_indices.global_to_local(gi);
686  const std::pair<unsigned int, size_type> bj =
688 
689  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
690  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
691 
692  matrix.add(bi.second, bj.second, value);
693 }
694 
695 
696 template <typename MatrixType>
697 template <typename number>
698 inline void
699 MatrixBlock<MatrixType>::add(const std::vector<size_type> &r_indices,
700  const std::vector<size_type> &c_indices,
701  const FullMatrix<number> & values,
702  const bool elide_zero_values)
703 {
706 
707  AssertDimension(r_indices.size(), values.m());
708  AssertDimension(c_indices.size(), values.n());
709 
710  for (size_type i = 0; i < row_indices.size(); ++i)
711  add(r_indices[i],
712  c_indices.size(),
713  c_indices.data(),
714  &values(i, 0),
715  elide_zero_values);
716 }
717 
718 
719 template <typename MatrixType>
720 template <typename number>
721 inline void
723  const size_type n_cols,
724  const size_type *col_indices,
725  const number * values,
726  const bool,
727  const bool)
728 {
731 
732  const std::pair<unsigned int, size_type> bi =
734 
735  // In debug mode, we check whether
736  // all indices are in the correct
737  // block.
738 
739  // Actually, for the time being, we
740  // leave it at this. While it may
741  // not be the most efficient way,
742  // it is at least thread safe.
743  //#ifdef DEBUG
744  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
745 
746  for (size_type j = 0; j < n_cols; ++j)
747  {
748  const std::pair<unsigned int, size_type> bj =
749  column_indices.global_to_local(col_indices[j]);
750  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
751 
752  matrix.add(bi.second, bj.second, values[j]);
753  }
754  //#endif
755 }
756 
757 
758 template <typename MatrixType>
759 template <typename number>
760 inline void
761 MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
762  const FullMatrix<number> & values,
763  const bool elide_zero_values)
764 {
767 
768  AssertDimension(indices.size(), values.m());
769  Assert(values.n() == values.m(), ExcNotQuadratic());
770 
771  for (size_type i = 0; i < indices.size(); ++i)
772  add(indices[i],
773  indices.size(),
774  indices.data(),
775  &values(i, 0),
776  elide_zero_values);
777 }
778 
779 
780 
781 template <typename MatrixType>
782 template <typename number>
783 inline void
785  const std::vector<size_type> &col_indices,
786  const std::vector<number> & values,
787  const bool elide_zero_values)
788 {
791 
792  AssertDimension(col_indices.size(), values.size());
793  add(row,
794  col_indices.size(),
795  col_indices.data(),
796  values.data(),
797  elide_zero_values);
798 }
799 
800 
801 template <typename MatrixType>
802 template <class VectorType>
803 inline void
805 {
806  matrix.vmult(w, v);
807 }
808 
809 
810 template <typename MatrixType>
811 template <class VectorType>
812 inline void
814 {
815  matrix.vmult_add(w, v);
816 }
817 
818 
819 template <typename MatrixType>
820 template <class VectorType>
821 inline void
823 {
824  matrix.Tvmult(w, v);
825 }
826 
827 
828 template <typename MatrixType>
829 template <class VectorType>
830 inline void
832 {
833  matrix.Tvmult_add(w, v);
834 }
835 
836 
837 template <typename MatrixType>
838 inline std::size_t
840 {
841  return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
842  sizeof(matrix));
843 }
844 
845 //----------------------------------------------------------------------//
846 
847 template <typename MatrixType>
848 inline void
851  const std::string &name)
852 {
853  ptr_type p(new value_type(row, column));
854  AnyData::add(p, name);
855 }
856 
857 
858 template <typename MatrixType>
859 inline void
861 {
862  for (size_type i = 0; i < this->size(); ++i)
863  {
864  block(i).reinit(sparsity);
865  }
866 }
867 
868 
869 template <typename MatrixType>
870 inline void
872 {
873  if (really_clean)
874  {
875  Assert(false, ExcNotImplemented());
876  }
877  else
878  {
879  for (size_type i = 0; i < this->size(); ++i)
880  matrix(i).clear();
881  }
882 }
883 
884 
885 
886 template <typename MatrixType>
887 inline const MatrixBlock<MatrixType> &
889 {
890  return *this->read<ptr_type>(i);
891 }
892 
893 
894 template <typename MatrixType>
897 {
898  return *this->entry<ptr_type>(i);
899 }
900 
901 
902 template <typename MatrixType>
903 inline MatrixType &
905 {
906  return this->entry<ptr_type>(i)->matrix;
907 }
908 
909 
910 
911 //----------------------------------------------------------------------//
912 
913 template <typename MatrixType>
915  const bool f)
916  : edge_matrices(e)
917  , edge_flux_matrices(f)
918 {}
919 
920 
921 template <typename MatrixType>
922 inline unsigned int
924 {
925  return matrices.size();
926 }
927 
928 
929 template <typename MatrixType>
930 inline void
932  size_type column,
933  const std::string &name)
934 {
936  p[0].row = row;
937  p[0].column = column;
938 
939  matrices.add(p, name);
940  if (edge_matrices)
941  {
942  matrices_in.add(p, name);
943  matrices_out.add(p, name);
944  }
945  if (edge_flux_matrices)
946  {
947  flux_matrices_up.add(p, name);
948  flux_matrices_down.add(p, name);
949  }
950 }
951 
952 
953 template <typename MatrixType>
956 {
957  return *matrices.read<const MGLevelObject<MatrixType> *>(i);
958 }
959 
960 
961 template <typename MatrixType>
964 {
966 }
967 
968 
969 template <typename MatrixType>
972 {
973  return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
974 }
975 
976 
977 template <typename MatrixType>
980 {
982 }
983 
984 
985 template <typename MatrixType>
988 {
989  return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
990 }
991 
992 
993 template <typename MatrixType>
996 {
998 }
999 
1000 
1001 template <typename MatrixType>
1004 {
1005  return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
1006 }
1007 
1008 
1009 template <typename MatrixType>
1012 {
1014 }
1015 
1016 
1017 template <typename MatrixType>
1020 {
1022 }
1023 
1024 
1025 template <typename MatrixType>
1028 {
1030 }
1031 
1032 
1033 template <typename MatrixType>
1034 inline void
1036  const MGLevelObject<BlockSparsityPattern> &sparsity)
1037 {
1038  for (size_type i = 0; i < this->size(); ++i)
1039  {
1041  const size_type row = o[o.min_level()].row;
1042  const size_type col = o[o.min_level()].column;
1043 
1044  o.resize(sparsity.min_level(), sparsity.max_level());
1045  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1046  {
1047  o[level].row = row;
1048  o[level].column = col;
1049  internal::reinit(o[level], sparsity[level]);
1050  }
1051  }
1052 }
1053 
1054 
1055 template <typename MatrixType>
1056 inline void
1058  const MGLevelObject<BlockSparsityPattern> &sparsity)
1059 {
1060  for (size_type i = 0; i < this->size(); ++i)
1061  {
1063  const size_type row = o[o.min_level()].row;
1064  const size_type col = o[o.min_level()].column;
1065 
1066  block_in(i).resize(sparsity.min_level(), sparsity.max_level());
1067  block_out(i).resize(sparsity.min_level(), sparsity.max_level());
1068  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1069  {
1070  block_in(i)[level].row = row;
1071  block_in(i)[level].column = col;
1072  internal::reinit(block_in(i)[level], sparsity[level]);
1073  block_out(i)[level].row = row;
1074  block_out(i)[level].column = col;
1075  internal::reinit(block_out(i)[level], sparsity[level]);
1076  }
1077  }
1078 }
1079 
1080 
1081 template <typename MatrixType>
1082 inline void
1084  const MGLevelObject<BlockSparsityPattern> &sparsity)
1085 {
1086  for (size_type i = 0; i < this->size(); ++i)
1087  {
1089  const size_type row = o[o.min_level()].row;
1090  const size_type col = o[o.min_level()].column;
1091 
1092  block_up(i).resize(sparsity.min_level(), sparsity.max_level());
1093  block_down(i).resize(sparsity.min_level(), sparsity.max_level());
1094  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1095  {
1096  block_up(i)[level].row = row;
1097  block_up(i)[level].column = col;
1098  internal::reinit(block_up(i)[level], sparsity[level]);
1099  block_down(i)[level].row = row;
1100  block_down(i)[level].column = col;
1101  internal::reinit(block_down(i)[level], sparsity[level]);
1102  }
1103  }
1104 }
1105 
1106 
1107 template <typename MatrixType>
1108 inline void
1110 {
1111  for (size_type i = 0; i < mo.size(); ++i)
1112  {
1115  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1116  o[level].matrix.clear();
1117  }
1118 }
1119 
1120 
1121 template <typename MatrixType>
1122 inline void
1124 {
1125  if (really_clean)
1126  {
1127  Assert(false, ExcNotImplemented());
1128  }
1129  else
1130  {
1136  }
1137 }
1138 
1139 
1140 
1142 
1143 #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:200
const bool edge_matrices
Flag for storing matrices_in and matrices_out.
Definition: matrix_block.h:597
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:621
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:603
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
Contents is actually a matrix.
AnyData matrices_out
The matrix from the refinement edge to the interior of a level.
Definition: matrix_block.h:607
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:904
const value_type & block(size_type i) const
Definition: matrix_block.h:888
void Tvmult(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:822
void vmult(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:804
const value_type & block(size_type i) const
Definition: matrix_block.h:955
void unsubscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
Definition: subscriptor.cc:154
AnyData matrices_in
The matrix from the interior of a level to the refinement edge.
Definition: matrix_block.h:605
static ::ExceptionBase & ExcNotInitialized()
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
Definition: matrix_block.h:678
size_type column
Definition: matrix_block.h:312
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:600
void clear(bool really_clean=false)
Definition: matrix_block.h:871
const value_type & block_out(size_type i) const
Definition: matrix_block.h:987
void add(size_type row, size_type column, const std::string &name)
Definition: matrix_block.h:849
void reinit_edge(const MGLevelObject< BlockSparsityPattern > &sparsity)
MGMatrixBlockVector(const bool edge_matrices=false, const bool edge_flux_matrices=false)
Definition: matrix_block.h:914
size_type n() const
void subscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
Definition: subscriptor.cc:134
#define Assert(cond, exc)
Definition: exceptions.h:1419
AnyData flux_matrices_down
The DG flux from a level to the lower level.
Definition: matrix_block.h:609
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:630
const BlockIndices & get_row_indices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
unsigned int level
Definition: grid_out.cc:4355
std::size_t memory_consumption() const
Definition: matrix_block.h:839
BlockIndices column_indices
Definition: matrix_block.h:330
const value_type & block_in(size_type i) const
Definition: matrix_block.h:971
void vmult_add(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:813
std::shared_ptr< value_type > ptr_type
Definition: matrix_block.h:371
BlockIndices row_indices
Definition: matrix_block.h:324
typename FullMatrix< number > ::value_type value_type
Definition: matrix_block.h:123
MatrixType matrix
Definition: matrix_block.h:317
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:931
static ::ExceptionBase & ExcBlockIndexMismatch(size_type arg1, size_type arg2)
void clear(bool really_clean=false)
unsigned int size() const
Definition: matrix_block.h:923
void reinit_edge_flux(const MGLevelObject< BlockSparsityPattern > &sparsity)
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:860
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:358
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:656
AnyData flux_matrices_up
The DG flux from the lower level to a level.
Definition: matrix_block.h:611
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
static const bool value
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:831
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:307