Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\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 // 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
47  reinit(MatrixBlock<MatrixType> &v, const BlockSparsityPattern &p);
48 
49  template <typename number>
50  void
51  reinit(MatrixBlock<::SparseMatrix<number>> &v,
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 
134 
140 
148  void
149  reinit(const BlockSparsityPattern &sparsity);
150 
151  operator MatrixType &();
152  operator const MatrixType &() const;
153 
158  void
159  add(const size_type i,
160  const size_type j,
161  const typename MatrixType::value_type value);
162 
178  template <typename number>
179  void
180  add(const std::vector<size_type> &indices,
181  const FullMatrix<number> & full_matrix,
182  const bool elide_zero_values = true);
183 
198  template <typename number>
199  void
200  add(const std::vector<size_type> &row_indices,
201  const std::vector<size_type> &col_indices,
202  const FullMatrix<number> & full_matrix,
203  const bool elide_zero_values = true);
204 
221  template <typename number>
222  void
223  add(const size_type row_index,
224  const std::vector<size_type> &col_indices,
225  const std::vector<number> & values,
226  const bool elide_zero_values = true);
227 
237  template <typename number>
238  void
239  add(const size_type row,
240  const size_type n_cols,
241  const size_type *col_indices,
242  const number * values,
243  const bool elide_zero_values = true,
244  const bool col_indices_are_sorted = false);
245 
251  template <class VectorType>
252  void
253  vmult(VectorType &w, const VectorType &v) const;
254 
260  template <class VectorType>
261  void
262  vmult_add(VectorType &w, const VectorType &v) const;
263 
269  template <class VectorType>
270  void
271  Tvmult(VectorType &w, const VectorType &v) const;
272 
278  template <class VectorType>
279  void
280  Tvmult_add(VectorType &w, const VectorType &v) const;
281 
285  std::size_t
286  memory_consumption() const;
287 
292  DeclException2(ExcBlockIndexMismatch,
293  size_type,
294  size_type,
295  << "Block index " << arg1 << " does not match " << arg2);
296 
307 
311  MatrixType matrix;
312 
313 private:
325 
326  template <class OTHER_MatrixType>
327  friend void
328  ::internal::reinit(MatrixBlock<OTHER_MatrixType> &,
329  const BlockSparsityPattern &);
330 
331  template <typename number>
332  friend void
333  internal::reinit(MatrixBlock<::SparseMatrix<number>> &v,
334  const BlockSparsityPattern & p);
335 };
336 
337 
347 template <typename MatrixType>
348 class MatrixBlockVector : private AnyData
349 {
350 public:
355 
360 
365  using ptr_type = std::shared_ptr<value_type>;
366 
371  void
372  add(size_type row, size_type column, const std::string &name);
373 
378  void
379  reinit(const BlockSparsityPattern &sparsity);
380 
391  void
392  clear(bool really_clean = false);
393 
397  std::size_t
398  memory_consumption() const;
399 
403  const value_type &
404  block(size_type i) const;
405 
409  value_type &
410  block(size_type i);
411 
415  MatrixType &
416  matrix(size_type i);
417 
421  using AnyData::name;
422  using AnyData::size;
423  using AnyData::subscribe;
424  using AnyData::unsubscribe;
425 };
426 
427 
437 template <typename MatrixType>
439 {
440 public:
445 
459  MGMatrixBlockVector(const bool edge_matrices = false,
460  const bool edge_flux_matrices = false);
461 
465  unsigned int
466  size() const;
467 
473  void
474  add(size_type row, size_type column, const std::string &name);
475 
482  void
483  reinit_matrix(const MGLevelObject<BlockSparsityPattern> &sparsity);
491  void
492  reinit_edge(const MGLevelObject<BlockSparsityPattern> &sparsity);
499  void
500  reinit_edge_flux(const MGLevelObject<BlockSparsityPattern> &sparsity);
501 
512  void
513  clear(bool really_clean = false);
514 
518  const value_type &
519  block(size_type i) const;
520 
524  value_type &
525  block(size_type i);
526 
531  const value_type &
532  block_in(size_type i) const;
533 
537  value_type &
538  block_in(size_type i);
539 
544  const value_type &
545  block_out(size_type i) const;
546 
550  value_type &
551  block_out(size_type i);
552 
557  const value_type &
558  block_up(size_type i) const;
559 
563  value_type &
564  block_up(size_type i);
565 
570  const value_type &
571  block_down(size_type i) const;
572 
576  value_type &
577  block_down(size_type i);
578 
582  std::size_t
583  memory_consumption() const;
584 
585 private:
587  void
588  clear_object(AnyData &);
589 
591  const bool edge_matrices;
592 
594  const bool edge_flux_matrices;
595 
606 };
607 
608 
609 //----------------------------------------------------------------------//
610 
611 namespace internal
612 {
613  template <typename MatrixType>
614  void
615  reinit(MatrixBlock<MatrixType> &v, const BlockSparsityPattern &p)
616  {
619  }
620 
621 
622  template <typename number>
623  void
624  reinit(MatrixBlock<::SparseMatrix<number>> &v,
625  const BlockSparsityPattern & p)
626  {
629  v.matrix.reinit(p.block(v.row, v.column));
630  }
631 } // namespace internal
632 
633 
634 template <typename MatrixType>
636  : row(numbers::invalid_size_type)
637  , column(numbers::invalid_size_type)
638 {}
639 
640 
641 template <typename MatrixType>
643  : Subscriptor()
644  , row(M.row)
645  , column(M.column)
646  , matrix(M.matrix)
649 {}
650 
651 
652 template <typename MatrixType>
654  : row(i)
655  , column(j)
656 {}
657 
658 
659 template <typename MatrixType>
660 inline void
662 {
663  internal::reinit(*this, sparsity);
664 }
665 
666 
667 template <typename MatrixType>
668 inline MatrixBlock<MatrixType>::operator MatrixType &()
669 {
670  return matrix;
671 }
672 
673 
674 template <typename MatrixType>
675 inline MatrixBlock<MatrixType>::operator const MatrixType &() const
676 {
677  return matrix;
678 }
679 
680 
681 template <typename MatrixType>
682 inline void
684  const size_type gj,
685  const typename MatrixType::value_type value)
686 {
689 
690  const std::pair<unsigned int, size_type> bi = row_indices.global_to_local(gi);
691  const std::pair<unsigned int, size_type> bj =
693 
694  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
695  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
696 
697  matrix.add(bi.second, bj.second, value);
698 }
699 
700 
701 template <typename MatrixType>
702 template <typename number>
703 inline void
704 MatrixBlock<MatrixType>::add(const std::vector<size_type> &r_indices,
705  const std::vector<size_type> &c_indices,
706  const FullMatrix<number> & values,
707  const bool elide_zero_values)
708 {
711 
712  AssertDimension(r_indices.size(), values.m());
713  AssertDimension(c_indices.size(), values.n());
714 
715  for (size_type i = 0; i < row_indices.size(); ++i)
716  add(r_indices[i],
717  c_indices.size(),
718  c_indices.data(),
719  &values(i, 0),
720  elide_zero_values);
721 }
722 
723 
724 template <typename MatrixType>
725 template <typename number>
726 inline void
728  const size_type n_cols,
729  const size_type *col_indices,
730  const number * values,
731  const bool,
732  const bool)
733 {
736 
737  const std::pair<unsigned int, size_type> bi =
739 
740  // In debug mode, we check whether
741  // all indices are in the correct
742  // block.
743 
744  // Actually, for the time being, we
745  // leave it at this. While it may
746  // not be the most efficient way,
747  // it is at least thread safe.
748  //#ifdef DEBUG
749  Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
750 
751  for (size_type j = 0; j < n_cols; ++j)
752  {
753  const std::pair<unsigned int, size_type> bj =
754  column_indices.global_to_local(col_indices[j]);
755  Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
756 
757  matrix.add(bi.second, bj.second, values[j]);
758  }
759  //#endif
760 }
761 
762 
763 template <typename MatrixType>
764 template <typename number>
765 inline void
766 MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
767  const FullMatrix<number> & values,
768  const bool elide_zero_values)
769 {
772 
773  AssertDimension(indices.size(), values.m());
774  Assert(values.n() == values.m(), ExcNotQuadratic());
775 
776  for (size_type i = 0; i < indices.size(); ++i)
777  add(indices[i],
778  indices.size(),
779  indices.data(),
780  &values(i, 0),
781  elide_zero_values);
782 }
783 
784 
785 
786 template <typename MatrixType>
787 template <typename number>
788 inline void
790  const std::vector<size_type> &col_indices,
791  const std::vector<number> & values,
792  const bool elide_zero_values)
793 {
796 
797  AssertDimension(col_indices.size(), values.size());
798  add(row,
799  col_indices.size(),
800  col_indices.data(),
801  values.data(),
802  elide_zero_values);
803 }
804 
805 
806 template <typename MatrixType>
807 template <class VectorType>
808 inline void
809 MatrixBlock<MatrixType>::vmult(VectorType &w, const VectorType &v) const
810 {
811  matrix.vmult(w, v);
812 }
813 
814 
815 template <typename MatrixType>
816 template <class VectorType>
817 inline void
818 MatrixBlock<MatrixType>::vmult_add(VectorType &w, const VectorType &v) const
819 {
820  matrix.vmult_add(w, v);
821 }
822 
823 
824 template <typename MatrixType>
825 template <class VectorType>
826 inline void
827 MatrixBlock<MatrixType>::Tvmult(VectorType &w, const VectorType &v) const
828 {
829  matrix.Tvmult(w, v);
830 }
831 
832 
833 template <typename MatrixType>
834 template <class VectorType>
835 inline void
836 MatrixBlock<MatrixType>::Tvmult_add(VectorType &w, const VectorType &v) const
837 {
838  matrix.Tvmult_add(w, v);
839 }
840 
841 
842 template <typename MatrixType>
843 inline std::size_t
845 {
846  return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
847  sizeof(matrix));
848 }
849 
850 //----------------------------------------------------------------------//
851 
852 template <typename MatrixType>
853 inline void
856  const std::string &name)
857 {
858  ptr_type p(new value_type(row, column));
859  AnyData::add(p, name);
860 }
861 
862 
863 template <typename MatrixType>
864 inline void
866 {
867  for (size_type i = 0; i < this->size(); ++i)
868  {
869  block(i).reinit(sparsity);
870  }
871 }
872 
873 
874 template <typename MatrixType>
875 inline void
877 {
878  if (really_clean)
879  {
880  Assert(false, ExcNotImplemented());
881  }
882  else
883  {
884  for (size_type i = 0; i < this->size(); ++i)
885  matrix(i).clear();
886  }
887 }
888 
889 
890 
891 template <typename MatrixType>
892 inline const MatrixBlock<MatrixType> &
894 {
895  return *this->read<ptr_type>(i);
896 }
897 
898 
899 template <typename MatrixType>
902 {
903  return *this->entry<ptr_type>(i);
904 }
905 
906 
907 template <typename MatrixType>
908 inline MatrixType &
910 {
911  return this->entry<ptr_type>(i)->matrix;
912 }
913 
914 
915 
916 //----------------------------------------------------------------------//
917 
918 template <typename MatrixType>
920  const bool f)
921  : edge_matrices(e)
922  , edge_flux_matrices(f)
923 {}
924 
925 
926 template <typename MatrixType>
927 inline unsigned int
929 {
930  return matrices.size();
931 }
932 
933 
934 template <typename MatrixType>
935 inline void
937  size_type column,
938  const std::string &name)
939 {
941  p[0].row = row;
942  p[0].column = column;
943 
944  matrices.add(p, name);
945  if (edge_matrices)
946  {
947  matrices_in.add(p, name);
948  matrices_out.add(p, name);
949  }
950  if (edge_flux_matrices)
951  {
952  flux_matrices_up.add(p, name);
953  flux_matrices_down.add(p, name);
954  }
955 }
956 
957 
958 template <typename MatrixType>
961 {
962  return *matrices.read<const MGLevelObject<MatrixType> *>(i);
963 }
964 
965 
966 template <typename MatrixType>
969 {
971 }
972 
973 
974 template <typename MatrixType>
977 {
978  return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
979 }
980 
981 
982 template <typename MatrixType>
985 {
987 }
988 
989 
990 template <typename MatrixType>
993 {
994  return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
995 }
996 
997 
998 template <typename MatrixType>
1001 {
1003 }
1004 
1005 
1006 template <typename MatrixType>
1009 {
1010  return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
1011 }
1012 
1013 
1014 template <typename MatrixType>
1017 {
1019 }
1020 
1021 
1022 template <typename MatrixType>
1025 {
1027 }
1028 
1029 
1030 template <typename MatrixType>
1033 {
1035 }
1036 
1037 
1038 template <typename MatrixType>
1039 inline void
1041  const MGLevelObject<BlockSparsityPattern> &sparsity)
1042 {
1043  for (size_type i = 0; i < this->size(); ++i)
1044  {
1046  const size_type row = o[o.min_level()].row;
1047  const size_type col = o[o.min_level()].column;
1048 
1049  o.resize(sparsity.min_level(), sparsity.max_level());
1050  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1051  {
1052  o[level].row = row;
1053  o[level].column = col;
1054  internal::reinit(o[level], sparsity[level]);
1055  }
1056  }
1057 }
1058 
1059 
1060 template <typename MatrixType>
1061 inline void
1063  const MGLevelObject<BlockSparsityPattern> &sparsity)
1064 {
1065  for (size_type i = 0; i < this->size(); ++i)
1066  {
1068  const size_type row = o[o.min_level()].row;
1069  const size_type col = o[o.min_level()].column;
1070 
1071  block_in(i).resize(sparsity.min_level(), sparsity.max_level());
1072  block_out(i).resize(sparsity.min_level(), sparsity.max_level());
1073  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1074  {
1075  block_in(i)[level].row = row;
1076  block_in(i)[level].column = col;
1077  internal::reinit(block_in(i)[level], sparsity[level]);
1078  block_out(i)[level].row = row;
1079  block_out(i)[level].column = col;
1080  internal::reinit(block_out(i)[level], sparsity[level]);
1081  }
1082  }
1083 }
1084 
1085 
1086 template <typename MatrixType>
1087 inline void
1089  const MGLevelObject<BlockSparsityPattern> &sparsity)
1090 {
1091  for (size_type i = 0; i < this->size(); ++i)
1092  {
1094  const size_type row = o[o.min_level()].row;
1095  const size_type col = o[o.min_level()].column;
1096 
1097  block_up(i).resize(sparsity.min_level(), sparsity.max_level());
1098  block_down(i).resize(sparsity.min_level(), sparsity.max_level());
1099  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1100  {
1101  block_up(i)[level].row = row;
1102  block_up(i)[level].column = col;
1103  internal::reinit(block_up(i)[level], sparsity[level]);
1104  block_down(i)[level].row = row;
1105  block_down(i)[level].column = col;
1106  internal::reinit(block_down(i)[level], sparsity[level]);
1107  }
1108  }
1109 }
1110 
1111 
1112 template <typename MatrixType>
1113 inline void
1115 {
1116  for (size_type i = 0; i < mo.size(); ++i)
1117  {
1120  for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1121  o[level].matrix.clear();
1122  }
1123 }
1124 
1125 
1126 template <typename MatrixType>
1127 inline void
1129 {
1130  if (really_clean)
1131  {
1132  Assert(false, ExcNotImplemented());
1133  }
1134  else
1135  {
1141  }
1142 }
1143 
1144 
1145 
1146 DEAL_II_NAMESPACE_CLOSE
1147 
1148 #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:591
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:597
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
AnyData matrices_out
The matrix from the refinement edge to the interior of a level.
Definition: matrix_block.h:601
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:909
const value_type & block(size_type i) const
Definition: matrix_block.h:893
void Tvmult(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:827
void vmult(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:809
const value_type & block(size_type i) const
Definition: matrix_block.h:960
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:599
static ::ExceptionBase & ExcNotInitialized()
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
Definition: matrix_block.h:683
size_type column
Definition: matrix_block.h:306
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:594
void clear(bool really_clean=false)
Definition: matrix_block.h:876
const value_type & block_out(size_type i) const
Definition: matrix_block.h:992
void add(size_type row, size_type column, const std::string &name)
Definition: matrix_block.h:854
types::global_dof_index size_type
Definition: matrix_block.h:118
void reinit_edge(const MGLevelObject< BlockSparsityPattern > &sparsity)
MGMatrixBlockVector(const bool edge_matrices=false, const bool edge_flux_matrices=false)
Definition: matrix_block.h:919
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:1407
AnyData flux_matrices_down
The DG flux from a level to the lower level.
Definition: matrix_block.h:603
const BlockIndices & get_row_indices() const
std::size_t memory_consumption() const
Definition: matrix_block.h:844
BlockIndices column_indices
Definition: matrix_block.h:324
const value_type & block_in(size_type i) const
Definition: matrix_block.h:976
void vmult_add(VectorType &w, const VectorType &v) const
Definition: matrix_block.h:818
std::shared_ptr< value_type > ptr_type
Definition: matrix_block.h:365
BlockIndices row_indices
Definition: matrix_block.h:318
typename FullMatrix< number > ::value_type value_type
Definition: matrix_block.h:123
types::global_dof_index size_type
Definition: matrix_block.h:444
MatrixType matrix
Definition: matrix_block.h:311
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:936
static ::ExceptionBase & ExcBlockIndexMismatch(size_type arg1, size_type arg2)
void clear(bool really_clean=false)
unsigned int size() const
Definition: matrix_block.h:928
void reinit_edge_flux(const MGLevelObject< BlockSparsityPattern > &sparsity)
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:865
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
void reinit(const BlockSparsityPattern &sparsity)
Definition: matrix_block.h:661
types::global_dof_index size_type
Definition: matrix_block.h:354
AnyData flux_matrices_up
The DG flux from the lower level to a level.
Definition: matrix_block.h:605
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:836
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:301