Reference documentation for deal.II version 9.4.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\}}\)
block_sparsity_pattern.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2000 - 2022 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
31
33
34// Forward declarations
35#ifndef DOXYGEN
36template <typename number>
39#endif
40
77template <typename SparsityPatternType>
79{
80public:
85
95
102
109 const size_type n_block_columns);
110
118
122 ~BlockSparsityPatternBase() override;
123
137 void
138 reinit(const size_type n_block_rows, const size_type n_block_columns);
139
147
154 void
156
160 SparsityPatternType &
161 block(const size_type row, const size_type column);
162
163
168 const SparsityPatternType &
169 block(const size_type row, const size_type column) const;
170
175 const BlockIndices &
177
182 const BlockIndices &
184
189 void
190 compress();
191
197
203
210 bool
211 empty() const;
212
219 max_entries_per_row() const;
220
230 void
231 add(const size_type i, const size_type j);
232
243 template <typename ForwardIterator>
244 void
246 ForwardIterator begin,
247 ForwardIterator end,
248 const bool indices_are_sorted = false);
249
255 n_rows() const;
256
263 n_cols() const;
264
268 bool
269 exists(const size_type i, const size_type j) const;
270
275 unsigned int
276 row_length(const size_type row) const;
277
290 n_nonzero_elements() const;
291
297 void
298 print(std::ostream &out) const;
299
307 void
308 print_gnuplot(std::ostream &out) const;
309
315 void
316 print_svg(std::ostream &out) const;
317
327 int,
328 int,
329 int,
330 int,
331 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
332 << ',' << arg4 << "] have differing row numbers.");
337 int,
338 int,
339 int,
340 int,
341 << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
342 << ',' << arg4 << "] have differing column numbers.");
344
345protected:
350
355
359 Table<2,
360 SmartPointer<SparsityPatternType,
363
369
375
376private:
381 std::vector<size_type> counter_within_block;
382
387 std::vector<std::vector<size_type>> block_column_indices;
388
389 // Make the block sparse matrix a friend, so that it can use our
390 // #row_indices and #column_indices objects.
391 template <typename number>
392 friend class BlockSparseMatrix;
393};
394
395
396
406class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
407{
408public:
415
421 BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
422
426 void
427 reinit(const size_type n_block_rows, const size_type n_block_columns);
428
441 void
443 const BlockIndices & col_indices,
444 const std::vector<std::vector<unsigned int>> &row_lengths);
445
446
451 bool
452 is_compressed() const;
453
458 std::size_t
459 memory_consumption() const;
460
466 void
468};
469
470
471
525 : public BlockSparsityPatternBase<DynamicSparsityPattern>
526{
527public:
534
541 const size_type n_columns);
542
549 BlockDynamicSparsityPattern(const std::vector<size_type> &row_block_sizes,
550 const std::vector<size_type> &col_block_sizes);
551
560 BlockDynamicSparsityPattern(const std::vector<IndexSet> &partitioning);
561
567 const BlockIndices &col_indices);
568
569
578 void
579 reinit(const std::vector<size_type> &row_block_sizes,
580 const std::vector<size_type> &col_block_sizes);
581
586 void
587 reinit(const std::vector<IndexSet> &partitioning);
588
594 void
595 reinit(const BlockIndices &row_indices, const BlockIndices &col_indices);
596
602 column_number(const size_type row, const unsigned int index) const;
603
608};
609
613#ifdef DEAL_II_WITH_TRILINOS
614
615
616namespace TrilinosWrappers
617{
640 : public ::BlockSparsityPatternBase<SparsityPattern>
641 {
642 public:
649
655 BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
656
663 BlockSparsityPattern(const std::vector<size_type> &row_block_sizes,
664 const std::vector<size_type> &col_block_sizes);
665
673 BlockSparsityPattern(const std::vector<IndexSet> &parallel_partitioning,
674 const MPI_Comm &communicator = MPI_COMM_WORLD);
675
688 const std::vector<IndexSet> &row_parallel_partitioning,
689 const std::vector<IndexSet> &column_parallel_partitioning,
690 const std::vector<IndexSet> &writeable_rows,
691 const MPI_Comm & communicator = MPI_COMM_WORLD);
692
702 void
703 reinit(const std::vector<size_type> &row_block_sizes,
704 const std::vector<size_type> &col_block_sizes);
705
710 void
711 reinit(const std::vector<IndexSet> &parallel_partitioning,
712 const MPI_Comm & communicator = MPI_COMM_WORLD);
713
719 void
720 reinit(const std::vector<IndexSet> &row_parallel_partitioning,
721 const std::vector<IndexSet> &column_parallel_partitioning,
722 const MPI_Comm & communicator = MPI_COMM_WORLD);
723
732 void
733 reinit(const std::vector<IndexSet> &row_parallel_partitioning,
734 const std::vector<IndexSet> &column_parallel_partitioning,
735 const std::vector<IndexSet> &writeable_rows,
736 const MPI_Comm & communicator = MPI_COMM_WORLD);
737
742 };
743
746} /* namespace TrilinosWrappers */
747
748#endif
749
750/*--------------------- Template functions ----------------------------------*/
751
752
753
754template <typename SparsityPatternType>
755inline SparsityPatternType &
757 const size_type column)
758{
759 AssertIndexRange(row, rows);
760 AssertIndexRange(column, columns);
761 return *sub_objects[row][column];
762}
763
764
765
766template <typename SparsityPatternType>
767inline const SparsityPatternType &
769 const size_type row,
770 const size_type column) const
771{
772 AssertIndexRange(row, rows);
773 AssertIndexRange(column, columns);
774 return *sub_objects[row][column];
775}
776
777
778
779template <typename SparsityPatternType>
780inline const BlockIndices &
782{
783 return row_indices;
784}
785
786
787
788template <typename SparsityPatternType>
789inline const BlockIndices &
791{
792 return column_indices;
793}
794
795
796
797template <typename SparsityPatternType>
798inline void
800 const size_type j)
801{
802 // if you get an error here, are
803 // you sure you called
804 // <tt>collect_sizes()</tt> before?
805 const std::pair<size_type, size_type> row_index =
806 row_indices.global_to_local(i),
807 col_index =
808 column_indices.global_to_local(j);
809 sub_objects[row_index.first][col_index.first]->add(row_index.second,
810 col_index.second);
811}
812
813
814
815template <typename SparsityPatternType>
816template <typename ForwardIterator>
817void
819 const size_type row,
820 ForwardIterator begin,
821 ForwardIterator end,
822 const bool indices_are_sorted)
823{
824 // Resize scratch arrays
825 if (block_column_indices.size() < this->n_block_cols())
826 {
827 block_column_indices.resize(this->n_block_cols());
828 counter_within_block.resize(this->n_block_cols());
829 }
830
831 const size_type n_cols = static_cast<size_type>(end - begin);
832
833 // Resize sub-arrays to n_cols. This
834 // is a bit wasteful, but we resize
835 // only a few times (then the maximum
836 // row length won't increase that
837 // much any more). At least we know
838 // that all arrays are going to be of
839 // the same size, so we can check
840 // whether the size of one is large
841 // enough before actually going
842 // through all of them.
843 if (block_column_indices[0].size() < n_cols)
844 for (size_type i = 0; i < this->n_block_cols(); ++i)
845 block_column_indices[i].resize(n_cols);
846
847 // Reset the number of added elements
848 // in each block to zero.
849 for (size_type i = 0; i < this->n_block_cols(); ++i)
850 counter_within_block[i] = 0;
851
852 // Go through the column indices to
853 // find out which portions of the
854 // values should be set in which
855 // block of the matrix. We need to
856 // touch all the data, since we can't
857 // be sure that the data of one block
858 // is stored contiguously (in fact,
859 // indices will be intermixed when it
860 // comes from an element matrix).
861 for (ForwardIterator it = begin; it != end; ++it)
862 {
863 const size_type col = *it;
864
865 const std::pair<size_type, size_type> col_index =
866 this->column_indices.global_to_local(col);
867
868 const size_type local_index = counter_within_block[col_index.first]++;
869
870 block_column_indices[col_index.first][local_index] = col_index.second;
871 }
872
873#ifdef DEBUG
874 // If in debug mode, do a check whether
875 // the right length has been obtained.
876 size_type length = 0;
877 for (size_type i = 0; i < this->n_block_cols(); ++i)
878 length += counter_within_block[i];
879 Assert(length == n_cols, ExcInternalError());
880#endif
881
882 // Now we found out about where the
883 // individual columns should start and
884 // where we should start reading out
885 // data. Now let's write the data into
886 // the individual blocks!
887 const std::pair<size_type, size_type> row_index =
888 this->row_indices.global_to_local(row);
889 for (size_type block_col = 0; block_col < n_block_cols(); ++block_col)
890 {
891 if (counter_within_block[block_col] == 0)
892 continue;
893 sub_objects[row_index.first][block_col]->add_entries(
894 row_index.second,
895 block_column_indices[block_col].begin(),
896 block_column_indices[block_col].begin() +
897 counter_within_block[block_col],
898 indices_are_sorted);
899 }
900}
901
902
903
904template <typename SparsityPatternType>
905inline bool
907 const size_type j) const
908{
909 // if you get an error here, are
910 // you sure you called
911 // <tt>collect_sizes()</tt> before?
912 const std::pair<size_type, size_type> row_index =
913 row_indices.global_to_local(i),
914 col_index =
915 column_indices.global_to_local(j);
916 return sub_objects[row_index.first][col_index.first]->exists(
917 row_index.second, col_index.second);
918}
919
920
921
922template <typename SparsityPatternType>
923inline unsigned int
925 const size_type row) const
926{
927 const std::pair<size_type, size_type> row_index =
928 row_indices.global_to_local(row);
929
930 unsigned int c = 0;
931
932 for (size_type b = 0; b < rows; ++b)
933 c += sub_objects[row_index.first][b]->row_length(row_index.second);
934
935 return c;
936}
937
938
939
940template <typename SparsityPatternType>
943{
944 return columns;
945}
946
947
948
949template <typename SparsityPatternType>
952{
953 return rows;
954}
955
956
959 const unsigned int index) const
960{
961 // .first= ith block, .second = jth row in that block
962 const std::pair<size_type, size_type> row_index =
964
965 AssertIndexRange(index, row_length(row));
966
967 size_type c = 0;
968 size_type block_columns = 0; // sum of n_cols for all blocks to the left
969 for (unsigned int b = 0; b < columns; ++b)
970 {
971 unsigned int rowlen =
972 sub_objects[row_index.first][b]->row_length(row_index.second);
973 if (index < c + rowlen)
974 return block_columns +
975 sub_objects[row_index.first][b]->column_number(row_index.second,
976 index - c);
977 c += rowlen;
978 block_columns += sub_objects[row_index.first][b]->n_cols();
979 }
980
981 Assert(false, ExcInternalError());
982 return 0;
983}
984
985
986inline void
988 const size_type n_block_columns)
989{
991 n_block_columns);
992}
993
994
996
997#endif
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
static const size_type invalid_entry
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
std::vector< size_type > counter_within_block
size_type column_number(const size_type row, const unsigned int index) const
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
std::vector< std::vector< size_type > > block_column_indices
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)
Definition: exceptions.h:1473
void copy_from(const BlockDynamicSparsityPattern &dsp)
Table< 2, SmartPointer< SparsityPatternType, BlockSparsityPatternBase< SparsityPatternType > > > sub_objects
void reinit(const size_type n_block_rows, const size_type n_block_columns)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
BlockSparsityPattern()=default
std::size_t memory_consumption() const
static const size_type invalid_entry
const SparsityPatternType & block(const size_type row, const size_type column) const
void print_gnuplot(std::ostream &out) const
void print(std::ostream &out) const
SparsityPatternType & block(const size_type row, const size_type column)
types::global_dof_index size_type
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
const BlockIndices & get_column_indices() const
const BlockIndices & get_row_indices() const
unsigned int row_length(const size_type row) const
void add(const size_type i, const size_type j)
BlockSparsityPatternBase & operator=(const BlockSparsityPatternBase &)
void print_svg(std::ostream &out) const
bool exists(const size_type i, const size_type j) const
void reinit(const size_type n_block_rows, const size_type n_block_columns)
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
unsigned int global_dof_index
Definition: types.h:76