deal.II version GIT relicensing-2877-g3bfda9281d 2025-03-20 21:30:00+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
91template <typename MatrixType>
93{
94public:
99
103 using value_type = typename MatrixType::value_type;
104
109
114
120
126
134 void
135 reinit(const BlockSparsityPattern &sparsity);
136
137 operator MatrixType &();
138 operator const MatrixType &() const;
139
144 void
145 add(const size_type i,
146 const size_type j,
147 const typename MatrixType::value_type value);
148
164 template <typename number>
165 void
166 add(const std::vector<size_type> &indices,
167 const FullMatrix<number> &full_matrix,
168 const bool elide_zero_values = true);
169
184 template <typename number>
185 void
186 add(const std::vector<size_type> &row_indices,
187 const std::vector<size_type> &col_indices,
188 const FullMatrix<number> &full_matrix,
189 const bool elide_zero_values = true);
190
207 template <typename number>
208 void
209 add(const size_type row_index,
210 const std::vector<size_type> &col_indices,
211 const std::vector<number> &values,
212 const bool elide_zero_values = true);
213
223 template <typename number>
224 void
226 const size_type n_cols,
227 const size_type *col_indices,
228 const number *values,
229 const bool elide_zero_values = true,
230 const bool col_indices_are_sorted = false);
231
237 template <typename VectorType>
238 void
239 vmult(VectorType &w, const VectorType &v) const;
240
246 template <typename VectorType>
247 void
248 vmult_add(VectorType &w, const VectorType &v) const;
249
255 template <typename VectorType>
256 void
257 Tvmult(VectorType &w, const VectorType &v) const;
258
264 template <typename VectorType>
265 void
266 Tvmult_add(VectorType &w, const VectorType &v) const;
267
271 std::size_t
273
279 size_type,
280 size_type,
281 << "Block index " << arg1 << " does not match " << arg2);
282
293
297 MatrixType matrix;
298
299private:
311
312 template <class OTHER_MatrixType>
313 friend void
314 ::internal::reinit(MatrixBlock<OTHER_MatrixType> &,
315 const BlockSparsityPattern &);
316
317 template <typename number>
318 friend void
320 const BlockSparsityPattern &p);
321};
322
323
332template <typename MatrixType>
334{
335public:
340
345
350 using ptr_type = std::shared_ptr<value_type>;
351
356 void
357 add(size_type row, size_type column, const std::string &name);
358
363 void
364 reinit(const BlockSparsityPattern &sparsity);
365
376 void
377 clear(bool really_clean = false);
378
382 std::size_t
384
388 const value_type &
389 block(size_type i) const;
390
394 value_type &
395 block(size_type i);
396
400 MatrixType &
401 matrix(size_type i);
402
406 using AnyData::name;
407 using AnyData::size;
408};
409
410
419template <typename MatrixType>
421{
422public:
427
441 MGMatrixBlockVector(const bool edge_matrices = false,
442 const bool edge_flux_matrices = false);
443
447 unsigned int
448 size() const;
449
455 void
456 add(size_type row, size_type column, const std::string &name);
457
464 void
473 void
481 void
483
494 void
495 clear(bool really_clean = false);
496
500 const value_type &
501 block(size_type i) const;
502
506 value_type &
507 block(size_type i);
508
513 const value_type &
514 block_in(size_type i) const;
515
519 value_type &
521
526 const value_type &
527 block_out(size_type i) const;
528
532 value_type &
534
539 const value_type &
540 block_up(size_type i) const;
541
545 value_type &
547
552 const value_type &
553 block_down(size_type i) const;
554
558 value_type &
560
564 std::size_t
566
567private:
569 void
571
573 const bool edge_matrices;
574
577
588};
589
590
591//----------------------------------------------------------------------//
592
593namespace internal
594{
595 template <typename MatrixType>
596 void
602
603
604 template <typename number>
605 void
607 const BlockSparsityPattern &p)
608 {
609 v.row_indices = p.get_row_indices();
610 v.column_indices = p.get_column_indices();
611 v.matrix.reinit(p.block(v.row, v.column));
612 }
613} // namespace internal
614
615
616template <typename MatrixType>
618 : row(numbers::invalid_size_type)
619 , column(numbers::invalid_size_type)
620{}
621
622
623template <typename MatrixType>
625 : row(i)
626 , column(j)
627{}
628
629
630template <typename MatrixType>
631inline void
633{
634 internal::reinit(*this, sparsity);
635}
636
637
638template <typename MatrixType>
640{
641 return matrix;
642}
643
644
645template <typename MatrixType>
646inline MatrixBlock<MatrixType>::operator const MatrixType &() const
647{
648 return matrix;
649}
650
651
652template <typename MatrixType>
653inline void
655 const size_type gj,
656 const typename MatrixType::value_type value)
657{
658 Assert(row_indices.size() != 0, ExcNotInitialized());
659 Assert(column_indices.size() != 0, ExcNotInitialized());
660
661 const std::pair<unsigned int, size_type> bi = row_indices.global_to_local(gi);
662 const std::pair<unsigned int, size_type> bj =
663 column_indices.global_to_local(gj);
664
665 Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
666 Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
667
668 matrix.add(bi.second, bj.second, value);
669}
670
671
672template <typename MatrixType>
673template <typename number>
674inline void
675MatrixBlock<MatrixType>::add(const std::vector<size_type> &r_indices,
676 const std::vector<size_type> &c_indices,
677 const FullMatrix<number> &values,
678 const bool elide_zero_values)
679{
680 Assert(row_indices.size() != 0, ExcNotInitialized());
681 Assert(column_indices.size() != 0, ExcNotInitialized());
682
683 AssertDimension(r_indices.size(), values.m());
684 AssertDimension(c_indices.size(), values.n());
685
686 for (size_type i = 0; i < row_indices.size(); ++i)
687 add(r_indices[i],
688 c_indices.size(),
689 c_indices.data(),
690 &values(i, 0),
691 elide_zero_values);
692}
693
694
695template <typename MatrixType>
696template <typename number>
697inline void
699 const size_type n_cols,
700 const size_type *col_indices,
701 const number *values,
702 const bool,
703 const bool)
704{
705 Assert(row_indices.size() != 0, ExcNotInitialized());
706 Assert(column_indices.size() != 0, ExcNotInitialized());
707
708 const std::pair<unsigned int, size_type> bi =
709 row_indices.global_to_local(b_row);
710
711 // In debug mode, we check whether
712 // all indices are in the correct
713 // block.
714
715 // Actually, for the time being, we
716 // leave it at this. While it may
717 // not be the most efficient way,
718 // it is at least thread safe.
719 // if constexpr (running_in_debug_mode()){
720 Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
721
722 for (size_type j = 0; j < n_cols; ++j)
723 {
724 const std::pair<unsigned int, size_type> bj =
725 column_indices.global_to_local(col_indices[j]);
726 Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
727
728 matrix.add(bi.second, bj.second, values[j]);
729 }
730 // }
731}
732
733
734template <typename MatrixType>
735template <typename number>
736inline void
737MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
738 const FullMatrix<number> &values,
739 const bool elide_zero_values)
740{
741 Assert(row_indices.size() != 0, ExcNotInitialized());
742 Assert(column_indices.size() != 0, ExcNotInitialized());
743
744 AssertDimension(indices.size(), values.m());
745 Assert(values.n() == values.m(), ExcNotQuadratic());
746
747 for (size_type i = 0; i < indices.size(); ++i)
748 add(indices[i],
749 indices.size(),
750 indices.data(),
751 &values(i, 0),
752 elide_zero_values);
753}
754
755
756
757template <typename MatrixType>
758template <typename number>
759inline void
761 const std::vector<size_type> &col_indices,
762 const std::vector<number> &values,
763 const bool elide_zero_values)
764{
765 Assert(row_indices.size() != 0, ExcNotInitialized());
766 Assert(column_indices.size() != 0, ExcNotInitialized());
767
768 AssertDimension(col_indices.size(), values.size());
769 add(row,
770 col_indices.size(),
771 col_indices.data(),
772 values.data(),
773 elide_zero_values);
774}
775
776
777template <typename MatrixType>
778template <typename VectorType>
779inline void
780MatrixBlock<MatrixType>::vmult(VectorType &w, const VectorType &v) const
781{
782 matrix.vmult(w, v);
783}
784
785
786template <typename MatrixType>
787template <typename VectorType>
788inline void
789MatrixBlock<MatrixType>::vmult_add(VectorType &w, const VectorType &v) const
790{
791 matrix.vmult_add(w, v);
792}
793
794
795template <typename MatrixType>
796template <typename VectorType>
797inline void
798MatrixBlock<MatrixType>::Tvmult(VectorType &w, const VectorType &v) const
799{
800 matrix.Tvmult(w, v);
801}
802
803
804template <typename MatrixType>
805template <typename VectorType>
806inline void
807MatrixBlock<MatrixType>::Tvmult_add(VectorType &w, const VectorType &v) const
808{
809 matrix.Tvmult_add(w, v);
810}
811
812
813template <typename MatrixType>
814inline std::size_t
816{
817 return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
818 sizeof(matrix));
819}
820
821//----------------------------------------------------------------------//
822
823template <typename MatrixType>
824inline void
826 size_type column,
827 const std::string &name)
828{
829 ptr_type p(new value_type(row, column));
830 AnyData::add(p, name);
831}
832
833
834template <typename MatrixType>
835inline void
837{
838 for (size_type i = 0; i < this->size(); ++i)
839 {
840 block(i).reinit(sparsity);
841 }
842}
843
844
845template <typename MatrixType>
846inline void
848{
849 if (really_clean)
850 {
852 }
853 else
854 {
855 for (size_type i = 0; i < this->size(); ++i)
856 matrix(i).clear();
857 }
858}
859
860
861
862template <typename MatrixType>
863inline const MatrixBlock<MatrixType> &
865{
866 return *this->read<ptr_type>(i);
867}
868
869
870template <typename MatrixType>
873{
874 return *this->entry<ptr_type>(i);
875}
876
877
878template <typename MatrixType>
879inline MatrixType &
881{
882 return this->entry<ptr_type>(i)->matrix;
883}
884
885
886
887//----------------------------------------------------------------------//
888
889template <typename MatrixType>
891 const bool f)
892 : edge_matrices(e)
893 , edge_flux_matrices(f)
894{}
895
896
897template <typename MatrixType>
898inline unsigned int
900{
901 return matrices.size();
902}
903
904
905template <typename MatrixType>
906inline void
908 size_type column,
909 const std::string &name)
910{
912 p[0].row = row;
913 p[0].column = column;
914
915 matrices.add(p, name);
916 if (edge_matrices)
917 {
918 matrices_in.add(p, name);
919 matrices_out.add(p, name);
920 }
921 if (edge_flux_matrices)
922 {
923 flux_matrices_up.add(p, name);
924 flux_matrices_down.add(p, name);
925 }
926}
927
928
929template <typename MatrixType>
932{
933 return *matrices.read<const MGLevelObject<MatrixType> *>(i);
934}
935
936
937template <typename MatrixType>
943
944
945template <typename MatrixType>
948{
949 return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
950}
951
952
953template <typename MatrixType>
956{
957 return *matrices_in.entry<MGLevelObject<MatrixType> *>(i);
958}
959
960
961template <typename MatrixType>
964{
965 return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
966}
967
968
969template <typename MatrixType>
972{
973 return *matrices_out.entry<MGLevelObject<MatrixType> *>(i);
974}
975
976
977template <typename MatrixType>
980{
981 return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
982}
983
984
985template <typename MatrixType>
988{
989 return *flux_matrices_up.entry<MGLevelObject<MatrixType> *>(i);
990}
991
992
993template <typename MatrixType>
996{
997 return *flux_matrices_down.read<const MGLevelObject<MatrixType> *>(i);
998}
999
1000
1001template <typename MatrixType>
1004{
1005 return *flux_matrices_down.entry<MGLevelObject<MatrixType> *>(i);
1006}
1007
1008
1009template <typename MatrixType>
1010inline void
1013{
1014 for (size_type i = 0; i < this->size(); ++i)
1015 {
1017 const size_type row = o[o.min_level()].row;
1018 const size_type col = o[o.min_level()].column;
1019
1020 o.resize(sparsity.min_level(), sparsity.max_level());
1021 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1022 {
1023 o[level].row = row;
1024 o[level].column = col;
1025 internal::reinit(o[level], sparsity[level]);
1026 }
1027 }
1028}
1029
1030
1031template <typename MatrixType>
1032inline void
1035{
1036 for (size_type i = 0; i < this->size(); ++i)
1037 {
1039 const size_type row = o[o.min_level()].row;
1040 const size_type col = o[o.min_level()].column;
1041
1042 block_in(i).resize(sparsity.min_level(), sparsity.max_level());
1043 block_out(i).resize(sparsity.min_level(), sparsity.max_level());
1044 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1045 {
1046 block_in(i)[level].row = row;
1047 block_in(i)[level].column = col;
1048 internal::reinit(block_in(i)[level], sparsity[level]);
1049 block_out(i)[level].row = row;
1050 block_out(i)[level].column = col;
1051 internal::reinit(block_out(i)[level], sparsity[level]);
1052 }
1053 }
1054}
1055
1056
1057template <typename MatrixType>
1058inline void
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_up(i).resize(sparsity.min_level(), sparsity.max_level());
1069 block_down(i).resize(sparsity.min_level(), sparsity.max_level());
1070 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1071 {
1072 block_up(i)[level].row = row;
1073 block_up(i)[level].column = col;
1074 internal::reinit(block_up(i)[level], sparsity[level]);
1075 block_down(i)[level].row = row;
1076 block_down(i)[level].column = col;
1077 internal::reinit(block_down(i)[level], sparsity[level]);
1078 }
1079 }
1080}
1081
1082
1083template <typename MatrixType>
1084inline void
1086{
1087 for (size_type i = 0; i < mo.size(); ++i)
1088 {
1091 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1092 o[level].matrix.clear();
1093 }
1094}
1095
1096
1097template <typename MatrixType>
1098inline void
1100{
1101 if (really_clean)
1102 {
1104 }
1105 else
1106 {
1107 clear_object(matrices);
1108 clear_object(matrices_in);
1109 clear_object(matrices_out);
1110 clear_object(flux_matrices_up);
1111 clear_object(flux_matrices_down);
1112 }
1113}
1114
1115
1116
1118
1119#endif
const std::string & name(const unsigned int i) const
Name of object at index.
Definition any_data.h:310
type entry(const std::string &name)
Access to stored data object by name.
Definition any_data.h:350
void add(type entry, const std::string &name)
Add a new data object.
Definition any_data.h:431
unsigned int size() const
Number of stored data objects.
Definition any_data.h:222
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:310
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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
#define DEAL_II_NOT_IMPLEMENTED()
unsigned int level
Definition grid_out.cc:4635
static ::ExceptionBase & ExcBlockIndexMismatch(size_type arg1, size_type arg2)
#define Assert(cond, exc)
#define DeclException2(Exception2, type1, type2, outsequence)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
std::size_t size
Definition mpi.cc:745
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:90