Reference documentation for deal.II version 8.4.1
dynamic_sparsity_pattern.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2011 - 2015 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 at
12 // the top level of the deal.II distribution.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii__dynamic_sparsity_pattern_h
17 #define dealii__dynamic_sparsity_pattern_h
18 
19 
20 #include <deal.II/base/config.h>
21 #include <deal.II/base/subscriptor.h>
22 #include <deal.II/base/utilities.h>
23 #include <deal.II/lac/exceptions.h>
24 #include <deal.II/base/index_set.h>
25 
26 #include <vector>
27 #include <algorithm>
28 #include <iostream>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
32 template <typename number> class SparseMatrix;
34 
35 
45 {
46  // forward declaration
47  class Iterator;
48 
53 
65  class Accessor
66  {
67  public:
72  const size_type row,
73  const unsigned int index_within_row);
74 
78  Accessor (const DynamicSparsityPattern *sparsity_pattern);
79 
83  size_type row () const;
84 
88  size_type index () const;
89 
93  size_type column () const;
94 
98  bool operator == (const Accessor &) const;
99 
107  bool operator < (const Accessor &) const;
108 
109  protected:
114 
118  size_type current_row;
119 
124  std::vector<size_type>::const_iterator current_entry;
125 
132  std::vector<size_type>::const_iterator end_of_row;
133 
137  void advance ();
138 
142  friend class Iterator;
143  };
144 
145 
146 
171  class Iterator
172  {
173  public:
179  Iterator (const DynamicSparsityPattern *sp,
180  const size_type row,
181  const unsigned int index_within_row);
182 
187  Iterator (const DynamicSparsityPattern *sp);
188 
193 
197  Iterator operator++ (int);
198 
202  const Accessor &operator* () const;
203 
207  const Accessor *operator-> () const;
208 
212  bool operator == (const Iterator &) const;
213 
217  bool operator != (const Iterator &) const;
218 
226  bool operator < (const Iterator &) const;
227 
234  int operator - (const Iterator &p) const;
235 
236  private:
241  };
242 }
243 
244 
291 {
292 public:
297 
305  typedef
308 
313  typedef
316 
323 
334 
341  DynamicSparsityPattern (const size_type m,
342  const size_type n,
343  const IndexSet &rowset = IndexSet());
344 
348  DynamicSparsityPattern (const IndexSet &indexset);
349 
353  DynamicSparsityPattern (const size_type n);
354 
361 
368  void reinit (const size_type m,
369  const size_type n,
370  const IndexSet &rowset = IndexSet());
371 
377  void compress ();
378 
383  bool empty () const;
384 
389  size_type max_entries_per_row () const;
390 
394  void add (const size_type i,
395  const size_type j);
396 
401  template <typename ForwardIterator>
402  void add_entries (const size_type row,
403  ForwardIterator begin,
404  ForwardIterator end,
405  const bool indices_are_unique_and_sorted = false);
406 
410  bool exists (const size_type i,
411  const size_type j) const;
412 
420  void symmetrize ();
421 
427  void print (std::ostream &out) const;
428 
442  void print_gnuplot (std::ostream &out) const;
443 
447  size_type n_rows () const;
448 
453  size_type n_cols () const;
454 
459  size_type row_length (const size_type row) const;
460 
465  size_type column_number (const size_type row,
466  const size_type index) const;
467 
471 // @{
472 
486  iterator begin () const;
487 
491  iterator end () const;
492 
509  iterator begin (const size_type r) const;
510 
519  iterator end (const size_type r) const;
520 
521 // @}
522 
528  size_type bandwidth () const;
529 
534  size_type n_nonzero_elements () const;
535 
541  const IndexSet &row_index_set () const;
542 
553  static
555 
560  size_type memory_consumption () const;
561 
562 private:
566  size_type rows;
567 
571  size_type cols;
572 
578 
579 
586  struct Line
587  {
588  public:
593  std::vector<size_type> entries;
594 
598  Line ();
599 
603  void add (const size_type col_num);
604 
608  template <typename ForwardIterator>
609  void add_entries (ForwardIterator begin,
610  ForwardIterator end,
611  const bool indices_are_sorted);
612 
616  size_type memory_consumption () const;
617  };
618 
619 
623  std::vector<Line> lines;
624 
625  // make the accessor class a friend
627 };
628 
630 /*---------------------- Inline functions -----------------------------------*/
631 
632 
634 {
635  inline
636  Accessor::
637  Accessor (const DynamicSparsityPattern *sparsity_pattern,
638  const size_type row,
639  const unsigned int index_within_row)
640  :
641  sparsity_pattern(sparsity_pattern),
642  current_row (row),
643  current_entry(((sparsity_pattern->rowset.size()==0)
644  ?
645  sparsity_pattern->lines[current_row].entries.begin()
646  :
647  sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.begin())
648  +
649  index_within_row),
650  end_of_row((sparsity_pattern->rowset.size()==0)
651  ?
652  sparsity_pattern->lines[current_row].entries.end()
653  :
654  sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.end())
655  {
656  AssertIndexRange(current_row, sparsity_pattern->n_rows());
657  Assert ((sparsity_pattern->rowset.size()==0)
658  ||
659  sparsity_pattern->rowset.is_element(current_row),
660  ExcMessage ("You can't create an iterator into a "
661  "DynamicSparsityPattern's row that is not "
662  "actually stored by that sparsity pattern "
663  "based on the IndexSet argument to it."));
664  AssertIndexRange(index_within_row,
665  ((sparsity_pattern->rowset.size()==0)
666  ?
667  sparsity_pattern->lines[current_row].entries.size()
668  :
669  sparsity_pattern->lines[sparsity_pattern->rowset.index_within_set(current_row)].entries.size()));
670  }
671 
672 
673  inline
674  Accessor::
675  Accessor (const DynamicSparsityPattern *sparsity_pattern)
676  :
677  sparsity_pattern(sparsity_pattern),
678  current_row(numbers::invalid_size_type),
679  current_entry(),
680  end_of_row()
681  {}
682 
683 
684 
685  inline
686  size_type
688  {
689  Assert (current_row < sparsity_pattern->n_rows(),
690  ExcInternalError());
691 
692  return current_row;
693  }
694 
695 
696  inline
697  size_type
699  {
700  Assert (current_row < sparsity_pattern->n_rows(),
701  ExcInternalError());
702 
703  return *current_entry;
704  }
705 
706 
707  inline
708  size_type
710  {
711  Assert (current_row < sparsity_pattern->n_rows(),
712  ExcInternalError());
713 
714  return (current_entry -
715  ((sparsity_pattern->rowset.size()==0)
716  ?
717  sparsity_pattern->lines[current_row].entries.begin()
718  :
720  }
721 
722 
723 
724 
725  inline
726  bool
727  Accessor::operator == (const Accessor &other) const
728  {
729  // compare the sparsity pattern the iterator points into, the
730  // current row, and the location within this row. ignore the
731  // latter if the row is past-the-end because in that case the
732  // current_entry field may not point to a deterministic location
733  return (sparsity_pattern == other.sparsity_pattern &&
734  current_row == other.current_row &&
736  || (current_entry == other.current_entry)));
737  }
738 
739 
740 
741  inline
742  bool
743  Accessor::operator < (const Accessor &other) const
744  {
746  ExcInternalError());
747 
748  // if *this is past-the-end, then it is less than no one
750  return (false);
751  // now *this should be an valid value
752  Assert (current_row < sparsity_pattern->n_rows(),
753  ExcInternalError());
754 
755  // if other is past-the-end
757  return (true);
758  // now other should be an valid value
760  ExcInternalError());
761 
762  // both iterators are not one-past-the-end
763  return ((current_row < other.current_row) ||
764  ((current_row == other.current_row) &&
765  (current_entry < other.current_entry)));
766  }
767 
768 
769  inline
770  void
772  {
773  Assert (current_row < sparsity_pattern->n_rows(),
774  ExcInternalError());
775 
776  // move to the next element in this row
777  ++current_entry;
778 
779  // if this moves us beyond the end of the row, go to the next row
780  // if possible, or set the iterator to an invalid state if not.
781  //
782  // going to the next row is a bit complicated because we may have
783  // to skip over empty rows, and because we also have to avoid rows
784  // that aren't listed in a possibly passed IndexSet argument of
785  // the sparsity pattern. consequently, rather than trying to
786  // duplicate code here, just call the begin() function of the
787  // sparsity pattern itself
788  if (current_entry == end_of_row)
789  {
791  *this = *sparsity_pattern->begin(current_row+1);
792  else
793  *this = Accessor(sparsity_pattern); // invalid object
794  }
795  }
796 
797 
798 
799  inline
800  Iterator::Iterator (const DynamicSparsityPattern *sparsity_pattern,
801  const size_type row,
802  const unsigned int index_within_row)
803  :
804  accessor(sparsity_pattern, row, index_within_row)
805  {}
806 
807 
808 
809  inline
810  Iterator::Iterator (const DynamicSparsityPattern *sparsity_pattern)
811  :
812  accessor(sparsity_pattern)
813  {}
814 
815 
816 
817  inline
818  Iterator &
820  {
821  accessor.advance ();
822  return *this;
823  }
824 
825 
826 
827  inline
828  Iterator
830  {
831  const Iterator iter = *this;
832  accessor.advance ();
833  return iter;
834  }
835 
836 
837 
838  inline
839  const Accessor &
841  {
842  return accessor;
843  }
844 
845 
846 
847  inline
848  const Accessor *
850  {
851  return &accessor;
852  }
853 
854 
855  inline
856  bool
857  Iterator::operator == (const Iterator &other) const
858  {
859  return (accessor == other.accessor);
860  }
861 
862 
863 
864  inline
865  bool
866  Iterator::operator != (const Iterator &other) const
867  {
868  return ! (*this == other);
869  }
870 
871 
872  inline
873  bool
874  Iterator::operator < (const Iterator &other) const
875  {
876  return accessor < other.accessor;
877  }
878 
879 
880  inline
881  int
882  Iterator::operator - (const Iterator &other) const
883  {
884  (void)other;
886  ExcInternalError());
887  Assert (false, ExcNotImplemented());
888 
889  return 0;
890  }
891 }
892 
893 
894 inline
895 void
897 {
898  // first check the last element (or if line is still empty)
899  if ( (entries.size()==0) || ( entries.back() < j) )
900  {
901  entries.push_back(j);
902  return;
903  }
904 
905  // do a binary search to find the place where to insert:
906  std::vector<size_type>::iterator
907  it = Utilities::lower_bound(entries.begin(),
908  entries.end(),
909  j);
910 
911  // If this entry is a duplicate, exit immediately
912  if (*it == j)
913  return;
914 
915  // Insert at the right place in the vector. Vector grows automatically to
916  // fit elements. Always doubles its size.
917  entries.insert(it, j);
918 }
919 
920 
921 
922 inline
925 {
926  return rows;
927 }
928 
929 
930 
931 inline
934 {
935  return cols;
936 }
937 
938 
939 
940 inline
941 void
942 DynamicSparsityPattern::add (const size_type i,
943  const size_type j)
944 {
945  Assert (i<rows, ExcIndexRangeType<size_type>(i, 0, rows));
946  Assert (j<cols, ExcIndexRangeType<size_type>(j, 0, cols));
947 
948  if (rowset.size() > 0 && !rowset.is_element(i))
949  return;
950 
951  const size_type rowindex =
952  rowset.size()==0 ? i : rowset.index_within_set(i);
953  lines[rowindex].add (j);
954 }
955 
956 
957 
958 template <typename ForwardIterator>
959 inline
960 void
962  ForwardIterator begin,
963  ForwardIterator end,
964  const bool indices_are_sorted)
965 {
966  Assert (row < rows, ExcIndexRangeType<size_type> (row, 0, rows));
967 
968  if (rowset.size() > 0 && !rowset.is_element(row))
969  return;
970 
971  const size_type rowindex =
972  rowset.size()==0 ? row : rowset.index_within_set(row);
973  lines[rowindex].add_entries (begin, end, indices_are_sorted);
974 }
975 
976 
977 
978 inline
980 {}
981 
982 
983 
984 inline
986 DynamicSparsityPattern::row_length (const size_type row) const
987 {
988  Assert (row < n_rows(), ExcIndexRangeType<size_type> (row, 0, n_rows()));
989  if (rowset.size() > 0 && !rowset.is_element(row))
990  return 0;
991 
992  const size_type rowindex =
993  rowset.size()==0 ? row : rowset.index_within_set(row);
994  return lines[rowindex].entries.size();
995 }
996 
997 
998 
999 inline
1002  const size_type index) const
1003 {
1004  Assert (row < n_rows(), ExcIndexRangeType<size_type> (row, 0, n_rows()));
1005  Assert( rowset.size() == 0 || rowset.is_element(row), ExcInternalError());
1006 
1007  const size_type local_row = rowset.size() ? rowset.index_within_set(row) : row;
1008  Assert (index < lines[local_row].entries.size(),
1009  ExcIndexRangeType<size_type> (index, 0, lines[local_row].entries.size()));
1010  return lines[local_row].entries[index];
1011 }
1012 
1013 
1014 
1015 inline
1018 {
1019  return begin(0);
1020 }
1021 
1022 
1023 inline
1026 {
1027  return iterator(this);
1028 }
1029 
1030 
1031 
1032 inline
1034 DynamicSparsityPattern::begin (const size_type r) const
1035 {
1036  Assert (r<n_rows(), ExcIndexRangeType<size_type>(r,0,n_rows()));
1037 
1038  // find the first row starting at r that has entries and return the
1039  // begin iterator to it. also skip rows for which we do not have
1040  // store anything based on the IndexSet given to the sparsity
1041  // pattern
1042  //
1043  // note: row_length(row) returns zero if the row is not locally stored
1044  //
1045  // TODO: this is way too slow when used in parallel, so do not use it on
1046  // non-owned rows
1047  size_type row = r;
1048  while ((row<n_rows())
1049  &&
1050  (row_length(row)==0))
1051  ++row;
1052 
1053  if (row == n_rows())
1054  return iterator(this);
1055  else
1056  return iterator(this, row, 0);
1057 }
1058 
1059 
1060 
1061 inline
1063 DynamicSparsityPattern::end (const size_type r) const
1064 {
1065  Assert (r<n_rows(), ExcIndexRangeType<size_type>(r,0,n_rows()));
1066 
1067  // find the first row after r that has entries and return the begin
1068  // iterator to it. also skip rows for which we do not have
1069  // store anything based on the IndexSet given to the sparsity
1070  // pattern
1071  //
1072  // note: row_length(row) returns zero if the row is not locally stored
1073  unsigned int row = r+1;
1074  while ((row<n_rows())
1075  &&
1076  (row_length(row)==0))
1077  ++row;
1078 
1079  if (row == n_rows())
1080  return iterator(this);
1081  else
1082  return iterator(this, row, 0);
1083 }
1084 
1085 
1086 
1087 inline
1088 const IndexSet &
1090 {
1091  return rowset;
1092 }
1093 
1094 
1095 
1096 inline
1097 bool
1099 {
1100  return true;
1101 }
1102 
1103 
1104 DEAL_II_NAMESPACE_CLOSE
1105 
1106 #endif
size_type row_length(const size_type row) const
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:604
const types::global_dof_index invalid_size_type
Definition: types.h:173
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
std::vector< size_type >::const_iterator end_of_row
DynamicSparsityPatternIterators::Iterator const_iterator
bool exists(const size_type i, const size_type j) const
DynamicSparsityPatternIterators::Iterator iterator
Iterator(const DynamicSparsityPattern *sp, const size_type row, const unsigned int index_within_row)
Accessor(const DynamicSparsityPattern *sparsity_pattern, const size_type row, const unsigned int index_within_row)
void add(const size_type i, const size_type j)
::ExceptionBase & ExcMessage(std::string arg1)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1081
void print(std::ostream &out) const
size_type size() const
Definition: index_set.h:1247
types::global_dof_index size_type
const IndexSet & row_index_set() const
unsigned int global_dof_index
Definition: types.h:88
size_type memory_consumption() const
#define Assert(cond, exc)
Definition: exceptions.h:294
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)
void print_gnuplot(std::ostream &out) const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
size_type max_entries_per_row() const
size_type column_number(const size_type row, const size_type index) const
std::vector< size_type >::const_iterator current_entry
DynamicSparsityPattern & operator=(const DynamicSparsityPattern &)
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1475
bool is_element(const size_type index) const
Definition: index_set.h:1317
void add(const size_type col_num)
size_type n_nonzero_elements() const