deal.II version GIT relicensing-2206-gaa53ff9447 2024-12-02 09:10: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 using AnyData::subscribe;
410};
411
412
421template <typename MatrixType>
423{
424public:
429
443 MGMatrixBlockVector(const bool edge_matrices = false,
444 const bool edge_flux_matrices = false);
445
449 unsigned int
450 size() const;
451
457 void
458 add(size_type row, size_type column, const std::string &name);
459
466 void
475 void
483 void
485
496 void
497 clear(bool really_clean = false);
498
502 const value_type &
503 block(size_type i) const;
504
508 value_type &
509 block(size_type i);
510
515 const value_type &
516 block_in(size_type i) const;
517
521 value_type &
523
528 const value_type &
529 block_out(size_type i) const;
530
534 value_type &
536
541 const value_type &
542 block_up(size_type i) const;
543
547 value_type &
549
554 const value_type &
555 block_down(size_type i) const;
556
560 value_type &
562
566 std::size_t
568
569private:
571 void
573
575 const bool edge_matrices;
576
579
590};
591
592
593//----------------------------------------------------------------------//
594
595namespace internal
596{
597 template <typename MatrixType>
598 void
604
605
606 template <typename number>
607 void
609 const BlockSparsityPattern &p)
610 {
611 v.row_indices = p.get_row_indices();
612 v.column_indices = p.get_column_indices();
613 v.matrix.reinit(p.block(v.row, v.column));
614 }
615} // namespace internal
616
617
618template <typename MatrixType>
620 : row(numbers::invalid_size_type)
621 , column(numbers::invalid_size_type)
622{}
623
624
625template <typename MatrixType>
627 : row(i)
628 , column(j)
629{}
630
631
632template <typename MatrixType>
633inline void
635{
636 internal::reinit(*this, sparsity);
637}
638
639
640template <typename MatrixType>
642{
643 return matrix;
644}
645
646
647template <typename MatrixType>
648inline MatrixBlock<MatrixType>::operator const MatrixType &() const
649{
650 return matrix;
651}
652
653
654template <typename MatrixType>
655inline void
657 const size_type gj,
658 const typename MatrixType::value_type value)
659{
660 Assert(row_indices.size() != 0, ExcNotInitialized());
661 Assert(column_indices.size() != 0, ExcNotInitialized());
662
663 const std::pair<unsigned int, size_type> bi = row_indices.global_to_local(gi);
664 const std::pair<unsigned int, size_type> bj =
665 column_indices.global_to_local(gj);
666
667 Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
668 Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
669
670 matrix.add(bi.second, bj.second, value);
671}
672
673
674template <typename MatrixType>
675template <typename number>
676inline void
677MatrixBlock<MatrixType>::add(const std::vector<size_type> &r_indices,
678 const std::vector<size_type> &c_indices,
679 const FullMatrix<number> &values,
680 const bool elide_zero_values)
681{
682 Assert(row_indices.size() != 0, ExcNotInitialized());
683 Assert(column_indices.size() != 0, ExcNotInitialized());
684
685 AssertDimension(r_indices.size(), values.m());
686 AssertDimension(c_indices.size(), values.n());
687
688 for (size_type i = 0; i < row_indices.size(); ++i)
689 add(r_indices[i],
690 c_indices.size(),
691 c_indices.data(),
692 &values(i, 0),
693 elide_zero_values);
694}
695
696
697template <typename MatrixType>
698template <typename number>
699inline void
701 const size_type n_cols,
702 const size_type *col_indices,
703 const number *values,
704 const bool,
705 const bool)
706{
707 Assert(row_indices.size() != 0, ExcNotInitialized());
708 Assert(column_indices.size() != 0, ExcNotInitialized());
709
710 const std::pair<unsigned int, size_type> bi =
711 row_indices.global_to_local(b_row);
712
713 // In debug mode, we check whether
714 // all indices are in the correct
715 // block.
716
717 // Actually, for the time being, we
718 // leave it at this. While it may
719 // not be the most efficient way,
720 // it is at least thread safe.
721 // #ifdef DEBUG
722 Assert(bi.first == row, ExcBlockIndexMismatch(bi.first, row));
723
724 for (size_type j = 0; j < n_cols; ++j)
725 {
726 const std::pair<unsigned int, size_type> bj =
727 column_indices.global_to_local(col_indices[j]);
728 Assert(bj.first == column, ExcBlockIndexMismatch(bj.first, column));
729
730 matrix.add(bi.second, bj.second, values[j]);
731 }
732 // #endif
733}
734
735
736template <typename MatrixType>
737template <typename number>
738inline void
739MatrixBlock<MatrixType>::add(const std::vector<size_type> &indices,
740 const FullMatrix<number> &values,
741 const bool elide_zero_values)
742{
743 Assert(row_indices.size() != 0, ExcNotInitialized());
744 Assert(column_indices.size() != 0, ExcNotInitialized());
745
746 AssertDimension(indices.size(), values.m());
747 Assert(values.n() == values.m(), ExcNotQuadratic());
748
749 for (size_type i = 0; i < indices.size(); ++i)
750 add(indices[i],
751 indices.size(),
752 indices.data(),
753 &values(i, 0),
754 elide_zero_values);
755}
756
757
758
759template <typename MatrixType>
760template <typename number>
761inline void
763 const std::vector<size_type> &col_indices,
764 const std::vector<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(col_indices.size(), values.size());
771 add(row,
772 col_indices.size(),
773 col_indices.data(),
774 values.data(),
775 elide_zero_values);
776}
777
778
779template <typename MatrixType>
780template <typename VectorType>
781inline void
782MatrixBlock<MatrixType>::vmult(VectorType &w, const VectorType &v) const
783{
784 matrix.vmult(w, v);
785}
786
787
788template <typename MatrixType>
789template <typename VectorType>
790inline void
791MatrixBlock<MatrixType>::vmult_add(VectorType &w, const VectorType &v) const
792{
793 matrix.vmult_add(w, v);
794}
795
796
797template <typename MatrixType>
798template <typename VectorType>
799inline void
800MatrixBlock<MatrixType>::Tvmult(VectorType &w, const VectorType &v) const
801{
802 matrix.Tvmult(w, v);
803}
804
805
806template <typename MatrixType>
807template <typename VectorType>
808inline void
809MatrixBlock<MatrixType>::Tvmult_add(VectorType &w, const VectorType &v) const
810{
811 matrix.Tvmult_add(w, v);
812}
813
814
815template <typename MatrixType>
816inline std::size_t
818{
819 return (sizeof(*this) + MemoryConsumption::memory_consumption(matrix) -
820 sizeof(matrix));
821}
822
823//----------------------------------------------------------------------//
824
825template <typename MatrixType>
826inline void
828 size_type column,
829 const std::string &name)
830{
831 ptr_type p(new value_type(row, column));
832 AnyData::add(p, name);
833}
834
835
836template <typename MatrixType>
837inline void
839{
840 for (size_type i = 0; i < this->size(); ++i)
841 {
842 block(i).reinit(sparsity);
843 }
844}
845
846
847template <typename MatrixType>
848inline void
850{
851 if (really_clean)
852 {
854 }
855 else
856 {
857 for (size_type i = 0; i < this->size(); ++i)
858 matrix(i).clear();
859 }
860}
861
862
863
864template <typename MatrixType>
865inline const MatrixBlock<MatrixType> &
867{
868 return *this->read<ptr_type>(i);
869}
870
871
872template <typename MatrixType>
875{
876 return *this->entry<ptr_type>(i);
877}
878
879
880template <typename MatrixType>
881inline MatrixType &
883{
884 return this->entry<ptr_type>(i)->matrix;
885}
886
887
888
889//----------------------------------------------------------------------//
890
891template <typename MatrixType>
893 const bool f)
894 : edge_matrices(e)
895 , edge_flux_matrices(f)
896{}
897
898
899template <typename MatrixType>
900inline unsigned int
902{
903 return matrices.size();
904}
905
906
907template <typename MatrixType>
908inline void
910 size_type column,
911 const std::string &name)
912{
914 p[0].row = row;
915 p[0].column = column;
916
917 matrices.add(p, name);
918 if (edge_matrices)
919 {
920 matrices_in.add(p, name);
921 matrices_out.add(p, name);
922 }
923 if (edge_flux_matrices)
924 {
925 flux_matrices_up.add(p, name);
926 flux_matrices_down.add(p, name);
927 }
928}
929
930
931template <typename MatrixType>
934{
935 return *matrices.read<const MGLevelObject<MatrixType> *>(i);
936}
937
938
939template <typename MatrixType>
945
946
947template <typename MatrixType>
950{
951 return *matrices_in.read<const MGLevelObject<MatrixType> *>(i);
952}
953
954
955template <typename MatrixType>
958{
959 return *matrices_in.entry<MGLevelObject<MatrixType> *>(i);
960}
961
962
963template <typename MatrixType>
966{
967 return *matrices_out.read<const MGLevelObject<MatrixType> *>(i);
968}
969
970
971template <typename MatrixType>
974{
975 return *matrices_out.entry<MGLevelObject<MatrixType> *>(i);
976}
977
978
979template <typename MatrixType>
982{
983 return *flux_matrices_up.read<const MGLevelObject<MatrixType> *>(i);
984}
985
986
987template <typename MatrixType>
990{
991 return *flux_matrices_up.entry<MGLevelObject<MatrixType> *>(i);
992}
993
994
995template <typename MatrixType>
998{
999 return *flux_matrices_down.read<const MGLevelObject<MatrixType> *>(i);
1000}
1001
1002
1003template <typename MatrixType>
1006{
1007 return *flux_matrices_down.entry<MGLevelObject<MatrixType> *>(i);
1008}
1009
1010
1011template <typename MatrixType>
1012inline void
1015{
1016 for (size_type i = 0; i < this->size(); ++i)
1017 {
1019 const size_type row = o[o.min_level()].row;
1020 const size_type col = o[o.min_level()].column;
1021
1022 o.resize(sparsity.min_level(), sparsity.max_level());
1023 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1024 {
1025 o[level].row = row;
1026 o[level].column = col;
1027 internal::reinit(o[level], sparsity[level]);
1028 }
1029 }
1030}
1031
1032
1033template <typename MatrixType>
1034inline void
1037{
1038 for (size_type i = 0; i < this->size(); ++i)
1039 {
1041 const size_type row = o[o.min_level()].row;
1042 const size_type col = o[o.min_level()].column;
1043
1044 block_in(i).resize(sparsity.min_level(), sparsity.max_level());
1045 block_out(i).resize(sparsity.min_level(), sparsity.max_level());
1046 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1047 {
1048 block_in(i)[level].row = row;
1049 block_in(i)[level].column = col;
1050 internal::reinit(block_in(i)[level], sparsity[level]);
1051 block_out(i)[level].row = row;
1052 block_out(i)[level].column = col;
1053 internal::reinit(block_out(i)[level], sparsity[level]);
1054 }
1055 }
1056}
1057
1058
1059template <typename MatrixType>
1060inline void
1063{
1064 for (size_type i = 0; i < this->size(); ++i)
1065 {
1067 const size_type row = o[o.min_level()].row;
1068 const size_type col = o[o.min_level()].column;
1069
1070 block_up(i).resize(sparsity.min_level(), sparsity.max_level());
1071 block_down(i).resize(sparsity.min_level(), sparsity.max_level());
1072 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1073 {
1074 block_up(i)[level].row = row;
1075 block_up(i)[level].column = col;
1076 internal::reinit(block_up(i)[level], sparsity[level]);
1077 block_down(i)[level].row = row;
1078 block_down(i)[level].column = col;
1079 internal::reinit(block_down(i)[level], sparsity[level]);
1080 }
1081 }
1082}
1083
1084
1085template <typename MatrixType>
1086inline void
1088{
1089 for (size_type i = 0; i < mo.size(); ++i)
1090 {
1093 for (size_type level = o.min_level(); level <= o.max_level(); ++level)
1094 o[level].matrix.clear();
1095 }
1096}
1097
1098
1099template <typename MatrixType>
1100inline void
1102{
1103 if (really_clean)
1104 {
1106 }
1107 else
1108 {
1109 clear_object(matrices);
1110 clear_object(matrices_in);
1111 clear_object(matrices_out);
1112 clear_object(flux_matrices_up);
1113 clear_object(flux_matrices_down);
1114 }
1115}
1116
1117
1118
1120
1121#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 subscribe(std::atomic< bool > *const validity, const std::string &identifier="") const
void unsubscribe(std::atomic< bool > *const validity, const std::string &identifier="") 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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
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:534
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcNotInitialized()
static ::ExceptionBase & ExcNotQuadratic()
#define DEAL_II_NOT_IMPLEMENTED()
std::size_t size
Definition mpi.cc:734
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