Reference documentation for deal.II version Git ede8f93e86 2020-12-03 14:59:20 -0700
\(\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 - 2020 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
36 template <typename number>
37 class BlockSparseMatrix;
39 #endif
40 
77 template <typename SparsityPatternType>
79 {
80 public:
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
155  collect_sizes();
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 &
176  get_row_indices() const;
177 
182  const BlockIndices &
183  get_column_indices() const;
184 
189  void
190  compress();
191 
195  size_type
196  n_block_rows() const;
197 
201  size_type
202  n_block_cols() const;
203 
210  bool
211  empty() const;
212 
218  size_type
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
245  add_entries(const size_type row,
246  ForwardIterator begin,
247  ForwardIterator end,
248  const bool indices_are_sorted = false);
249 
254  size_type
255  n_rows() const;
256 
262  size_type
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 
289  size_type
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 
345 protected:
350 
355 
359  Table<2,
360  SmartPointer<SparsityPatternType,
363 
369 
375 
376 private:
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 
406 class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
407 {
408 public:
414  BlockSparsityPattern() = default;
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
467  copy_from(const BlockDynamicSparsityPattern &dsp);
468 };
469 
470 
471 
525  : public BlockSparsityPatternBase<DynamicSparsityPattern>
526 {
527 public:
533  BlockDynamicSparsityPattern() = default;
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 
601  size_type
602  column_number(const size_type row, const unsigned int index) const;
603 
608 };
609 
613 #ifdef DEAL_II_WITH_TRILINOS
614 
615 
616 namespace TrilinosWrappers
617 {
640  : public ::BlockSparsityPatternBase<SparsityPattern>
641  {
642  public:
648  BlockSparsityPattern() = default;
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 
754 template <typename SparsityPatternType>
755 inline 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 
766 template <typename SparsityPatternType>
767 inline 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 
779 template <typename SparsityPatternType>
780 inline const BlockIndices &
782 {
783  return row_indices;
784 }
785 
786 
787 
788 template <typename SparsityPatternType>
789 inline const BlockIndices &
791 {
792  return column_indices;
793 }
794 
795 
796 
797 template <typename SparsityPatternType>
798 inline 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 =
807  col_index =
809  sub_objects[row_index.first][col_index.first]->add(row_index.second,
810  col_index.second);
811 }
812 
813 
814 
815 template <typename SparsityPatternType>
816 template <typename ForwardIterator>
817 void
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 
904 template <typename SparsityPatternType>
905 inline 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 =
914  col_index =
916  return sub_objects[row_index.first][col_index.first]->exists(
917  row_index.second, col_index.second);
918 }
919 
920 
921 
922 template <typename SparsityPatternType>
923 inline unsigned int
925  const size_type row) const
926 {
927  const std::pair<size_type, size_type> row_index =
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 
940 template <typename SparsityPatternType>
943 {
944  return columns;
945 }
946 
947 
948 
949 template <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 
986 inline void
988  const size_type n_block_columns)
989 {
991  n_block_columns);
992 }
993 
994 
996 
997 #endif
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
size_type column_number(const size_type row, const unsigned int index) const
std::vector< std::vector< size_type > > block_column_indices
std::vector< size_type > counter_within_block
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
SparsityPatternType & block(const size_type row, const size_type column)
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 size_type n_block_rows, const size_type n_block_columns)
static const size_type invalid_entry
void add(const size_type i, const size_type j)
void print(std::ostream &out) const
#define Assert(cond, exc)
Definition: exceptions.h:1466
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
const BlockIndices & get_row_indices() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
VectorType::value_type * end(VectorType &V)
BlockSparsityPatternBase & operator=(const BlockSparsityPatternBase &)
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void print_svg(std::ostream &out) const
unsigned int global_dof_index
Definition: types.h:76
Table< 2, SmartPointer< SparsityPatternType, BlockSparsityPatternBase< SparsityPatternType > > > sub_objects
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:584
unsigned int row_length(const size_type row) const
void print_gnuplot(std::ostream &out) const
Definition: table.h:687
static const size_type invalid_entry
types::global_dof_index size_type
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
const BlockIndices & get_column_indices() const
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
static ::ExceptionBase & ExcInternalError()