Reference documentation for deal.II version GIT relicensing-426-g7976cfd195 2024-04-18 21:10:01+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
block_sparsity_pattern.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 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_block_sparsity_pattern_h
16#define dealii_block_sparsity_pattern_h
17
18
19#include <deal.II/base/config.h>
20
24#include <deal.II/base/table.h>
25
31
33
34// Forward declarations
35#ifndef DOXYGEN
36template <typename number>
39#endif
40
78template <typename SparsityPatternType>
80{
81public:
86
96
103
110 const size_type n_block_columns);
111
119
133 void
134 reinit(const size_type n_block_rows, const size_type n_block_columns);
135
143
150 void
152
156 SparsityPatternType &
157 block(const size_type row, const size_type column);
158
163 const SparsityPatternType &
164 block(const size_type row, const size_type column) const;
165
170 const BlockIndices &
172
177 const BlockIndices &
179
184 void
185 compress();
186
192
198
205 bool
206 empty() const;
207
214 max_entries_per_row() const;
215
225 void
226 add(const size_type i, const size_type j);
227
238 template <typename ForwardIterator>
239 void
241 ForwardIterator begin,
242 ForwardIterator end,
243 const bool indices_are_sorted = false);
244
250 virtual void
252 const ArrayView<const size_type> &columns,
253 const bool indices_are_sorted = false) override;
254
256
262
269
273 bool
274 exists(const size_type i, const size_type j) const;
275
280 unsigned int
281 row_length(const size_type row) const;
282
295 n_nonzero_elements() const;
296
302 void
303 print(std::ostream &out) const;
304
312 void
313 print_gnuplot(std::ostream &out) const;
314
320 void
321 print_svg(std::ostream &out) const;
322
327 std::size_t
328 memory_consumption() const;
329
340 "The number of rows and columns (returned by n_rows() and n_cols()) does "
341 "not match their directly computed values. This typically means that a "
342 "call to collect_sizes() is missing.");
343
348 int,
349 int,
350 int,
351 int,
352 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
353 << ',' << arg4 << "] have differing row numbers.");
358 int,
359 int,
360 int,
361 int,
362 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
363 << ',' << arg4 << "] have differing column numbers.");
366protected:
371
376
381
387
393
394private:
399 compute_n_rows() const;
400
405 compute_n_cols() const;
406
411 std::vector<size_type> counter_within_block;
412
417 std::vector<std::vector<size_type>> block_column_indices;
418
419 // Make the block sparse matrix a friend, so that it can use our
420 // #row_indices and #column_indices objects.
421 template <typename number>
422 friend class BlockSparseMatrix;
423};
424
425
426
436class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
437{
438public:
445
451 BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
452
456 void
457 reinit(const size_type n_block_rows, const size_type n_block_columns);
458
471 void
473 const BlockIndices &col_indices,
474 const std::vector<std::vector<unsigned int>> &row_lengths);
475
476
481 bool
482 is_compressed() const;
483
489 void
491};
492
493
494
548 : public BlockSparsityPatternBase<DynamicSparsityPattern>
549{
550public:
557
564 const size_type n_columns);
565
572 BlockDynamicSparsityPattern(const std::vector<size_type> &row_block_sizes,
573 const std::vector<size_type> &col_block_sizes);
574
583 BlockDynamicSparsityPattern(const std::vector<IndexSet> &partitioning);
584
590 const BlockIndices &col_indices);
591
592
601 void
602 reinit(const std::vector<size_type> &row_block_sizes,
603 const std::vector<size_type> &col_block_sizes);
604
609 void
610 reinit(const std::vector<IndexSet> &partitioning);
611
617 void
618 reinit(const BlockIndices &row_indices, const BlockIndices &col_indices);
619
625 column_number(const size_type row, const unsigned int index) const;
626
631};
632
636#ifdef DEAL_II_WITH_TRILINOS
637
638
639namespace TrilinosWrappers
640{
664 : public ::BlockSparsityPatternBase<SparsityPattern>
665 {
666 public:
673
679 BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
680
687 BlockSparsityPattern(const std::vector<size_type> &row_block_sizes,
688 const std::vector<size_type> &col_block_sizes);
689
697 BlockSparsityPattern(const std::vector<IndexSet> &parallel_partitioning,
698 const MPI_Comm communicator = MPI_COMM_WORLD);
699
712 const std::vector<IndexSet> &row_parallel_partitioning,
713 const std::vector<IndexSet> &column_parallel_partitioning,
714 const std::vector<IndexSet> &writeable_rows,
715 const MPI_Comm communicator = MPI_COMM_WORLD);
716
726 void
727 reinit(const std::vector<size_type> &row_block_sizes,
728 const std::vector<size_type> &col_block_sizes);
729
734 void
735 reinit(const std::vector<IndexSet> &parallel_partitioning,
736 const MPI_Comm communicator = MPI_COMM_WORLD);
737
743 void
744 reinit(const std::vector<IndexSet> &row_parallel_partitioning,
745 const std::vector<IndexSet> &column_parallel_partitioning,
746 const MPI_Comm communicator = MPI_COMM_WORLD);
747
756 void
757 reinit(const std::vector<IndexSet> &row_parallel_partitioning,
758 const std::vector<IndexSet> &column_parallel_partitioning,
759 const std::vector<IndexSet> &writeable_rows,
760 const MPI_Comm communicator = MPI_COMM_WORLD);
761
766 };
767
770} /* namespace TrilinosWrappers */
771
772#endif
773
774/*--------------------- Template functions ----------------------------------*/
775
776
777
778template <typename SparsityPatternType>
779inline SparsityPatternType &
781 const size_type column)
782{
783 AssertIndexRange(row, n_block_rows());
784 AssertIndexRange(column, n_block_cols());
785 return *sub_objects(row, column);
786}
787
788
789
790template <typename SparsityPatternType>
791inline const SparsityPatternType &
793 const size_type row,
794 const size_type column) const
795{
796 AssertIndexRange(row, n_block_rows());
797 AssertIndexRange(column, n_block_cols());
798 return *sub_objects(row, column);
799}
800
801
802
803template <typename SparsityPatternType>
804inline const BlockIndices &
809
810
811
812template <typename SparsityPatternType>
813inline const BlockIndices &
815{
816 return column_indices;
817}
818
819
820
821template <typename SparsityPatternType>
822inline void
824 const size_type j)
825{
826 // if you get an error here, are
827 // you sure you called
828 // <tt>collect_sizes()</tt> before?
829 const std::pair<size_type, size_type> row_index =
830 row_indices.global_to_local(i),
831 col_index =
832 column_indices.global_to_local(j);
833 sub_objects[row_index.first][col_index.first]->add(row_index.second,
834 col_index.second);
835}
836
837
838
839template <typename SparsityPatternType>
840template <typename ForwardIterator>
841void
843 const size_type row,
844 ForwardIterator begin,
845 ForwardIterator end,
846 const bool indices_are_sorted)
847{
848 // In debug mode, verify that collect_sizes() was called by redoing the
849 // calculation
850 Assert(n_rows() == compute_n_rows(), ExcNeedsCollectSizes());
851 Assert(n_cols() == compute_n_cols(), ExcNeedsCollectSizes());
852
853 const size_type n_cols = static_cast<size_type>(end - begin);
854
855 if (indices_are_sorted && n_cols > 0)
856 {
857 block_column_indices[0].resize(0);
858
859 const std::pair<size_type, size_type> row_index =
860 this->row_indices.global_to_local(row);
861 const auto n_blocks = column_indices.size();
862
863 // Assume we start with the first block: since we assemble sparsity
864 // patterns one cell at a time this should always be true
865 size_type current_block = 0;
866 size_type current_start_index = column_indices.block_start(current_block);
867 size_type next_start_index =
868 current_block == n_blocks - 1 ?
870 column_indices.block_start(current_block + 1);
871
872 for (auto it = begin; it < end; ++it)
873 {
874 // BlockIndices::global_to_local() is a bit slow so instead we just
875 // keep track of which block we are in - as the indices are sorted we
876 // know that the block number can only increase.
877 if (*it >= next_start_index)
878 {
879 // we found a column outside the present block: write the present
880 // block entries and continue to the next block
881 sub_objects[row_index.first][current_block]->add_entries(
882 row_index.second,
883 block_column_indices[0].begin(),
884 block_column_indices[0].end(),
885 true);
886 block_column_indices[0].clear();
887
888 auto block_and_col = column_indices.global_to_local(*it);
889 current_block = block_and_col.first;
890 current_start_index = column_indices.block_start(current_block);
891 next_start_index =
892 current_block == n_blocks - 1 ?
894 column_indices.block_start(current_block + 1);
895 }
896 const size_type local_index = *it - current_start_index;
897 block_column_indices[0].push_back(local_index);
898
899 // Check that calculation:
900#ifdef DEBUG
901 {
902 auto check_block_and_col = column_indices.global_to_local(*it);
903 Assert(current_block == check_block_and_col.first,
905 Assert(local_index == check_block_and_col.second,
907 }
908#endif
909 }
910 // add whatever is left over:
911 sub_objects[row_index.first][current_block]->add_entries(
912 row_index.second,
913 block_column_indices[0].begin(),
914 block_column_indices[0].end(),
915 true);
916
917 return;
918 }
919 else
920 {
921 // Resize sub-arrays to n_cols. This
922 // is a bit wasteful, but we resize
923 // only a few times (then the maximum
924 // row length won't increase that
925 // much any more). At least we know
926 // that all arrays are going to be of
927 // the same size, so we can check
928 // whether the size of one is large
929 // enough before actually going
930 // through all of them.
931 if (block_column_indices[0].size() < n_cols)
932 for (size_type i = 0; i < this->n_block_cols(); ++i)
933 block_column_indices[i].resize(n_cols);
934
935 // Reset the number of added elements
936 // in each block to zero.
937 for (size_type i = 0; i < this->n_block_cols(); ++i)
938 counter_within_block[i] = 0;
939
940 // Go through the column indices to
941 // find out which portions of the
942 // values should be set in which
943 // block of the matrix. We need to
944 // touch all the data, since we can't
945 // be sure that the data of one block
946 // is stored contiguously (in fact,
947 // indices will be intermixed when it
948 // comes from an element matrix).
949 for (ForwardIterator it = begin; it != end; ++it)
950 {
951 const size_type col = *it;
952
953 const std::pair<size_type, size_type> col_index =
954 this->column_indices.global_to_local(col);
955
956 const size_type local_index = counter_within_block[col_index.first]++;
957
958 block_column_indices[col_index.first][local_index] = col_index.second;
959 }
960
961 // Now we found out about where the
962 // individual columns should start and
963 // where we should start reading out
964 // data. Now let's write the data into
965 // the individual blocks!
966 const std::pair<size_type, size_type> row_index =
967 this->row_indices.global_to_local(row);
968 for (size_type block_col = 0; block_col < n_block_cols(); ++block_col)
969 {
970 if (counter_within_block[block_col] == 0)
971 continue;
972 sub_objects[row_index.first][block_col]->add_entries(
973 row_index.second,
974 block_column_indices[block_col].begin(),
975 block_column_indices[block_col].begin() +
976 counter_within_block[block_col],
977 indices_are_sorted);
978 }
979 }
980}
981
982
983
984template <typename SparsityPatternType>
985void
987 const size_type &row,
988 const ArrayView<const size_type> &columns,
989 const bool indices_are_sorted)
990{
991 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
992}
993
994
995
996template <typename SparsityPatternType>
997inline bool
999 const size_type j) const
1000{
1001 // if you get an error here, are
1002 // you sure you called
1003 // <tt>collect_sizes()</tt> before?
1004 const std::pair<size_type, size_type> row_index =
1005 row_indices.global_to_local(i),
1006 col_index =
1007 column_indices.global_to_local(j);
1008 return sub_objects[row_index.first][col_index.first]->exists(
1009 row_index.second, col_index.second);
1010}
1011
1012
1013
1014template <typename SparsityPatternType>
1015inline unsigned int
1017 const size_type row) const
1018{
1019 const std::pair<size_type, size_type> row_index =
1020 row_indices.global_to_local(row);
1021
1022 unsigned int c = 0;
1023
1024 for (size_type b = 0; b < n_block_rows(); ++b)
1025 c += sub_objects[row_index.first][b]->row_length(row_index.second);
1026
1027 return c;
1028}
1029
1030
1031
1032template <typename SparsityPatternType>
1035{
1036 return block_columns;
1037}
1038
1039
1040
1041template <typename SparsityPatternType>
1044{
1045 return block_rows;
1046}
1047
1048
1051 const unsigned int index) const
1052{
1053 // .first= ith block, .second = jth row in that block
1054 const std::pair<size_type, size_type> row_index =
1056
1057 AssertIndexRange(index, row_length(row));
1058
1059 size_type c = 0;
1060 size_type block_columns = 0; // sum of n_cols for all blocks to the left
1061 for (unsigned int b = 0; b < this->n_block_cols(); ++b)
1062 {
1063 unsigned int rowlen =
1064 sub_objects[row_index.first][b]->row_length(row_index.second);
1065 if (index < c + rowlen)
1066 return block_columns +
1067 sub_objects[row_index.first][b]->column_number(row_index.second,
1068 index - c);
1069 c += rowlen;
1070 block_columns += sub_objects[row_index.first][b]->n_cols();
1071 }
1072
1074 return 0;
1075}
1076
1077
1078inline void
1080 const size_type new_block_columns)
1081{
1083 new_block_columns);
1084}
1085
1086
1088
1089#endif
iterator begin() const
Definition array_view.h:702
iterator end() const
Definition array_view.h:711
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
size_type column_number(const size_type row, const unsigned int index) const
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
void print_gnuplot(std::ostream &out) const
static const size_type invalid_entry
const SparsityPatternType & block(const size_type row, const size_type column) const
std::vector< size_type > counter_within_block
Table< 2, std::unique_ptr< SparsityPatternType > > sub_objects
SparsityPatternType & block(const size_type row, const size_type column)
std::vector< std::vector< size_type > > block_column_indices
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
const BlockIndices & get_column_indices() const
void print(std::ostream &out) const
std::size_t memory_consumption() const
void print_svg(std::ostream &out) const
void reinit(const size_type n_block_rows, const size_type n_block_columns)
const BlockIndices & get_row_indices() const
unsigned int row_length(const size_type row) const
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
void add(const size_type i, const size_type j)
bool exists(const size_type i, const size_type j) const
BlockSparsityPatternBase & operator=(const BlockSparsityPatternBase &)
void copy_from(const BlockDynamicSparsityPattern &dsp)
void reinit(const size_type n_block_rows, const size_type n_block_columns)
BlockSparsityPattern()=default
size_type n_rows() const
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type > > &entries)
size_type n_cols() const
static constexpr size_type invalid_entry
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNeedsCollectSizes()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:585
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
#define DeclExceptionMsg(Exception, defaulttext)
Definition exceptions.h:494
static ::ExceptionBase & ExcInternalError()
#define DEAL_II_ASSERT_UNREACHABLE()
const types::global_dof_index invalid_dof_index
Definition types.h:252
unsigned int global_dof_index
Definition types.h:81