Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
matrix_block.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2007 - 2023 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
39template <typename MatrixType>
40class MatrixBlock;
41#endif
42
43namespace internal
44{
45 template <typename MatrixType>
46 void
48
49 template <typename number>
50 void
52 const BlockSparsityPattern & p);
53} // namespace internal
54
110template <typename MatrixType>
112{
113public:
118
122 using value_type = typename MatrixType::value_type;
123
128
133
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
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
292
298 size_type,
299 size_type,
300 << "Block index " << arg1 << " does not match " << arg2);
301
312
316 MatrixType matrix;
317
318private:
330
331 template <class OTHER_MatrixType>
332 friend void
333 ::internal::reinit(MatrixBlock<OTHER_MatrixType> &,
334 const BlockSparsityPattern &);
335
336 template <typename number>
337 friend void
339 const BlockSparsityPattern & p);
340};
341
342
351template <typename MatrixType>
353{
354public:
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
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;
429};
430
431
440template <typename MatrixType>
442{
443public:
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
494 void
502 void
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 &
542
547 const value_type &
548 block_out(size_type i) const;
549
553 value_type &
555
560 const value_type &
561 block_up(size_type i) const;
562
566 value_type &
568
573 const value_type &
574 block_down(size_type i) const;
575
579 value_type &
581
585 std::size_t
587
588private:
590 void
592
594 const bool edge_matrices;
595
598
609};
610
611
612//----------------------------------------------------------------------//
613
614namespace 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
637template <typename MatrixType>
639 : row(numbers::invalid_size_type)
640 , column(numbers::invalid_size_type)
641{}
642
643
644template <typename MatrixType>
646 : row(i)
647 , column(j)
648{}
649
650
651template <typename MatrixType>
652inline void
654{
655 internal::reinit(*this, sparsity);
656}
657
658
659template <typename MatrixType>
661{
662 return matrix;
663}
664
665
666template <typename MatrixType>
667inline MatrixBlock<MatrixType>::operator const MatrixType &() const
668{
669 return matrix;
670}
671
672
673template <typename MatrixType>
674inline void
676 const size_type gj,
677 const typename MatrixType::value_type value)
678{
679 Assert(row_indices.size() != 0, ExcNotInitialized());
680 Assert(column_indices.size() != 0, ExcNotInitialized());
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 =
684 column_indices.global_to_local(gj);
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
693template <typename MatrixType>
694template <typename number>
695inline void
696MatrixBlock<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{
701 Assert(row_indices.size() != 0, ExcNotInitialized());
702 Assert(column_indices.size() != 0, ExcNotInitialized());
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
716template <typename MatrixType>
717template <typename number>
718inline void
720 const size_type n_cols,
721 const size_type *col_indices,
722 const number * values,
723 const bool,
724 const bool)
725{
726 Assert(row_indices.size() != 0, ExcNotInitialized());
727 Assert(column_indices.size() != 0, ExcNotInitialized());
728
729 const std::pair<unsigned int, size_type> bi =
730 row_indices.global_to_local(b_row);
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
755template <typename MatrixType>
756template <typename number>
757inline void
758MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
759 const FullMatrix<number> & values,
760 const bool elide_zero_values)
761{
762 Assert(row_indices.size() != 0, ExcNotInitialized());
763 Assert(column_indices.size() != 0, ExcNotInitialized());
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
778template <typename MatrixType>
779template <typename number>
780inline void
782 const std::vector<size_type> &col_indices,
783 const std::vector<number> & values,
784 const bool elide_zero_values)
785{
786 Assert(row_indices.size() != 0, ExcNotInitialized());
787 Assert(column_indices.size() != 0, ExcNotInitialized());
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
798template <typename MatrixType>
799template <class VectorType>
800inline void
801MatrixBlock<MatrixType>::vmult(VectorType &w, const VectorType &v) const
802{
803 matrix.vmult(w, v);
804}
805
806
807template <typename MatrixType>
808template <class VectorType>
809inline void
810MatrixBlock<MatrixType>::vmult_add(VectorType &w, const VectorType &v) const
811{
812 matrix.vmult_add(w, v);
813}
814
815
816template <typename MatrixType>
817template <class VectorType>
818inline void
819MatrixBlock<MatrixType>::Tvmult(VectorType &w, const VectorType &v) const
820{
821 matrix.Tvmult(w, v);
822}
823
824
825template <typename MatrixType>
826template <class VectorType>
827inline void
828MatrixBlock<MatrixType>::Tvmult_add(VectorType &w, const VectorType &v) const
829{
830 matrix.Tvmult_add(w, v);
831}
832
833
834template <typename MatrixType>
835inline std::size_t
837{
838 return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
839 sizeof(matrix));
840}
841
842//----------------------------------------------------------------------//
843
844template <typename MatrixType>
845inline void
847 size_type column,
848 const std::string &name)
849{
850 ptr_type p(new value_type(row, column));
851 AnyData::add(p, name);
852}
853
854
855template <typename MatrixType>
856inline void
858{
859 for (size_type i = 0; i < this->size(); ++i)
860 {
861 block(i).reinit(sparsity);
862 }
863}
864
865
866template <typename MatrixType>
867inline 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
883template <typename MatrixType>
884inline const MatrixBlock<MatrixType> &
886{
887 return *this->read<ptr_type>(i);
888}
889
890
891template <typename MatrixType>
894{
895 return *this->entry<ptr_type>(i);
896}
897
898
899template <typename MatrixType>
900inline MatrixType &
902{
903 return this->entry<ptr_type>(i)->matrix;
904}
905
906
907
908//----------------------------------------------------------------------//
909
910template <typename MatrixType>
912 const bool f)
913 : edge_matrices(e)
914 , edge_flux_matrices(f)
915{}
916
917
918template <typename MatrixType>
919inline unsigned int
921{
922 return matrices.size();
923}
924
925
926template <typename MatrixType>
927inline 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
950template <typename MatrixType>
953{
954 return *matrices.read<const MGLevelObject<MatrixType> *>(i);
955}
956
957
958template <typename MatrixType>
961{
962 return *matrices.entry<MGLevelObject<MatrixType> *>(i);
963}
964
965
966template <typename MatrixType>
969{
970 return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
971}
972
973
974template <typename MatrixType>
977{
978 return *matrices_in.entry<MGLevelObject<MatrixType> *>(i);
979}
980
981
982template <typename MatrixType>
985{
986 return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
987}
988
989
990template <typename MatrixType>
993{
994 return *matrices_out.entry<MGLevelObject<MatrixType> *>(i);
995}
996
997
998template <typename MatrixType>
1001{
1002 return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
1003}
1004
1005
1006template <typename MatrixType>
1009{
1010 return *flux_matrices_up.entry<MGLevelObject<MatrixType> *>(i);
1011}
1012
1013
1014template <typename MatrixType>
1017{
1018 return *flux_matrices_down.read<const MGLevelObject<MatrixType> *>(i);
1019}
1020
1021
1022template <typename MatrixType>
1025{
1026 return *flux_matrices_down.entry<MGLevelObject<MatrixType> *>(i);
1027}
1028
1029
1030template <typename MatrixType>
1031inline void
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
1052template <typename MatrixType>
1053inline void
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
1078template <typename MatrixType>
1079inline void
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
1104template <typename MatrixType>
1105inline 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
1118template <typename MatrixType>
1119inline void
1121{
1122 if (really_clean)
1123 {
1124 Assert(false, ExcNotImplemented());
1125 }
1126 else
1127 {
1128 clear_object(matrices);
1129 clear_object(matrices_in);
1130 clear_object(matrices_out);
1131 clear_object(flux_matrices_up);
1132 clear_object(flux_matrices_down);
1133 }
1134}
1135
1136
1137
1139
1140#endif
const std::string & name(const unsigned int i) const
Name of object at index.
Definition any_data.h:311
type entry(const std::string &name)
Access to stored data object by name.
Definition any_data.h:351
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 unsigned int n_blocks, const size_type n_elements_per_block)
SparsityPatternType & block(const size_type row, const size_type column)
const BlockIndices & get_column_indices() const
const BlockIndices & get_row_indices() const
void resize(const unsigned int new_minlevel, const unsigned int new_maxlevel, Args &&...args)
unsigned int max_level() const
unsigned int min_level() const
AnyData flux_matrices_up
The DG flux from the lower level to a level.
const value_type & block_down(size_type i) const
void clear_object(AnyData &)
Clear one of the matrix objects.
const value_type & block(size_type i) const
void clear(bool really_clean=false)
void reinit_edge_flux(const MGLevelObject< BlockSparsityPattern > &sparsity)
const bool edge_flux_matrices
Flag for storing flux_matrices_up and flux_matrices_down.
AnyData matrices_out
The matrix from the refinement edge to the interior of a level.
void reinit_matrix(const MGLevelObject< BlockSparsityPattern > &sparsity)
const value_type & block_in(size_type i) const
MGMatrixBlockVector(const bool edge_matrices=false, const bool edge_flux_matrices=false)
const value_type & block_out(size_type i) const
std::size_t memory_consumption() const
AnyData flux_matrices_down
The DG flux from a level to the lower level.
AnyData matrices
The level matrices.
void reinit_edge(const MGLevelObject< BlockSparsityPattern > &sparsity)
AnyData matrices_in
The matrix from the interior of a level to the refinement edge.
const bool edge_matrices
Flag for storing matrices_in and matrices_out.
void add(size_type row, size_type column, const std::string &name)
unsigned int size() const
const value_type & block_up(size_type i) const
const std::string & name(const unsigned int i) const
Definition any_data.h:311
std::shared_ptr< value_type > ptr_type
void add(size_type row, size_type column, const std::string &name)
const value_type & block(size_type i) const
void clear(bool really_clean=false)
std::size_t memory_consumption() const
void reinit(const BlockSparsityPattern &sparsity)
MatrixType & matrix(size_type i)
MatrixBlock< MatrixType > & operator=(const MatrixBlock< MatrixType > &)=default
void vmult_add(VectorType &w, const VectorType &v) const
void reinit(const BlockSparsityPattern &sparsity)
BlockIndices column_indices
MatrixBlock(size_type i, size_type j)
void add(const size_type i, const size_type j, const typename MatrixType::value_type value)
std::size_t memory_consumption() const
MatrixType matrix
void add(const std::vector< size_type > &row_indices, const std::vector< size_type > &col_indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
void Tvmult_add(VectorType &w, const VectorType &v) const
void add(const std::vector< size_type > &indices, const FullMatrix< number > &full_matrix, const bool elide_zero_values=true)
void add(const size_type row_index, const std::vector< size_type > &col_indices, const std::vector< number > &values, const bool elide_zero_values=true)
void vmult(VectorType &w, const VectorType &v) const
void Tvmult(VectorType &w, const VectorType &v) const
BlockIndices row_indices
typename MatrixType::value_type value_type
size_type row
MatrixBlock(const MatrixBlock< MatrixType > &M)=default
size_type column
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
friend void internal::reinit(MatrixBlock<::SparseMatrix< number > > &v, const BlockSparsityPattern &p)
void unsubscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
void subscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
unsigned int level
Definition grid_out.cc:4618
static ::ExceptionBase & ExcBlockIndexMismatch(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
unsigned int global_dof_index
Definition types.h:82