Reference documentation for deal.II version GIT relicensing-870-g40dcf97519 2024-06-20 23:40:02+00:00
\(\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// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2009 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_matrix_block_h
16#define dealii_matrix_block_h
17
18#include <deal.II/base/config.h>
19
21
25
30
31#include <memory>
32
33
35
36// Forward declarations
37#ifndef DOXYGEN
38template <typename MatrixType>
39class MatrixBlock;
40#endif
41
42namespace internal
43{
44 template <typename MatrixType>
45 void
47
48 template <typename number>
49 void
51 const BlockSparsityPattern &p);
52} // namespace internal
53
109template <typename MatrixType>
111{
112public:
117
121 using value_type = typename MatrixType::value_type;
122
127
132
138
144
152 void
153 reinit(const BlockSparsityPattern &sparsity);
154
155 operator MatrixType &();
156 operator const MatrixType &() const;
157
162 void
163 add(const size_type i,
164 const size_type j,
165 const typename MatrixType::value_type value);
166
182 template <typename number>
183 void
184 add(const std::vector<size_type> &indices,
185 const FullMatrix<number> &full_matrix,
186 const bool elide_zero_values = true);
187
202 template <typename number>
203 void
204 add(const std::vector<size_type> &row_indices,
205 const std::vector<size_type> &col_indices,
206 const FullMatrix<number> &full_matrix,
207 const bool elide_zero_values = true);
208
225 template <typename number>
226 void
227 add(const size_type row_index,
228 const std::vector<size_type> &col_indices,
229 const std::vector<number> &values,
230 const bool elide_zero_values = true);
231
241 template <typename number>
242 void
244 const size_type n_cols,
245 const size_type *col_indices,
246 const number *values,
247 const bool elide_zero_values = true,
248 const bool col_indices_are_sorted = false);
249
255 template <typename VectorType>
256 void
257 vmult(VectorType &w, const VectorType &v) const;
258
264 template <typename VectorType>
265 void
266 vmult_add(VectorType &w, const VectorType &v) const;
267
273 template <typename VectorType>
274 void
275 Tvmult(VectorType &w, const VectorType &v) const;
276
282 template <typename VectorType>
283 void
284 Tvmult_add(VectorType &w, const VectorType &v) const;
285
289 std::size_t
291
297 size_type,
298 size_type,
299 << "Block index " << arg1 << " does not match " << arg2);
300
311
315 MatrixType matrix;
316
317private:
329
330 template <class OTHER_MatrixType>
331 friend void
332 ::internal::reinit(MatrixBlock<OTHER_MatrixType> &,
333 const BlockSparsityPattern &);
334
335 template <typename number>
336 friend void
338 const BlockSparsityPattern &p);
339};
340
341
350template <typename MatrixType>
352{
353public:
358
363
368 using ptr_type = std::shared_ptr<value_type>;
369
374 void
375 add(size_type row, size_type column, const std::string &name);
376
381 void
382 reinit(const BlockSparsityPattern &sparsity);
383
394 void
395 clear(bool really_clean = false);
396
400 std::size_t
402
406 const value_type &
407 block(size_type i) const;
408
412 value_type &
413 block(size_type i);
414
418 MatrixType &
419 matrix(size_type i);
420
424 using AnyData::name;
425 using AnyData::size;
426 using AnyData::subscribe;
428};
429
430
439template <typename MatrixType>
441{
442public:
447
461 MGMatrixBlockVector(const bool edge_matrices = false,
462 const bool edge_flux_matrices = false);
463
467 unsigned int
468 size() const;
469
475 void
476 add(size_type row, size_type column, const std::string &name);
477
484 void
493 void
501 void
503
514 void
515 clear(bool really_clean = false);
516
520 const value_type &
521 block(size_type i) const;
522
526 value_type &
527 block(size_type i);
528
533 const value_type &
534 block_in(size_type i) const;
535
539 value_type &
541
546 const value_type &
547 block_out(size_type i) const;
548
552 value_type &
554
559 const value_type &
560 block_up(size_type i) const;
561
565 value_type &
567
572 const value_type &
573 block_down(size_type i) const;
574
578 value_type &
580
584 std::size_t
586
587private:
589 void
591
593 const bool edge_matrices;
594
597
608};
609
610
611//----------------------------------------------------------------------//
612
613namespace internal
614{
615 template <typename MatrixType>
616 void
622
623
624 template <typename number>
625 void
627 const BlockSparsityPattern &p)
628 {
629 v.row_indices = p.get_row_indices();
630 v.column_indices = p.get_column_indices();
631 v.matrix.reinit(p.block(v.row, v.column));
632 }
633} // namespace internal
634
635
636template <typename MatrixType>
638 : row(numbers::invalid_size_type)
639 , column(numbers::invalid_size_type)
640{}
641
642
643template <typename MatrixType>
645 : row(i)
646 , column(j)
647{}
648
649
650template <typename MatrixType>
651inline void
653{
654 internal::reinit(*this, sparsity);
655}
656
657
658template <typename MatrixType>
660{
661 return matrix;
662}
663
664
665template <typename MatrixType>
666inline MatrixBlock<MatrixType>::operator const MatrixType &() const
667{
668 return matrix;
669}
670
671
672template <typename MatrixType>
673inline void
675 const size_type gj,
676 const typename MatrixType::value_type value)
677{
678 Assert(row_indices.size() != 0, ExcNotInitialized());
679 Assert(column_indices.size() != 0, ExcNotInitialized());
680
681 const std::pair<unsigned int, size_type> bi = row_indices.global_to_local(gi);
682 const std::pair<unsigned int, size_type> bj =
683 column_indices.global_to_local(gj);
684
685 Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
686 Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
687
688 matrix.add(bi.second, bj.second, value);
689}
690
691
692template <typename MatrixType>
693template <typename number>
694inline void
695MatrixBlock<MatrixType>::add(const std::vector<size_type> &r_indices,
696 const std::vector<size_type> &c_indices,
697 const FullMatrix<number> &values,
698 const bool elide_zero_values)
699{
700 Assert(row_indices.size() != 0, ExcNotInitialized());
701 Assert(column_indices.size() != 0, ExcNotInitialized());
702
703 AssertDimension(r_indices.size(), values.m());
704 AssertDimension(c_indices.size(), values.n());
705
706 for (size_type i = 0; i < row_indices.size(); ++i)
707 add(r_indices[i],
708 c_indices.size(),
709 c_indices.data(),
710 &values(i, 0),
711 elide_zero_values);
712}
713
714
715template <typename MatrixType>
716template <typename number>
717inline void
719 const size_type n_cols,
720 const size_type *col_indices,
721 const number *values,
722 const bool,
723 const bool)
724{
725 Assert(row_indices.size() != 0, ExcNotInitialized());
726 Assert(column_indices.size() != 0, ExcNotInitialized());
727
728 const std::pair<unsigned int, size_type> bi =
729 row_indices.global_to_local(b_row);
730
731 // In debug mode, we check whether
732 // all indices are in the correct
733 // block.
734
735 // Actually, for the time being, we
736 // leave it at this. While it may
737 // not be the most efficient way,
738 // it is at least thread safe.
739 // #ifdef DEBUG
740 Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
741
742 for (size_type j = 0; j < n_cols; ++j)
743 {
744 const std::pair<unsigned int, size_type> bj =
745 column_indices.global_to_local(col_indices[j]);
746 Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
747
748 matrix.add(bi.second, bj.second, values[j]);
749 }
750 // #endif
751}
752
753
754template <typename MatrixType>
755template <typename number>
756inline void
757MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
758 const FullMatrix<number> &values,
759 const bool elide_zero_values)
760{
761 Assert(row_indices.size() != 0, ExcNotInitialized());
762 Assert(column_indices.size() != 0, ExcNotInitialized());
763
764 AssertDimension(indices.size(), values.m());
765 Assert(values.n() == values.m(), ExcNotQuadratic());
766
767 for (size_type i = 0; i < indices.size(); ++i)
768 add(indices[i],
769 indices.size(),
770 indices.data(),
771 &values(i, 0),
772 elide_zero_values);
773}
774
775
776
777template <typename MatrixType>
778template <typename number>
779inline void
781 const std::vector<size_type> &col_indices,
782 const std::vector<number> &values,
783 const bool elide_zero_values)
784{
785 Assert(row_indices.size() != 0, ExcNotInitialized());
786 Assert(column_indices.size() != 0, ExcNotInitialized());
787
788 AssertDimension(col_indices.size(), values.size());
789 add(row,
790 col_indices.size(),
791 col_indices.data(),
792 values.data(),
793 elide_zero_values);
794}
795
796
797template <typename MatrixType>
798template <typename VectorType>
799inline void
800MatrixBlock<MatrixType>::vmult(VectorType &w, const VectorType &v) const
801{
802 matrix.vmult(w, v);
803}
804
805
806template <typename MatrixType>
807template <typename VectorType>
808inline void
809MatrixBlock<MatrixType>::vmult_add(VectorType &w, const VectorType &v) const
810{
811 matrix.vmult_add(w, v);
812}
813
814
815template <typename MatrixType>
816template <typename VectorType>
817inline void
818MatrixBlock<MatrixType>::Tvmult(VectorType &w, const VectorType &v) const
819{
820 matrix.Tvmult(w, v);
821}
822
823
824template <typename MatrixType>
825template <typename VectorType>
826inline void
827MatrixBlock<MatrixType>::Tvmult_add(VectorType &w, const VectorType &v) const
828{
829 matrix.Tvmult_add(w, v);
830}
831
832
833template <typename MatrixType>
834inline std::size_t
836{
837 return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
838 sizeof(matrix));
839}
840
841//----------------------------------------------------------------------//
842
843template <typename MatrixType>
844inline void
846 size_type column,
847 const std::string &name)
848{
849 ptr_type p(new value_type(row, column));
850 AnyData::add(p, name);
851}
852
853
854template <typename MatrixType>
855inline void
857{
858 for (size_type i = 0; i < this->size(); ++i)
859 {
860 block(i).reinit(sparsity);
861 }
862}
863
864
865template <typename MatrixType>
866inline void
868{
869 if (really_clean)
870 {
872 }
873 else
874 {
875 for (size_type i = 0; i < this->size(); ++i)
876 matrix(i).clear();
877 }
878}
879
880
881
882template <typename MatrixType>
883inline const MatrixBlock<MatrixType> &
885{
886 return *this->read<ptr_type>(i);
887}
888
889
890template <typename MatrixType>
893{
894 return *this->entry<ptr_type>(i);
895}
896
897
898template <typename MatrixType>
899inline MatrixType &
901{
902 return this->entry<ptr_type>(i)->matrix;
903}
904
905
906
907//----------------------------------------------------------------------//
908
909template <typename MatrixType>
911 const bool f)
912 : edge_matrices(e)
913 , edge_flux_matrices(f)
914{}
915
916
917template <typename MatrixType>
918inline unsigned int
920{
921 return matrices.size();
922}
923
924
925template <typename MatrixType>
926inline void
928 size_type column,
929 const std::string &name)
930{
932 p[0].row = row;
933 p[0].column = column;
934
935 matrices.add(p, name);
936 if (edge_matrices)
937 {
938 matrices_in.add(p, name);
939 matrices_out.add(p, name);
940 }
941 if (edge_flux_matrices)
942 {
943 flux_matrices_up.add(p, name);
944 flux_matrices_down.add(p, name);
945 }
946}
947
948
949template <typename MatrixType>
952{
953 return *matrices.read<const MGLevelObject<MatrixType> *>(i);
954}
955
956
957template <typename MatrixType>
963
964
965template <typename MatrixType>
968{
969 return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
970}
971
972
973template <typename MatrixType>
976{
977 return *matrices_in.entry<MGLevelObject<MatrixType> *>(i);
978}
979
980
981template <typename MatrixType>
984{
985 return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
986}
987
988
989template <typename MatrixType>
992{
993 return *matrices_out.entry<MGLevelObject<MatrixType> *>(i);
994}
995
996
997template <typename MatrixType>
1000{
1001 return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
1002}
1003
1004
1005template <typename MatrixType>
1008{
1009 return *flux_matrices_up.entry<MGLevelObject<MatrixType> *>(i);
1010}
1011
1012
1013template <typename MatrixType>
1016{
1017 return *flux_matrices_down.read<const MGLevelObject<MatrixType> *>(i);
1018}
1019
1020
1021template <typename MatrixType>
1024{
1025 return *flux_matrices_down.entry<MGLevelObject<MatrixType> *>(i);
1026}
1027
1028
1029template <typename MatrixType>
1030inline void
1033{
1034 for (size_type i = 0; i < this->size(); ++i)
1035 {
1037 const size_type row = o[o.min_level()].row;
1038 const size_type col = o[o.min_level()].column;
1039
1040 o.resize(sparsity.min_level(), sparsity.max_level());
1041 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1042 {
1043 o[level].row = row;
1044 o[level].column = col;
1045 internal::reinit(o[level], sparsity[level]);
1046 }
1047 }
1048}
1049
1050
1051template <typename MatrixType>
1052inline void
1055{
1056 for (size_type i = 0; i < this->size(); ++i)
1057 {
1059 const size_type row = o[o.min_level()].row;
1060 const size_type col = o[o.min_level()].column;
1061
1062 block_in(i).resize(sparsity.min_level(), sparsity.max_level());
1063 block_out(i).resize(sparsity.min_level(), sparsity.max_level());
1064 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1065 {
1066 block_in(i)[level].row = row;
1067 block_in(i)[level].column = col;
1068 internal::reinit(block_in(i)[level], sparsity[level]);
1069 block_out(i)[level].row = row;
1070 block_out(i)[level].column = col;
1071 internal::reinit(block_out(i)[level], sparsity[level]);
1072 }
1073 }
1074}
1075
1076
1077template <typename MatrixType>
1078inline void
1081{
1082 for (size_type i = 0; i < this->size(); ++i)
1083 {
1085 const size_type row = o[o.min_level()].row;
1086 const size_type col = o[o.min_level()].column;
1087
1088 block_up(i).resize(sparsity.min_level(), sparsity.max_level());
1089 block_down(i).resize(sparsity.min_level(), sparsity.max_level());
1090 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1091 {
1092 block_up(i)[level].row = row;
1093 block_up(i)[level].column = col;
1094 internal::reinit(block_up(i)[level], sparsity[level]);
1095 block_down(i)[level].row = row;
1096 block_down(i)[level].column = col;
1097 internal::reinit(block_down(i)[level], sparsity[level]);
1098 }
1099 }
1100}
1101
1102
1103template <typename MatrixType>
1104inline void
1106{
1107 for (size_type i = 0; i < mo.size(); ++i)
1108 {
1111 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1112 o[level].matrix.clear();
1113 }
1114}
1115
1116
1117template <typename MatrixType>
1118inline void
1120{
1121 if (really_clean)
1122 {
1124 }
1125 else
1126 {
1127 clear_object(matrices);
1128 clear_object(matrices_in);
1129 clear_object(matrices_out);
1130 clear_object(flux_matrices_up);
1131 clear_object(flux_matrices_down);
1132 }
1133}
1134
1135
1136
1138
1139#endif
const std::string & name(const unsigned int i) const
Name of object at index.
Definition any_data.h:309
type entry(const std::string &name)
Access to stored data object by name.
Definition any_data.h:349
void add(type entry, const std::string &name)
Add a new data object.
Definition any_data.h:430
unsigned int size() const
Number of stored data objects.
Definition any_data.h:221
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:309
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
unsigned int level
Definition grid_out.cc:4626
static ::ExceptionBase & ExcBlockIndexMismatch(size_type arg1, size_type arg2)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
#define DEAL_II_NOT_IMPLEMENTED()
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
unsigned int global_dof_index
Definition types.h:81