Reference documentation for deal.II version GIT 20f059c89a 2022-12-04 14:55:02+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\}}\)
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 
32 
34 
35 // Forward declarations
36 #ifndef DOXYGEN
37 template <typename number>
38 class BlockSparseMatrix;
40 #endif
41 
79 template <typename SparsityPatternType>
81 {
82 public:
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
152  collect_sizes();
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 
191  size_type
192  n_block_rows() const;
193 
197  size_type
198  n_block_cols() const;
199 
206  bool
207  empty() const;
208 
214  size_type
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 
295  size_type
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.");
367 protected:
372 
377 
382 
388 
394 
395 private:
399  size_type
400  compute_n_rows() const;
401 
405  size_type
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 
437 class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
438 {
439 public:
445  BlockSparsityPattern() = default;
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 {
551 public:
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 
625  size_type
626  column_number(const size_type row, const unsigned int index) const;
627 
632 };
633 
637 #ifdef DEAL_II_WITH_TRILINOS
638 
639 
640 namespace TrilinosWrappers
641 {
665  : public ::BlockSparsityPatternBase<SparsityPattern>
666  {
667  public:
673  BlockSparsityPattern() = default;
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 
779 template <typename SparsityPatternType>
780 inline 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 
791 template <typename SparsityPatternType>
792 inline 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 
804 template <typename SparsityPatternType>
805 inline const BlockIndices &
807 {
808  return row_indices;
809 }
810 
811 
812 
813 template <typename SparsityPatternType>
814 inline const BlockIndices &
816 {
817  return column_indices;
818 }
819 
820 
821 
822 template <typename SparsityPatternType>
823 inline 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 
840 template <typename SparsityPatternType>
841 template <typename ForwardIterator>
842 void
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,
905  ExcInternalError());
906  Assert(local_index == check_block_and_col.second,
907  ExcInternalError());
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 
985 template <typename SparsityPatternType>
986 void
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 
997 template <typename SparsityPatternType>
998 inline 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 
1015 template <typename SparsityPatternType>
1016 inline 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 
1033 template <typename SparsityPatternType>
1036 {
1037  return block_columns;
1038 }
1039 
1040 
1041 
1042 template <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 
1079 inline 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:585
iterator end() const
Definition: array_view.h:594
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
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
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
types::global_dof_index size_type
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:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
static ::ExceptionBase & ExcInternalError()
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:578
#define Assert(cond, exc)
Definition: exceptions.h:1501
#define AssertIndexRange(index, range)
Definition: exceptions.h:1760
#define DeclExceptionMsg(Exception, defaulttext)
Definition: exceptions.h:487
static ::ExceptionBase & ExcNeedsCollectSizes()
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
const types::global_dof_index invalid_dof_index
Definition: types.h:226
unsigned int global_dof_index
Definition: types.h:81