Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
matrix_block.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2007 - 2018 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 
21 #include <deal.II/algorithms/any_data.h>
22 
23 #include <deal.II/base/memory_consumption.h>
24 #include <deal.II/base/mg_level_object.h>
25 #include <deal.II/base/smartpointer.h>
26 
27 #include <deal.II/lac/block_indices.h>
28 #include <deal.II/lac/block_sparsity_pattern.h>
29 #include <deal.II/lac/full_matrix.h>
30 #include <deal.II/lac/sparse_matrix.h>
31 
32 #include <memory>
33 
34 
35 DEAL_II_NAMESPACE_OPEN
36 
37 template <typename MatrixType>
39 
40 namespace internal
41 {
42  template <typename MatrixType>
43  void
44  reinit(MatrixBlock<MatrixType> &v, const BlockSparsityPattern &p);
45 
46  template <typename number>
47  void
48  reinit(MatrixBlock<::SparseMatrix<number>> &v,
49  const BlockSparsityPattern & p);
50 } // namespace internal
51 
108 template <typename MatrixType>
109 class MatrixBlock : public Subscriptor
110 {
111 public:
116 
120  using value_type = typename MatrixType::value_type;
121 
125  MatrixBlock();
126 
131 
137 
145  void
146  reinit(const BlockSparsityPattern &sparsity);
147 
148  operator MatrixType &();
149  operator const MatrixType &() const;
150 
155  void
156  add(const size_type i,
157  const size_type j,
158  const typename MatrixType::value_type value);
159 
175  template <typename number>
176  void
177  add(const std::vector<size_type> &indices,
178  const FullMatrix<number> & full_matrix,
179  const bool elide_zero_values = true);
180 
195  template <typename number>
196  void
197  add(const std::vector<size_type> &row_indices,
198  const std::vector<size_type> &col_indices,
199  const FullMatrix<number> & full_matrix,
200  const bool elide_zero_values = true);
201 
218  template <typename number>
219  void
220  add(const size_type row_index,
221  const std::vector<size_type> &col_indices,
222  const std::vector<number> & values,
223  const bool elide_zero_values = true);
224 
234  template <typename number>
235  void
236  add(const size_type row,
237  const size_type n_cols,
238  const size_type *col_indices,
239  const number * values,
240  const bool elide_zero_values = true,
241  const bool col_indices_are_sorted = false);
242 
248  template <class VectorType>
249  void
250  vmult(VectorType &w, const VectorType &v) const;
251 
257  template <class VectorType>
258  void
259  vmult_add(VectorType &w, const VectorType &v) const;
260 
266  template <class VectorType>
267  void
268  Tvmult(VectorType &w, const VectorType &v) const;
269 
275  template <class VectorType>
276  void
277  Tvmult_add(VectorType &w, const VectorType &v) const;
278 
282  std::size_t
283  memory_consumption() const;
284 
290  size_type,
291  size_type,
292  << "Block index " << arg1 << " does not match " << arg2);
293 
304 
308  MatrixType matrix;
309 
310 private:
322 
323  template <class OTHER_MatrixType>
324  friend void
325  ::internal::reinit(MatrixBlock<OTHER_MatrixType> &,
326  const BlockSparsityPattern &);
327 
328  template <typename number>
329  friend void
330  internal::reinit(MatrixBlock<::SparseMatrix<number>> &v,
331  const BlockSparsityPattern & p);
332 };
333 
334 
344 template <typename MatrixType>
345 class MatrixBlockVector : private AnyData
346 {
347 public:
352 
357 
362  using ptr_type = std::shared_ptr<value_type>;
363 
368  void
369  add(size_type row, size_type column, const std::string &name);
370 
375  void
376  reinit(const BlockSparsityPattern &sparsity);
377 
388  void
389  clear(bool really_clean = false);
390 
394  std::size_t
395  memory_consumption() const;
396 
400  const value_type &
401  block(size_type i) const;
402 
406  value_type &
407  block(size_type i);
408 
412  MatrixType &
413  matrix(size_type i);
414 
418  using AnyData::name;
419  using AnyData::size;
420  using AnyData::subscribe;
421  using AnyData::unsubscribe;
422 };
423 
424 
434 template <typename MatrixType>
436 {
437 public:
442 
456  MGMatrixBlockVector(const bool edge_matrices = false,
457  const bool edge_flux_matrices = false);
458 
462  unsigned int
463  size() const;
464 
470  void
471  add(size_type row, size_type column, const std::string &name);
472 
479  void
488  void
496  void
498 
509  void
510  clear(bool really_clean = false);
511 
515  const value_type &
516  block(size_type i) const;
517 
521  value_type &
522  block(size_type i);
523 
528  const value_type &
529  block_in(size_type i) const;
530 
534  value_type &
535  block_in(size_type i);
536 
541  const value_type &
542  block_out(size_type i) const;
543 
547  value_type &
548  block_out(size_type i);
549 
554  const value_type &
555  block_up(size_type i) const;
556 
560  value_type &
561  block_up(size_type i);
562 
567  const value_type &
568  block_down(size_type i) const;
569 
573  value_type &
575 
579  std::size_t
580  memory_consumption() const;
581 
582 private:
584  void
586 
588  const bool edge_matrices;
589 
591  const bool edge_flux_matrices;
592 
603 };
604 
605 
606 //----------------------------------------------------------------------//
607 
608 namespace internal
609 {
610  template <typename MatrixType>
611  void
612  reinit(MatrixBlock<MatrixType> &v, const BlockSparsityPattern &p)
613  {
616  }
617 
618 
619  template <typename number>
620  void
621  reinit(MatrixBlock<::SparseMatrix<number>> &v,
622  const BlockSparsityPattern & p)
623  {
626  v.matrix.reinit(p.block(v.row, v.column));
627  }
628 } // namespace internal
629 
630 
631 template <typename MatrixType>
633  : row(numbers::invalid_size_type)
634  , column(numbers::invalid_size_type)
635 {}
636 
637 
638 template <typename MatrixType>
640  : Subscriptor()
641  , row(M.row)
642  , column(M.column)
643  , matrix(M.matrix)
644  , row_indices(M.row_indices)
645  , column_indices(M.column_indices)
646 {}
647 
648 
649 template <typename MatrixType>
651  : row(i)
652  , column(j)
653 {}
654 
655 
656 template <typename MatrixType>
657 inline void
659 {
660  internal::reinit(*this, sparsity);
661 }
662 
663 
664 template <typename MatrixType>
665 inline MatrixBlock<MatrixType>::operator MatrixType &()
666 {
667  return matrix;
668 }
669 
670 
671 template <typename MatrixType>
672 inline MatrixBlock<MatrixType>::operator const MatrixType &() const
673 {
674  return matrix;
675 }
676 
677 
678 template <typename MatrixType>
679 inline void
681  const size_type gj,
682  const typename MatrixType::value_type value)
683 {
684  Assert(row_indices.size() != 0, ExcNotInitialized());
685  Assert(column_indices.size() != 0, ExcNotInitialized());
686 
687  const std::pair<unsigned int, size_type> bi = row_indices.global_to_local(gi);
688  const std::pair<unsigned int, size_type> bj =
689  column_indices.global_to_local(gj);
690 
691  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
692  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
693 
694  matrix.add(bi.second, bj.second, value);
695 }
696 
697 
698 template <typename MatrixType>
699 template <typename number>
700 inline void
701 MatrixBlock<MatrixType>::add(const std::vector<size_type> &r_indices,
702  const std::vector<size_type> &c_indices,
703  const FullMatrix<number> & values,
704  const bool elide_zero_values)
705 {
706  Assert(row_indices.size() != 0, ExcNotInitialized());
707  Assert(column_indices.size() != 0, ExcNotInitialized());
708 
709  AssertDimension(r_indices.size(), values.m());
710  AssertDimension(c_indices.size(), values.n());
711 
712  for (size_type i = 0; i < row_indices.size(); ++i)
713  add(r_indices[i],
714  c_indices.size(),
715  c_indices.data(),
716  &values(i, 0),
717  elide_zero_values);
718 }
719 
720 
721 template <typename MatrixType>
722 template <typename number>
723 inline void
725  const size_type n_cols,
726  const size_type *col_indices,
727  const number * values,
728  const bool,
729  const bool)
730 {
731  Assert(row_indices.size() != 0, ExcNotInitialized());
732  Assert(column_indices.size() != 0, ExcNotInitialized());
733 
734  const std::pair<unsigned int, size_type> bi =
735  row_indices.global_to_local(b_row);
736 
737  // In debug mode, we check whether
738  // all indices are in the correct
739  // block.
740 
741  // Actually, for the time being, we
742  // leave it at this. While it may
743  // not be the most efficient way,
744  // it is at least thread safe.
745  //#ifdef DEBUG
746  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
747 
748  for (size_type j = 0; j < n_cols; ++j)
749  {
750  const std::pair<unsigned int, size_type> bj =
751  column_indices.global_to_local(col_indices[j]);
752  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
753 
754  matrix.add(bi.second, bj.second, values[j]);
755  }
756  //#endif
757 }
758 
759 
760 template <typename MatrixType>
761 template <typename number>
762 inline void
763 MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
764  const FullMatrix<number> & values,
765  const bool elide_zero_values)
766 {
767  Assert(row_indices.size() != 0, ExcNotInitialized());
768  Assert(column_indices.size() != 0, ExcNotInitialized());
769 
770  AssertDimension(indices.size(), values.m());
771  Assert(values.n() == values.m(), ExcNotQuadratic());
772 
773  for (size_type i = 0; i < indices.size(); ++i)
774  add(indices[i],
775  indices.size(),
776  indices.data(),
777  &values(i, 0),
778  elide_zero_values);
779 }
780 
781 
782 
783 template <typename MatrixType>
784 template <typename number>
785 inline void
787  const std::vector<size_type> &col_indices,
788  const std::vector<number> & values,
789  const bool elide_zero_values)
790 {
791  Assert(row_indices.size() != 0, ExcNotInitialized());
792  Assert(column_indices.size() != 0, ExcNotInitialized());
793 
794  AssertDimension(col_indices.size(), values.size());
795  add(row,
796  col_indices.size(),
797  col_indices.data(),
798  values.data(),
799  elide_zero_values);
800 }
801 
802 
803 template <typename MatrixType>
804 template <class VectorType>
805 inline void
806 MatrixBlock<MatrixType>::vmult(VectorType &w, const VectorType &v) const
807 {
808  matrix.vmult(w, v);
809 }
810 
811 
812 template <typename MatrixType>
813 template <class VectorType>
814 inline void
815 MatrixBlock<MatrixType>::vmult_add(VectorType &w, const VectorType &v) const
816 {
817  matrix.vmult_add(w, v);
818 }
819 
820 
821 template <typename MatrixType>
822 template <class VectorType>
823 inline void
824 MatrixBlock<MatrixType>::Tvmult(VectorType &w, const VectorType &v) const
825 {
826  matrix.Tvmult(w, v);
827 }
828 
829 
830 template <typename MatrixType>
831 template <class VectorType>
832 inline void
833 MatrixBlock<MatrixType>::Tvmult_add(VectorType &w, const VectorType &v) const
834 {
835  matrix.Tvmult_add(w, v);
836 }
837 
838 
839 template <typename MatrixType>
840 inline std::size_t
842 {
843  return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
844  sizeof(matrix));
845 }
846 
847 //----------------------------------------------------------------------//
848 
849 template <typename MatrixType>
850 inline void
852  size_type column,
853  const std::string &name)
854 {
855  ptr_type p(new value_type(row, column));
856  AnyData::add(p, name);
857 }
858 
859 
860 template <typename MatrixType>
861 inline void
863 {
864  for (size_type i = 0; i < this->size(); ++i)
865  {
866  block(i).reinit(sparsity);
867  }
868 }
869 
870 
871 template <typename MatrixType>
872 inline void
874 {
875  if (really_clean)
876  {
877  Assert(false, ExcNotImplemented());
878  }
879  else
880  {
881  for (size_type i = 0; i < this->size(); ++i)
882  matrix(i).clear();
883  }
884 }
885 
886 
887 
888 template <typename MatrixType>
889 inline const MatrixBlock<MatrixType> &
891 {
892  return *this->read<ptr_type>(i);
893 }
894 
895 
896 template <typename MatrixType>
899 {
900  return *this->entry<ptr_type>(i);
901 }
902 
903 
904 template <typename MatrixType>
905 inline MatrixType &
907 {
908  return this->entry<ptr_type>(i)->matrix;
909 }
910 
911 
912 
913 //----------------------------------------------------------------------//
914 
915 template <typename MatrixType>
917  const bool f)
918  : edge_matrices(e)
919  , edge_flux_matrices(f)
920 {}
921 
922 
923 template <typename MatrixType>
924 inline unsigned int
926 {
927  return matrices.size();
928 }
929 
930 
931 template <typename MatrixType>
932 inline void
934  size_type column,
935  const std::string &name)
936 {
938  p[0].row = row;
939  p[0].column = column;
940 
941  matrices.add(p, name);
942  if (edge_matrices)
943  {
944  matrices_in.add(p, name);
945  matrices_out.add(p, name);
946  }
947  if (edge_flux_matrices)
948  {
949  flux_matrices_up.add(p, name);
950  flux_matrices_down.add(p, name);
951  }
952 }
953 
954 
955 template <typename MatrixType>
958 {
959  return *matrices.read<const MGLevelObject<MatrixType> *>(i);
960 }
961 
962 
963 template <typename MatrixType>
966 {
967  return *matrices.entry<MGLevelObject<MatrixType> *>(i);
968 }
969 
970 
971 template <typename MatrixType>
974 {
975  return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
976 }
977 
978 
979 template <typename MatrixType>
982 {
983  return *matrices_in.entry<MGLevelObject<MatrixType> *>(i);
984 }
985 
986 
987 template <typename MatrixType>
990 {
991  return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
992 }
993 
994 
995 template <typename MatrixType>
998 {
999  return *matrices_out.entry<MGLevelObject<MatrixType> *>(i);
1000 }
1001 
1002 
1003 template <typename MatrixType>
1006 {
1007  return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
1008 }
1009 
1010 
1011 template <typename MatrixType>
1014 {
1015  return *flux_matrices_up.entry<MGLevelObject<MatrixType> *>(i);
1016 }
1017 
1018 
1019 template <typename MatrixType>
1022 {
1023  return *flux_matrices_down.read<const MGLevelObject<MatrixType> *>(i);
1024 }
1025 
1026 
1027 template <typename MatrixType>
1030 {
1031  return *flux_matrices_down.entry<MGLevelObject<MatrixType> *>(i);
1032 }
1033 
1034 
1035 template <typename MatrixType>
1036 inline void
1038  const MGLevelObject<BlockSparsityPattern> &sparsity)
1039 {
1040  for (size_type i = 0; i < this->size(); ++i)
1041  {
1043  const size_type row = o[o.min_level()].row;
1044  const size_type col = o[o.min_level()].column;
1045 
1046  o.resize(sparsity.min_level(), sparsity.max_level());
1047  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1048  {
1049  o[level].row = row;
1050  o[level].column = col;
1051  internal::reinit(o[level], sparsity[level]);
1052  }
1053  }
1054 }
1055 
1056 
1057 template <typename MatrixType>
1058 inline void
1060  const MGLevelObject<BlockSparsityPattern> &sparsity)
1061 {
1062  for (size_type i = 0; i < this->size(); ++i)
1063  {
1065  const size_type row = o[o.min_level()].row;
1066  const size_type col = o[o.min_level()].column;
1067 
1068  block_in(i).resize(sparsity.min_level(), sparsity.max_level());
1069  block_out(i).resize(sparsity.min_level(), sparsity.max_level());
1070  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1071  {
1072  block_in(i)[level].row = row;
1073  block_in(i)[level].column = col;
1074  internal::reinit(block_in(i)[level], sparsity[level]);
1075  block_out(i)[level].row = row;
1076  block_out(i)[level].column = col;
1077  internal::reinit(block_out(i)[level], sparsity[level]);
1078  }
1079  }
1080 }
1081 
1082 
1083 template <typename MatrixType>
1084 inline void
1086  const MGLevelObject<BlockSparsityPattern> &sparsity)
1087 {
1088  for (size_type i = 0; i < this->size(); ++i)
1089  {
1091  const size_type row = o[o.min_level()].row;
1092  const size_type col = o[o.min_level()].column;
1093 
1094  block_up(i).resize(sparsity.min_level(), sparsity.max_level());
1095  block_down(i).resize(sparsity.min_level(), sparsity.max_level());
1096  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1097  {
1098  block_up(i)[level].row = row;
1099  block_up(i)[level].column = col;
1100  internal::reinit(block_up(i)[level], sparsity[level]);
1101  block_down(i)[level].row = row;
1102  block_down(i)[level].column = col;
1103  internal::reinit(block_down(i)[level], sparsity[level]);
1104  }
1105  }
1106 }
1107 
1108 
1109 template <typename MatrixType>
1110 inline void
1112 {
1113  for (size_type i = 0; i < mo.size(); ++i)
1114  {
1117  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1118  o[level].matrix.clear();
1119  }
1120 }
1121 
1122 
1123 template <typename MatrixType>
1124 inline void
1126 {
1127  if (really_clean)
1128  {
1129  Assert(false, ExcNotImplemented());
1130  }
1131  else
1132  {
1133  clear_object(matrices);
1134  clear_object(matrices_in);
1135  clear_object(matrices_out);
1136  clear_object(flux_matrices_up);
1137  clear_object(flux_matrices_down);
1138  }
1139 }
1140 
1141 
1142 
1143 DEAL_II_NAMESPACE_CLOSE
1144 
1145 #endif
void reinit_matrix(const MGLevelObject< BlockSparsityPattern > &sparsity)
unsigned int max_level() const
size_type m() const
const bool edge_matrices
Flag for storing matrices_in and matrices_out.
Definition: matrix_block.h:588
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:594
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
Contents is actually a matrix.
AnyData matrices_out
The matrix from the refinement edge to the interior of a level.
Definition: matrix_block.h:598
MatrixType & matrix(size_type i)
Definition: matrix_block.h:906
const value_type & block(size_type i) const
Definition: matrix_block.h:890
void Tvmult(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:824
void vmult(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:806
const value_type & block(size_type i) const
Definition: matrix_block.h:957
void unsubscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
Definition: subscriptor.cc:150
AnyData matrices_in
The matrix from the interior of a level to the refinement edge.
Definition: matrix_block.h:596
static ::ExceptionBase & ExcNotInitialized()
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
Definition: matrix_block.h:680
size_type column
Definition: matrix_block.h:303
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:591
void clear(bool really_clean=false)
Definition: matrix_block.h:873
const value_type & block_out(size_type i) const
Definition: matrix_block.h:989
void add(size_type row, size_type column, const std::string &name)
Definition: matrix_block.h:851
types::global_dof_index size_type
Definition: matrix_block.h:115
void reinit_edge(const MGLevelObject< BlockSparsityPattern > &sparsity)
MGMatrixBlockVector(const bool edge_matrices=false, const bool edge_flux_matrices=false)
Definition: matrix_block.h:916
size_type n() const
void subscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
Definition: subscriptor.cc:130
#define Assert(cond, exc)
Definition: exceptions.h:1407
AnyData flux_matrices_down
The DG flux from a level to the lower level.
Definition: matrix_block.h:600
const BlockIndices & get_row_indices() const
std::size_t memory_consumption() const
Definition: matrix_block.h:841
BlockIndices column_indices
Definition: matrix_block.h:321
const value_type & block_in(size_type i) const
Definition: matrix_block.h:973
void vmult_add(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:815
std::shared_ptr< value_type > ptr_type
Definition: matrix_block.h:362
BlockIndices row_indices
Definition: matrix_block.h:315
typename FullMatrix< number > ::value_type value_type
Definition: matrix_block.h:120
types::global_dof_index size_type
Definition: matrix_block.h:441
MatrixType matrix
Definition: matrix_block.h:308
static ::ExceptionBase & ExcNotQuadratic()
unsigned int global_dof_index
Definition: types.h:89
void add(size_type row, size_type column, const std::string &name)
Definition: matrix_block.h:933
static ::ExceptionBase & ExcBlockIndexMismatch(size_type arg1, size_type arg2)
void clear(bool really_clean=false)
unsigned int size() const
Definition: matrix_block.h:925
void reinit_edge_flux(const MGLevelObject< BlockSparsityPattern > &sparsity)
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:862
const value_type & block_up(size_type i) const
std::size_t memory_consumption() 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
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:658
types::global_dof_index size_type
Definition: matrix_block.h:351
AnyData flux_matrices_up
The DG flux from the lower level to a level.
Definition: matrix_block.h:602
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.
const BlockIndices & get_column_indices() const
void Tvmult_add(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:833
const std::string & name(const unsigned int i) const
Name of object at index.
Definition: any_data.h:311
std::size_t memory_consumption() const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
size_type row
Definition: matrix_block.h:298