Reference documentation for deal.II version Git e7bb9ce7b3 2020-09-18 12:07:32 -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 
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 
319  int,
320  int,
321  int,
322  int,
323  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
324  << ',' << arg4 << "] have differing row numbers.");
329  int,
330  int,
331  int,
332  int,
333  << "The blocks [" << arg1 << ',' << arg2 << "] and [" << arg3
334  << ',' << arg4 << "] have differing column numbers.");
336 
337 protected:
342 
347 
351  Table<2,
352  SmartPointer<SparsityPatternType,
355 
361 
367 
368 private:
373  std::vector<size_type> counter_within_block;
374 
379  std::vector<std::vector<size_type>> block_column_indices;
380 
381  // Make the block sparse matrix a friend, so that it can use our
382  // #row_indices and #column_indices objects.
383  template <typename number>
384  friend class BlockSparseMatrix;
385 };
386 
387 
388 
398 class BlockSparsityPattern : public BlockSparsityPatternBase<SparsityPattern>
399 {
400 public:
406  BlockSparsityPattern() = default;
407 
413  BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
414 
418  void
419  reinit(const size_type n_block_rows, const size_type n_block_columns);
420 
433  void
435  const BlockIndices & col_indices,
436  const std::vector<std::vector<unsigned int>> &row_lengths);
437 
438 
443  bool
444  is_compressed() const;
445 
450  std::size_t
451  memory_consumption() const;
452 
458  void
459  copy_from(const BlockDynamicSparsityPattern &dsp);
460 };
461 
462 
463 
517  : public BlockSparsityPatternBase<DynamicSparsityPattern>
518 {
519 public:
525  BlockDynamicSparsityPattern() = default;
526 
533  const size_type n_columns);
534 
541  BlockDynamicSparsityPattern(const std::vector<size_type> &row_block_sizes,
542  const std::vector<size_type> &col_block_sizes);
543 
552  BlockDynamicSparsityPattern(const std::vector<IndexSet> &partitioning);
553 
559  const BlockIndices &col_indices);
560 
561 
570  void
571  reinit(const std::vector<size_type> &row_block_sizes,
572  const std::vector<size_type> &col_block_sizes);
573 
578  void
579  reinit(const std::vector<IndexSet> &partitioning);
580 
586  void
587  reinit(const BlockIndices &row_indices, const BlockIndices &col_indices);
588 
593  size_type
594  column_number(const size_type row, const unsigned int index) const;
595 
600 };
601 
605 #ifdef DEAL_II_WITH_TRILINOS
606 
607 
608 namespace TrilinosWrappers
609 {
632  : public ::BlockSparsityPatternBase<SparsityPattern>
633  {
634  public:
640  BlockSparsityPattern() = default;
641 
647  BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
648 
655  BlockSparsityPattern(const std::vector<size_type> &row_block_sizes,
656  const std::vector<size_type> &col_block_sizes);
657 
665  BlockSparsityPattern(const std::vector<IndexSet> &parallel_partitioning,
666  const MPI_Comm &communicator = MPI_COMM_WORLD);
667 
680  const std::vector<IndexSet> &row_parallel_partitioning,
681  const std::vector<IndexSet> &column_parallel_partitioning,
682  const std::vector<IndexSet> &writeable_rows,
683  const MPI_Comm & communicator = MPI_COMM_WORLD);
684 
694  void
695  reinit(const std::vector<size_type> &row_block_sizes,
696  const std::vector<size_type> &col_block_sizes);
697 
702  void
703  reinit(const std::vector<IndexSet> &parallel_partitioning,
704  const MPI_Comm & communicator = MPI_COMM_WORLD);
705 
711  void
712  reinit(const std::vector<IndexSet> &row_parallel_partitioning,
713  const std::vector<IndexSet> &column_parallel_partitioning,
714  const MPI_Comm & communicator = MPI_COMM_WORLD);
715 
724  void
725  reinit(const std::vector<IndexSet> &row_parallel_partitioning,
726  const std::vector<IndexSet> &column_parallel_partitioning,
727  const std::vector<IndexSet> &writeable_rows,
728  const MPI_Comm & communicator = MPI_COMM_WORLD);
729 
734  };
735 
738 } /* namespace TrilinosWrappers */
739 
740 #endif
741 
742 /*--------------------- Template functions ----------------------------------*/
743 
744 
745 
746 template <typename SparsityPatternType>
747 inline SparsityPatternType &
749  const size_type column)
750 {
751  AssertIndexRange(row, rows);
752  AssertIndexRange(column, columns);
753  return *sub_objects[row][column];
754 }
755 
756 
757 
758 template <typename SparsityPatternType>
759 inline const SparsityPatternType &
761  const size_type row,
762  const size_type column) const
763 {
764  AssertIndexRange(row, rows);
765  AssertIndexRange(column, columns);
766  return *sub_objects[row][column];
767 }
768 
769 
770 
771 template <typename SparsityPatternType>
772 inline const BlockIndices &
774 {
775  return row_indices;
776 }
777 
778 
779 
780 template <typename SparsityPatternType>
781 inline const BlockIndices &
783 {
784  return column_indices;
785 }
786 
787 
788 
789 template <typename SparsityPatternType>
790 inline void
792  const size_type j)
793 {
794  // if you get an error here, are
795  // you sure you called
796  // <tt>collect_sizes()</tt> before?
797  const std::pair<size_type, size_type> row_index =
799  col_index =
801  sub_objects[row_index.first][col_index.first]->add(row_index.second,
802  col_index.second);
803 }
804 
805 
806 
807 template <typename SparsityPatternType>
808 template <typename ForwardIterator>
809 void
811  const size_type row,
812  ForwardIterator begin,
813  ForwardIterator end,
814  const bool indices_are_sorted)
815 {
816  // Resize scratch arrays
817  if (block_column_indices.size() < this->n_block_cols())
818  {
819  block_column_indices.resize(this->n_block_cols());
820  counter_within_block.resize(this->n_block_cols());
821  }
822 
823  const size_type n_cols = static_cast<size_type>(end - begin);
824 
825  // Resize sub-arrays to n_cols. This
826  // is a bit wasteful, but we resize
827  // only a few times (then the maximum
828  // row length won't increase that
829  // much any more). At least we know
830  // that all arrays are going to be of
831  // the same size, so we can check
832  // whether the size of one is large
833  // enough before actually going
834  // through all of them.
835  if (block_column_indices[0].size() < n_cols)
836  for (size_type i = 0; i < this->n_block_cols(); ++i)
837  block_column_indices[i].resize(n_cols);
838 
839  // Reset the number of added elements
840  // in each block to zero.
841  for (size_type i = 0; i < this->n_block_cols(); ++i)
842  counter_within_block[i] = 0;
843 
844  // Go through the column indices to
845  // find out which portions of the
846  // values should be set in which
847  // block of the matrix. We need to
848  // touch all the data, since we can't
849  // be sure that the data of one block
850  // is stored contiguously (in fact,
851  // indices will be intermixed when it
852  // comes from an element matrix).
853  for (ForwardIterator it = begin; it != end; ++it)
854  {
855  const size_type col = *it;
856 
857  const std::pair<size_type, size_type> col_index =
858  this->column_indices.global_to_local(col);
859 
860  const size_type local_index = counter_within_block[col_index.first]++;
861 
862  block_column_indices[col_index.first][local_index] = col_index.second;
863  }
864 
865 #ifdef DEBUG
866  // If in debug mode, do a check whether
867  // the right length has been obtained.
868  size_type length = 0;
869  for (size_type i = 0; i < this->n_block_cols(); ++i)
870  length += counter_within_block[i];
871  Assert(length == n_cols, ExcInternalError());
872 #endif
873 
874  // Now we found out about where the
875  // individual columns should start and
876  // where we should start reading out
877  // data. Now let's write the data into
878  // the individual blocks!
879  const std::pair<size_type, size_type> row_index =
880  this->row_indices.global_to_local(row);
881  for (size_type block_col = 0; block_col < n_block_cols(); ++block_col)
882  {
883  if (counter_within_block[block_col] == 0)
884  continue;
885  sub_objects[row_index.first][block_col]->add_entries(
886  row_index.second,
887  block_column_indices[block_col].begin(),
888  block_column_indices[block_col].begin() +
889  counter_within_block[block_col],
890  indices_are_sorted);
891  }
892 }
893 
894 
895 
896 template <typename SparsityPatternType>
897 inline bool
899  const size_type j) const
900 {
901  // if you get an error here, are
902  // you sure you called
903  // <tt>collect_sizes()</tt> before?
904  const std::pair<size_type, size_type> row_index =
906  col_index =
908  return sub_objects[row_index.first][col_index.first]->exists(
909  row_index.second, col_index.second);
910 }
911 
912 
913 
914 template <typename SparsityPatternType>
915 inline unsigned int
917  const size_type row) const
918 {
919  const std::pair<size_type, size_type> row_index =
921 
922  unsigned int c = 0;
923 
924  for (size_type b = 0; b < rows; ++b)
925  c += sub_objects[row_index.first][b]->row_length(row_index.second);
926 
927  return c;
928 }
929 
930 
931 
932 template <typename SparsityPatternType>
935 {
936  return columns;
937 }
938 
939 
940 
941 template <typename SparsityPatternType>
944 {
945  return rows;
946 }
947 
948 
951  const unsigned int index) const
952 {
953  // .first= ith block, .second = jth row in that block
954  const std::pair<size_type, size_type> row_index =
956 
957  AssertIndexRange(index, row_length(row));
958 
959  size_type c = 0;
960  size_type block_columns = 0; // sum of n_cols for all blocks to the left
961  for (unsigned int b = 0; b < columns; ++b)
962  {
963  unsigned int rowlen =
964  sub_objects[row_index.first][b]->row_length(row_index.second);
965  if (index < c + rowlen)
966  return block_columns +
967  sub_objects[row_index.first][b]->column_number(row_index.second,
968  index - c);
969  c += rowlen;
970  block_columns += sub_objects[row_index.first][b]->n_cols();
971  }
972 
973  Assert(false, ExcInternalError());
974  return 0;
975 }
976 
977 
978 inline void
980  const size_type n_block_columns)
981 {
983  n_block_columns);
984 }
985 
986 
988 
989 #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:1636
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:1411
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:363
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:362
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()