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
block_sparsity_pattern.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 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_block_sparsity_pattern_h
17#define dealii_block_sparsity_pattern_h
18
19
20#include <deal.II/base/config.h>
21
25#include <deal.II/base/table.h>
26
32
34
35// Forward declarations
36#ifndef DOXYGEN
37template <typename number>
40#endif
41
79template <typename SparsityPatternType>
81{
82public:
87
97
104
111 const size_type n_block_columns);
112
120
134 void
135 reinit(const size_type n_block_rows, const size_type n_block_columns);
136
144
151 void
153
157 SparsityPatternType &
158 block(const size_type row, const size_type column);
159
164 const SparsityPatternType &
165 block(const size_type row, const size_type column) const;
166
171 const BlockIndices &
173
178 const BlockIndices &
180
185 void
186 compress();
187
193
199
206 bool
207 empty() const;
208
215 max_entries_per_row() const;
216
226 void
227 add(const size_type i, const size_type j);
228
239 template <typename ForwardIterator>
240 void
242 ForwardIterator begin,
243 ForwardIterator end,
244 const bool indices_are_sorted = false);
245
251 virtual void
253 const ArrayView<const size_type> &columns,
254 const bool indices_are_sorted = false) override;
255
257
263
270
274 bool
275 exists(const size_type i, const size_type j) const;
276
281 unsigned int
282 row_length(const size_type row) const;
283
296 n_nonzero_elements() const;
297
303 void
304 print(std::ostream &out) const;
305
313 void
314 print_gnuplot(std::ostream &out) const;
315
321 void
322 print_svg(std::ostream &out) const;
323
328 std::size_t
329 memory_consumption() const;
330
341 "The number of rows and columns (returned by n_rows() and n_cols()) does "
342 "not match their directly computed values. This typically means that a "
343 "call to collect_sizes() is missing.");
344
349 int,
350 int,
351 int,
352 int,
353 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
354 << ',' << arg4 << "] have differing row numbers.");
359 int,
360 int,
361 int,
362 int,
363 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
364 << ',' << arg4 << "] have differing column numbers.");
367protected:
372
377
382
388
394
395private:
400 compute_n_rows() const;
401
406 compute_n_cols() const;
407
412 std::vector<size_type> counter_within_block;
413
418 std::vector<std::vector<size_type>> block_column_indices;
419
420 // Make the block sparse matrix a friend, so that it can use our
421 // #row_indices and #column_indices objects.
422 template <typename number>
423 friend class BlockSparseMatrix;
424};
425
426
427
437class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
438{
439public:
446
452 BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
453
457 void
458 reinit(const size_type n_block_rows, const size_type n_block_columns);
459
472 void
474 const BlockIndices & col_indices,
475 const std::vector<std::vector<unsigned int>> &row_lengths);
476
477
482 bool
483 is_compressed() const;
484
490 void
492};
493
494
495
549 : public BlockSparsityPatternBase<DynamicSparsityPattern>
550{
551public:
558
565 const size_type n_columns);
566
573 BlockDynamicSparsityPattern(const std::vector<size_type> &row_block_sizes,
574 const std::vector<size_type> &col_block_sizes);
575
584 BlockDynamicSparsityPattern(const std::vector<IndexSet> &partitioning);
585
591 const BlockIndices &col_indices);
592
593
602 void
603 reinit(const std::vector<size_type> &row_block_sizes,
604 const std::vector<size_type> &col_block_sizes);
605
610 void
611 reinit(const std::vector<IndexSet> &partitioning);
612
618 void
619 reinit(const BlockIndices &row_indices, const BlockIndices &col_indices);
620
626 column_number(const size_type row, const unsigned int index) const;
627
632};
633
637#ifdef DEAL_II_WITH_TRILINOS
638
639
640namespace TrilinosWrappers
641{
665 : public ::BlockSparsityPatternBase<SparsityPattern>
666 {
667 public:
674
680 BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
681
688 BlockSparsityPattern(const std::vector<size_type> &row_block_sizes,
689 const std::vector<size_type> &col_block_sizes);
690
698 BlockSparsityPattern(const std::vector<IndexSet> &parallel_partitioning,
699 const MPI_Comm communicator = MPI_COMM_WORLD);
700
713 const std::vector<IndexSet> &row_parallel_partitioning,
714 const std::vector<IndexSet> &column_parallel_partitioning,
715 const std::vector<IndexSet> &writeable_rows,
716 const MPI_Comm communicator = MPI_COMM_WORLD);
717
727 void
728 reinit(const std::vector<size_type> &row_block_sizes,
729 const std::vector<size_type> &col_block_sizes);
730
735 void
736 reinit(const std::vector<IndexSet> &parallel_partitioning,
737 const MPI_Comm communicator = MPI_COMM_WORLD);
738
744 void
745 reinit(const std::vector<IndexSet> &row_parallel_partitioning,
746 const std::vector<IndexSet> &column_parallel_partitioning,
747 const MPI_Comm communicator = MPI_COMM_WORLD);
748
757 void
758 reinit(const std::vector<IndexSet> &row_parallel_partitioning,
759 const std::vector<IndexSet> &column_parallel_partitioning,
760 const std::vector<IndexSet> &writeable_rows,
761 const MPI_Comm communicator = MPI_COMM_WORLD);
762
767 };
768
771} /* namespace TrilinosWrappers */
772
773#endif
774
775/*--------------------- Template functions ----------------------------------*/
776
777
778
779template <typename SparsityPatternType>
780inline SparsityPatternType &
782 const size_type column)
783{
784 AssertIndexRange(row, n_block_rows());
785 AssertIndexRange(column, n_block_cols());
786 return *sub_objects(row, column);
787}
788
789
790
791template <typename SparsityPatternType>
792inline const SparsityPatternType &
794 const size_type row,
795 const size_type column) const
796{
797 AssertIndexRange(row, n_block_rows());
798 AssertIndexRange(column, n_block_cols());
799 return *sub_objects(row, column);
800}
801
802
803
804template <typename SparsityPatternType>
805inline const BlockIndices &
807{
808 return row_indices;
809}
810
811
812
813template <typename SparsityPatternType>
814inline const BlockIndices &
816{
817 return column_indices;
818}
819
820
821
822template <typename SparsityPatternType>
823inline void
825 const size_type j)
826{
827 // if you get an error here, are
828 // you sure you called
829 // <tt>collect_sizes()</tt> before?
830 const std::pair<size_type, size_type> row_index =
831 row_indices.global_to_local(i),
832 col_index =
833 column_indices.global_to_local(j);
834 sub_objects[row_index.first][col_index.first]->add(row_index.second,
835 col_index.second);
836}
837
838
839
840template <typename SparsityPatternType>
841template <typename ForwardIterator>
842void
844 const size_type row,
845 ForwardIterator begin,
846 ForwardIterator end,
847 const bool indices_are_sorted)
848{
849 // In debug mode, verify that collect_sizes() was called by redoing the
850 // calculation
851 Assert(n_rows() == compute_n_rows(), ExcNeedsCollectSizes());
852 Assert(n_cols() == compute_n_cols(), ExcNeedsCollectSizes());
853
854 const size_type n_cols = static_cast<size_type>(end - begin);
855
856 if (indices_are_sorted && n_cols > 0)
857 {
858 block_column_indices[0].resize(0);
859
860 const std::pair<size_type, size_type> row_index =
861 this->row_indices.global_to_local(row);
862 const auto n_blocks = column_indices.size();
863
864 // Assume we start with the first block: since we assemble sparsity
865 // patterns one cell at a time this should always be true
866 size_type current_block = 0;
867 size_type current_start_index = column_indices.block_start(current_block);
868 size_type next_start_index =
869 current_block == n_blocks - 1 ?
871 column_indices.block_start(current_block + 1);
872
873 for (auto it = begin; it < end; ++it)
874 {
875 // BlockIndices::global_to_local() is a bit slow so instead we just
876 // keep track of which block we are in - as the indices are sorted we
877 // know that the block number can only increase.
878 if (*it >= next_start_index)
879 {
880 // we found a column outside the present block: write the present
881 // block entries and continue to the next block
882 sub_objects[row_index.first][current_block]->add_entries(
883 row_index.second,
884 block_column_indices[0].begin(),
885 block_column_indices[0].end(),
886 true);
887 block_column_indices[0].clear();
888
889 auto block_and_col = column_indices.global_to_local(*it);
890 current_block = block_and_col.first;
891 current_start_index = column_indices.block_start(current_block);
892 next_start_index =
893 current_block == n_blocks - 1 ?
895 column_indices.block_start(current_block + 1);
896 }
897 const size_type local_index = *it - current_start_index;
898 block_column_indices[0].push_back(local_index);
899
900 // Check that calculation:
901#ifdef DEBUG
902 {
903 auto check_block_and_col = column_indices.global_to_local(*it);
904 Assert(current_block == check_block_and_col.first,
906 Assert(local_index == check_block_and_col.second,
908 }
909#endif
910 }
911 // add whatever is left over:
912 sub_objects[row_index.first][current_block]->add_entries(
913 row_index.second,
914 block_column_indices[0].begin(),
915 block_column_indices[0].end(),
916 true);
917
918 return;
919 }
920 else
921 {
922 // Resize sub-arrays to n_cols. This
923 // is a bit wasteful, but we resize
924 // only a few times (then the maximum
925 // row length won't increase that
926 // much any more). At least we know
927 // that all arrays are going to be of
928 // the same size, so we can check
929 // whether the size of one is large
930 // enough before actually going
931 // through all of them.
932 if (block_column_indices[0].size() < n_cols)
933 for (size_type i = 0; i < this->n_block_cols(); ++i)
934 block_column_indices[i].resize(n_cols);
935
936 // Reset the number of added elements
937 // in each block to zero.
938 for (size_type i = 0; i < this->n_block_cols(); ++i)
939 counter_within_block[i] = 0;
940
941 // Go through the column indices to
942 // find out which portions of the
943 // values should be set in which
944 // block of the matrix. We need to
945 // touch all the data, since we can't
946 // be sure that the data of one block
947 // is stored contiguously (in fact,
948 // indices will be intermixed when it
949 // comes from an element matrix).
950 for (ForwardIterator it = begin; it != end; ++it)
951 {
952 const size_type col = *it;
953
954 const std::pair<size_type, size_type> col_index =
955 this->column_indices.global_to_local(col);
956
957 const size_type local_index = counter_within_block[col_index.first]++;
958
959 block_column_indices[col_index.first][local_index] = col_index.second;
960 }
961
962 // Now we found out about where the
963 // individual columns should start and
964 // where we should start reading out
965 // data. Now let's write the data into
966 // the individual blocks!
967 const std::pair<size_type, size_type> row_index =
968 this->row_indices.global_to_local(row);
969 for (size_type block_col = 0; block_col < n_block_cols(); ++block_col)
970 {
971 if (counter_within_block[block_col] == 0)
972 continue;
973 sub_objects[row_index.first][block_col]->add_entries(
974 row_index.second,
975 block_column_indices[block_col].begin(),
976 block_column_indices[block_col].begin() +
977 counter_within_block[block_col],
978 indices_are_sorted);
979 }
980 }
981}
982
983
984
985template <typename SparsityPatternType>
986void
988 const size_type & row,
989 const ArrayView<const size_type> &columns,
990 const bool indices_are_sorted)
991{
992 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
993}
994
995
996
997template <typename SparsityPatternType>
998inline bool
1000 const size_type j) const
1001{
1002 // if you get an error here, are
1003 // you sure you called
1004 // <tt>collect_sizes()</tt> before?
1005 const std::pair<size_type, size_type> row_index =
1006 row_indices.global_to_local(i),
1007 col_index =
1008 column_indices.global_to_local(j);
1009 return sub_objects[row_index.first][col_index.first]->exists(
1010 row_index.second, col_index.second);
1011}
1012
1013
1014
1015template <typename SparsityPatternType>
1016inline unsigned int
1018 const size_type row) const
1019{
1020 const std::pair<size_type, size_type> row_index =
1021 row_indices.global_to_local(row);
1022
1023 unsigned int c = 0;
1024
1025 for (size_type b = 0; b < n_block_rows(); ++b)
1026 c += sub_objects[row_index.first][b]->row_length(row_index.second);
1027
1028 return c;
1029}
1030
1031
1032
1033template <typename SparsityPatternType>
1036{
1037 return block_columns;
1038}
1039
1040
1041
1042template <typename SparsityPatternType>
1045{
1046 return block_rows;
1047}
1048
1049
1052 const unsigned int index) const
1053{
1054 // .first= ith block, .second = jth row in that block
1055 const std::pair<size_type, size_type> row_index =
1057
1058 AssertIndexRange(index, row_length(row));
1059
1060 size_type c = 0;
1061 size_type block_columns = 0; // sum of n_cols for all blocks to the left
1062 for (unsigned int b = 0; b < this->n_block_cols(); ++b)
1063 {
1064 unsigned int rowlen =
1065 sub_objects[row_index.first][b]->row_length(row_index.second);
1066 if (index < c + rowlen)
1067 return block_columns +
1068 sub_objects[row_index.first][b]->column_number(row_index.second,
1069 index - c);
1070 c += rowlen;
1071 block_columns += sub_objects[row_index.first][b]->n_cols();
1072 }
1073
1074 Assert(false, ExcInternalError());
1075 return 0;
1076}
1077
1078
1079inline void
1081 const size_type new_block_columns)
1082{
1084 new_block_columns);
1085}
1086
1087
1089
1090#endif
iterator begin() const
Definition array_view.h:594
iterator end() const
Definition array_view.h:603
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)
types::global_dof_index size_type
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNeedsCollectSizes()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition exceptions.h:579
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:488
static ::ExceptionBase & ExcInternalError()
const types::global_dof_index invalid_dof_index
Definition types.h:233
unsigned int global_dof_index
Definition types.h:82