Reference documentation for deal.II version Git 92c1895bf8 2020-06-04 17:48:28 -0400
\(\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 
78 template <typename SparsityPatternType>
80 {
81 public:
86 
96 
103 
110  const size_type n_block_columns);
111 
119 
123  ~BlockSparsityPatternBase() override;
124 
138  void
139  reinit(const size_type n_block_rows, const size_type n_block_columns);
140 
148 
155  void
156  collect_sizes();
157 
161  SparsityPatternType &
162  block(const size_type row, const size_type column);
163 
164 
169  const SparsityPatternType &
170  block(const size_type row, const size_type column) const;
171 
176  const BlockIndices &
177  get_row_indices() const;
178 
183  const BlockIndices &
184  get_column_indices() const;
185 
190  void
191  compress();
192 
196  size_type
197  n_block_rows() const;
198 
202  size_type
203  n_block_cols() const;
204 
211  bool
212  empty() const;
213 
219  size_type
220  max_entries_per_row() const;
221 
231  void
232  add(const size_type i, const size_type j);
233 
244  template <typename ForwardIterator>
245  void
246  add_entries(const size_type row,
247  ForwardIterator begin,
248  ForwardIterator end,
249  const bool indices_are_sorted = false);
250 
255  size_type
256  n_rows() const;
257 
263  size_type
264  n_cols() const;
265 
269  bool
270  exists(const size_type i, const size_type j) const;
271 
276  unsigned int
277  row_length(const size_type row) const;
278 
290  size_type
291  n_nonzero_elements() const;
292 
298  void
299  print(std::ostream &out) const;
300 
308  void
309  print_gnuplot(std::ostream &out) const;
310 
320  int,
321  int,
322  int,
323  int,
324  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
325  << ',' << arg4 << "] have differing row numbers.");
330  int,
331  int,
332  int,
333  int,
334  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
335  << ',' << arg4 << "] have differing column numbers.");
337 
338 protected:
343 
348 
352  Table<2,
353  SmartPointer<SparsityPatternType,
356 
362 
368 
369 private:
374  std::vector<size_type> counter_within_block;
375 
380  std::vector<std::vector<size_type>> block_column_indices;
381 
382  // Make the block sparse matrix a friend, so that it can use our
383  // #row_indices and #column_indices objects.
384  template <typename number>
385  friend class BlockSparseMatrix;
386 };
387 
388 
389 
401 class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
402 {
403 public:
409  BlockSparsityPattern() = default;
410 
416  BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
417 
421  void
422  reinit(const size_type n_block_rows, const size_type n_block_columns);
423 
436  void
438  const BlockIndices & col_indices,
439  const std::vector<std::vector<unsigned int>> &row_lengths);
440 
441 
446  bool
447  is_compressed() const;
448 
453  std::size_t
454  memory_consumption() const;
455 
461  void
462  copy_from(const BlockDynamicSparsityPattern &dsp);
463 };
464 
465 
466 
522  : public BlockSparsityPatternBase<DynamicSparsityPattern>
523 {
524 public:
530  BlockDynamicSparsityPattern() = default;
531 
538  const size_type n_columns);
539 
546  BlockDynamicSparsityPattern(const std::vector<size_type> &row_block_sizes,
547  const std::vector<size_type> &col_block_sizes);
548 
557  BlockDynamicSparsityPattern(const std::vector<IndexSet> &partitioning);
558 
564  const BlockIndices &col_indices);
565 
566 
575  void
576  reinit(const std::vector<size_type> &row_block_sizes,
577  const std::vector<size_type> &col_block_sizes);
578 
583  void
584  reinit(const std::vector<IndexSet> &partitioning);
585 
591  void
592  reinit(const BlockIndices &row_indices, const BlockIndices &col_indices);
593 
598  size_type
599  column_number(const size_type row, const unsigned int index) const;
600 
605 };
606 
610 #ifdef DEAL_II_WITH_TRILINOS
611 
612 
613 namespace TrilinosWrappers
614 {
639  : public ::BlockSparsityPatternBase<SparsityPattern>
640  {
641  public:
647  BlockSparsityPattern() = default;
648 
654  BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
655 
662  BlockSparsityPattern(const std::vector<size_type> &row_block_sizes,
663  const std::vector<size_type> &col_block_sizes);
664 
672  BlockSparsityPattern(const std::vector<IndexSet> &parallel_partitioning,
673  const MPI_Comm &communicator = MPI_COMM_WORLD);
674 
687  const std::vector<IndexSet> &row_parallel_partitioning,
688  const std::vector<IndexSet> &column_parallel_partitioning,
689  const std::vector<IndexSet> &writeable_rows,
690  const MPI_Comm & communicator = MPI_COMM_WORLD);
691 
701  void
702  reinit(const std::vector<size_type> &row_block_sizes,
703  const std::vector<size_type> &col_block_sizes);
704 
709  void
710  reinit(const std::vector<IndexSet> &parallel_partitioning,
711  const MPI_Comm & communicator = MPI_COMM_WORLD);
712 
718  void
719  reinit(const std::vector<IndexSet> &row_parallel_partitioning,
720  const std::vector<IndexSet> &column_parallel_partitioning,
721  const MPI_Comm & communicator = MPI_COMM_WORLD);
722 
731  void
732  reinit(const std::vector<IndexSet> &row_parallel_partitioning,
733  const std::vector<IndexSet> &column_parallel_partitioning,
734  const std::vector<IndexSet> &writeable_rows,
735  const MPI_Comm & communicator = MPI_COMM_WORLD);
736 
741  };
742 
745 } /* namespace TrilinosWrappers */
746 
747 #endif
748 
749 /*--------------------- Template functions ----------------------------------*/
750 
751 
752 
753 template <typename SparsityPatternType>
754 inline SparsityPatternType &
756  const size_type column)
757 {
758  AssertIndexRange(row, rows);
759  AssertIndexRange(column, columns);
760  return *sub_objects[row][column];
761 }
762 
763 
764 
765 template <typename SparsityPatternType>
766 inline const SparsityPatternType &
768  const size_type row,
769  const size_type column) const
770 {
771  AssertIndexRange(row, rows);
772  AssertIndexRange(column, columns);
773  return *sub_objects[row][column];
774 }
775 
776 
777 
778 template <typename SparsityPatternType>
779 inline const BlockIndices &
781 {
782  return row_indices;
783 }
784 
785 
786 
787 template <typename SparsityPatternType>
788 inline const BlockIndices &
790 {
791  return column_indices;
792 }
793 
794 
795 
796 template <typename SparsityPatternType>
797 inline void
799  const size_type j)
800 {
801  // if you get an error here, are
802  // you sure you called
803  // <tt>collect_sizes()</tt> before?
804  const std::pair<size_type, size_type> row_index =
806  col_index =
808  sub_objects[row_index.first][col_index.first]->add(row_index.second,
809  col_index.second);
810 }
811 
812 
813 
814 template <typename SparsityPatternType>
815 template <typename ForwardIterator>
816 void
818  const size_type row,
819  ForwardIterator begin,
820  ForwardIterator end,
821  const bool indices_are_sorted)
822 {
823  // Resize scratch arrays
824  if (block_column_indices.size() < this->n_block_cols())
825  {
826  block_column_indices.resize(this->n_block_cols());
827  counter_within_block.resize(this->n_block_cols());
828  }
829 
830  const size_type n_cols = static_cast<size_type>(end - begin);
831 
832  // Resize sub-arrays to n_cols. This
833  // is a bit wasteful, but we resize
834  // only a few times (then the maximum
835  // row length won't increase that
836  // much any more). At least we know
837  // that all arrays are going to be of
838  // the same size, so we can check
839  // whether the size of one is large
840  // enough before actually going
841  // through all of them.
842  if (block_column_indices[0].size() < n_cols)
843  for (size_type i = 0; i < this->n_block_cols(); ++i)
844  block_column_indices[i].resize(n_cols);
845 
846  // Reset the number of added elements
847  // in each block to zero.
848  for (size_type i = 0; i < this->n_block_cols(); ++i)
849  counter_within_block[i] = 0;
850 
851  // Go through the column indices to
852  // find out which portions of the
853  // values should be set in which
854  // block of the matrix. We need to
855  // touch all the data, since we can't
856  // be sure that the data of one block
857  // is stored contiguously (in fact,
858  // indices will be intermixed when it
859  // comes from an element matrix).
860  for (ForwardIterator it = begin; it != end; ++it)
861  {
862  const size_type col = *it;
863 
864  const std::pair<size_type, size_type> col_index =
865  this->column_indices.global_to_local(col);
866 
867  const size_type local_index = counter_within_block[col_index.first]++;
868 
869  block_column_indices[col_index.first][local_index] = col_index.second;
870  }
871 
872 #ifdef DEBUG
873  // If in debug mode, do a check whether
874  // the right length has been obtained.
875  size_type length = 0;
876  for (size_type i = 0; i < this->n_block_cols(); ++i)
877  length += counter_within_block[i];
878  Assert(length == n_cols, ExcInternalError());
879 #endif
880 
881  // Now we found out about where the
882  // individual columns should start and
883  // where we should start reading out
884  // data. Now let's write the data into
885  // the individual blocks!
886  const std::pair<size_type, size_type> row_index =
887  this->row_indices.global_to_local(row);
888  for (size_type block_col = 0; block_col < n_block_cols(); ++block_col)
889  {
890  if (counter_within_block[block_col] == 0)
891  continue;
892  sub_objects[row_index.first][block_col]->add_entries(
893  row_index.second,
894  block_column_indices[block_col].begin(),
895  block_column_indices[block_col].begin() +
896  counter_within_block[block_col],
897  indices_are_sorted);
898  }
899 }
900 
901 
902 
903 template <typename SparsityPatternType>
904 inline bool
906  const size_type j) const
907 {
908  // if you get an error here, are
909  // you sure you called
910  // <tt>collect_sizes()</tt> before?
911  const std::pair<size_type, size_type> row_index =
913  col_index =
915  return sub_objects[row_index.first][col_index.first]->exists(
916  row_index.second, col_index.second);
917 }
918 
919 
920 
921 template <typename SparsityPatternType>
922 inline unsigned int
924  const size_type row) const
925 {
926  const std::pair<size_type, size_type> row_index =
928 
929  unsigned int c = 0;
930 
931  for (size_type b = 0; b < rows; ++b)
932  c += sub_objects[row_index.first][b]->row_length(row_index.second);
933 
934  return c;
935 }
936 
937 
938 
939 template <typename SparsityPatternType>
942 {
943  return columns;
944 }
945 
946 
947 
948 template <typename SparsityPatternType>
951 {
952  return rows;
953 }
954 
955 
958  const unsigned int index) const
959 {
960  // .first= ith block, .second = jth row in that block
961  const std::pair<size_type, size_type> row_index =
963 
964  AssertIndexRange(index, row_length(row));
965 
966  size_type c = 0;
967  size_type block_columns = 0; // sum of n_cols for all blocks to the left
968  for (unsigned int b = 0; b < columns; ++b)
969  {
970  unsigned int rowlen =
971  sub_objects[row_index.first][b]->row_length(row_index.second);
972  if (index < c + rowlen)
973  return block_columns +
974  sub_objects[row_index.first][b]->column_number(row_index.second,
975  index - c);
976  c += rowlen;
977  block_columns += sub_objects[row_index.first][b]->n_cols();
978  }
979 
980  Assert(false, ExcInternalError());
981  return 0;
982 }
983 
984 
985 inline void
987  const size_type n_block_columns)
988 {
990  n_block_columns);
991 }
992 
993 
995 
996 #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:1638
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:1408
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:361
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)
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:360
VectorType::value_type * begin(VectorType &V)
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Definition: exceptions.h:587
unsigned int row_length(const size_type row) const
void print_gnuplot(std::ostream &out) const
Definition: table.h:695
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()