Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
block_sparsity_pattern.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 2018 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 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/base/subscriptor.h>
25 #include <deal.II/base/table.h>
26 
27 #include <deal.II/lac/block_indices.h>
28 #include <deal.II/lac/dynamic_sparsity_pattern.h>
29 #include <deal.II/lac/sparsity_pattern.h>
30 #include <deal.II/lac/trilinos_sparsity_pattern.h>
31 
32 DEAL_II_NAMESPACE_OPEN
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 
510  : public BlockSparsityPatternBase<DynamicSparsityPattern>
511 {
512 public:
518  BlockDynamicSparsityPattern() = default;
519 
526  const size_type n_columns);
527 
534  BlockDynamicSparsityPattern(const std::vector<size_type> &row_block_sizes,
535  const std::vector<size_type> &col_block_sizes);
536 
545  BlockDynamicSparsityPattern(const std::vector<IndexSet> &partitioning);
546 
552  const BlockIndices &col_indices);
553 
554 
563  void
564  reinit(const std::vector<size_type> &row_block_sizes,
565  const std::vector<size_type> &col_block_sizes);
566 
571  void
572  reinit(const std::vector<IndexSet> &partitioning);
573 
579  void
580  reinit(const BlockIndices &row_indices, const BlockIndices &col_indices);
581 
586  size_type
587  column_number(const size_type row, const unsigned int index) const;
588 
593 };
594 
598 #ifdef DEAL_II_WITH_TRILINOS
599 
600 
601 namespace TrilinosWrappers
602 {
627  : public ::BlockSparsityPatternBase<SparsityPattern>
628  {
629  public:
635  BlockSparsityPattern() = default;
636 
642  BlockSparsityPattern(const size_type n_rows, const size_type n_columns);
643 
650  BlockSparsityPattern(const std::vector<size_type> &row_block_sizes,
651  const std::vector<size_type> &col_block_sizes);
652 
663  DEAL_II_DEPRECATED
664  BlockSparsityPattern(const std::vector<Epetra_Map> &parallel_partitioning);
665 
673  BlockSparsityPattern(const std::vector<IndexSet> &parallel_partitioning,
674  const MPI_Comm &communicator = MPI_COMM_WORLD);
675 
687  BlockSparsityPattern(
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 
713  DEAL_II_DEPRECATED
714  void
715  reinit(const std::vector<Epetra_Map> &parallel_partitioning);
716 
721  void
722  reinit(const std::vector<IndexSet> &parallel_partitioning,
723  const MPI_Comm & communicator = MPI_COMM_WORLD);
724 
730  void
731  reinit(const std::vector<IndexSet> &row_parallel_partitioning,
732  const std::vector<IndexSet> &column_parallel_partitioning,
733  const MPI_Comm & communicator = MPI_COMM_WORLD);
734 
743  void
744  reinit(const std::vector<IndexSet> &row_parallel_partitioning,
745  const std::vector<IndexSet> &column_parallel_partitioning,
746  const std::vector<IndexSet> &writeable_rows,
747  const MPI_Comm & communicator = MPI_COMM_WORLD);
748 
753  };
754 
757 } /* namespace TrilinosWrappers */
758 
759 #endif
760 
761 /*--------------------- Template functions ----------------------------------*/
762 
763 
764 
765 template <typename SparsityPatternType>
766 inline SparsityPatternType &
768  const size_type column)
769 {
770  Assert(row < rows, ExcIndexRange(row, 0, rows));
771  Assert(column < columns, ExcIndexRange(column, 0, columns));
772  return *sub_objects[row][column];
773 }
774 
775 
776 
777 template <typename SparsityPatternType>
778 inline const SparsityPatternType &
780  const size_type row,
781  const size_type column) const
782 {
783  Assert(row < rows, ExcIndexRange(row, 0, rows));
784  Assert(column < columns, ExcIndexRange(column, 0, columns));
785  return *sub_objects[row][column];
786 }
787 
788 
789 
790 template <typename SparsityPatternType>
791 inline const BlockIndices &
793 {
794  return row_indices;
795 }
796 
797 
798 
799 template <typename SparsityPatternType>
800 inline const BlockIndices &
802 {
803  return column_indices;
804 }
805 
806 
807 
808 template <typename SparsityPatternType>
809 inline void
811  const size_type j)
812 {
813  // if you get an error here, are
814  // you sure you called
815  // <tt>collect_sizes()</tt> before?
816  const std::pair<size_type, size_type> row_index =
818  col_index =
820  sub_objects[row_index.first][col_index.first]->add(row_index.second,
821  col_index.second);
822 }
823 
824 
825 
826 template <typename SparsityPatternType>
827 template <typename ForwardIterator>
828 void
830  const size_type row,
831  ForwardIterator begin,
832  ForwardIterator end,
833  const bool indices_are_sorted)
834 {
835  // Resize scratch arrays
836  if (block_column_indices.size() < this->n_block_cols())
837  {
838  block_column_indices.resize(this->n_block_cols());
839  counter_within_block.resize(this->n_block_cols());
840  }
841 
842  const size_type n_cols = static_cast<size_type>(end - begin);
843 
844  // Resize sub-arrays to n_cols. This
845  // is a bit wasteful, but we resize
846  // only a few times (then the maximum
847  // row length won't increase that
848  // much any more). At least we know
849  // that all arrays are going to be of
850  // the same size, so we can check
851  // whether the size of one is large
852  // enough before actually going
853  // through all of them.
854  if (block_column_indices[0].size() < n_cols)
855  for (size_type i = 0; i < this->n_block_cols(); ++i)
856  block_column_indices[i].resize(n_cols);
857 
858  // Reset the number of added elements
859  // in each block to zero.
860  for (size_type i = 0; i < this->n_block_cols(); ++i)
861  counter_within_block[i] = 0;
862 
863  // Go through the column indices to
864  // find out which portions of the
865  // values should be set in which
866  // block of the matrix. We need to
867  // touch all the data, since we can't
868  // be sure that the data of one block
869  // is stored contiguously (in fact,
870  // indices will be intermixed when it
871  // comes from an element matrix).
872  for (ForwardIterator it = begin; it != end; ++it)
873  {
874  const size_type col = *it;
875 
876  const std::pair<size_type, size_type> col_index =
877  this->column_indices.global_to_local(col);
878 
879  const size_type local_index = counter_within_block[col_index.first]++;
880 
881  block_column_indices[col_index.first][local_index] = col_index.second;
882  }
883 
884 #ifdef DEBUG
885  // If in debug mode, do a check whether
886  // the right length has been obtained.
887  size_type length = 0;
888  for (size_type i = 0; i < this->n_block_cols(); ++i)
889  length += counter_within_block[i];
890  Assert(length == n_cols, ExcInternalError());
891 #endif
892 
893  // Now we found out about where the
894  // individual columns should start and
895  // where we should start reading out
896  // data. Now let's write the data into
897  // the individual blocks!
898  const std::pair<size_type, size_type> row_index =
899  this->row_indices.global_to_local(row);
900  for (size_type block_col = 0; block_col < n_block_cols(); ++block_col)
901  {
902  if (counter_within_block[block_col] == 0)
903  continue;
904  sub_objects[row_index.first][block_col]->add_entries(
905  row_index.second,
906  block_column_indices[block_col].begin(),
907  block_column_indices[block_col].begin() +
908  counter_within_block[block_col],
909  indices_are_sorted);
910  }
911 }
912 
913 
914 
915 template <typename SparsityPatternType>
916 inline bool
918  const size_type j) const
919 {
920  // if you get an error here, are
921  // you sure you called
922  // <tt>collect_sizes()</tt> before?
923  const std::pair<size_type, size_type> row_index =
925  col_index =
927  return sub_objects[row_index.first][col_index.first]->exists(
928  row_index.second, col_index.second);
929 }
930 
931 
932 
933 template <typename SparsityPatternType>
934 inline unsigned int
936  const size_type row) const
937 {
938  const std::pair<size_type, size_type> row_index =
940 
941  unsigned int c = 0;
942 
943  for (size_type b = 0; b < rows; ++b)
944  c += sub_objects[row_index.first][b]->row_length(row_index.second);
945 
946  return c;
947 }
948 
949 
950 
951 template <typename SparsityPatternType>
954 {
955  return columns;
956 }
957 
958 
959 
960 template <typename SparsityPatternType>
963 {
964  return rows;
965 }
966 
967 
970  const unsigned int index) const
971 {
972  // .first= ith block, .second = jth row in that block
973  const std::pair<size_type, size_type> row_index =
975 
976  Assert(index < row_length(row), ExcIndexRange(index, 0, row_length(row)));
977 
978  size_type c = 0;
979  size_type block_columns = 0; // sum of n_cols for all blocks to the left
980  for (unsigned int b = 0; b < columns; ++b)
981  {
982  unsigned int rowlen =
983  sub_objects[row_index.first][b]->row_length(row_index.second);
984  if (index < c + rowlen)
985  return block_columns +
986  sub_objects[row_index.first][b]->column_number(row_index.second,
987  index - c);
988  c += rowlen;
989  block_columns += sub_objects[row_index.first][b]->n_cols();
990  }
991 
992  Assert(false, ExcInternalError());
993  return 0;
994 }
995 
996 
997 inline void
999  const size_type n_block_columns)
1000 {
1002  n_block_columns);
1003 }
1004 
1005 
1006 DEAL_II_NAMESPACE_CLOSE
1007 
1008 #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
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
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:1407
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
const BlockIndices & get_row_indices() const
BlockSparsityPatternBase & operator=(const BlockSparsityPatternBase &)
static ::ExceptionBase & ExcIncompatibleColNumbers(int arg1, int arg2, int arg3, int arg4)
unsigned int global_dof_index
Definition: types.h:89
Table< 2, SmartPointer< SparsityPatternType, BlockSparsityPatternBase< SparsityPatternType > > > sub_objects
#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:699
static const size_type invalid_entry
std::pair< unsigned int, size_type > global_to_local(const size_type i) const
const BlockIndices & get_column_indices() const
static ::ExceptionBase & ExcInternalError()