Reference documentation for deal.II version GIT 8e09676776 2023-03-27 21:15:01+00:00
\(\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\}}\)
dynamic_sparsity_pattern.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 2022 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 
17 #include <deal.II/base/utilities.h>
18 
21 
22 #include <algorithm>
23 #include <cmath>
24 #include <functional>
25 #include <numeric>
26 #include <set>
27 
29 
30 
31 
32 template <typename ForwardIterator>
33 void
35  ForwardIterator end,
36  const bool indices_are_sorted)
37 {
38  const int n_elements = end - begin;
39  if (n_elements <= 0)
40  return;
41 
42  const size_type stop_size = entries.size() + n_elements;
43 
44  if (indices_are_sorted == true && n_elements > 3)
45  {
46  // in debug mode, check whether the
47  // indices really are sorted.
48 #ifdef DEBUG
49  {
50  ForwardIterator test = begin, test1 = begin;
51  ++test1;
52  for (; test1 != end; ++test, ++test1)
53  Assert(*test1 > *test, ExcInternalError());
54  }
55 #endif
56 
57  if (entries.size() == 0 || entries.back() < *begin)
58  {
59  entries.insert(entries.end(), begin, end);
60  return;
61  }
62 
63  // find a possible insertion point for
64  // the first entry. check whether the
65  // first entry is a duplicate before
66  // actually doing something.
67  ForwardIterator my_it = begin;
68  size_type col = *my_it;
69  std::vector<size_type>::iterator it =
70  Utilities::lower_bound(entries.begin(), entries.end(), col);
71  while (*it == col)
72  {
73  ++my_it;
74  if (my_it == end)
75  break;
76  col = *my_it;
77  // check the very next entry in the
78  // current array
79  ++it;
80  if (it == entries.end())
81  break;
82  if (*it > col)
83  break;
84  if (*it == col)
85  continue;
86  // ok, it wasn't the very next one, do a
87  // binary search to find the insert point
88  it = Utilities::lower_bound(it, entries.end(), col);
89  if (it == entries.end())
90  break;
91  }
92  // all input entries were duplicates.
93  if (my_it == end)
94  return;
95 
96  // resize vector by just inserting the
97  // list
98  const size_type pos1 = it - entries.begin();
99  Assert(pos1 <= entries.size(), ExcInternalError());
100  entries.insert(it, my_it, end);
101  it = entries.begin() + pos1;
102  Assert(entries.size() >= static_cast<size_type>(it - entries.begin()),
103  ExcInternalError());
104 
105  // now merge the two lists.
106  std::vector<size_type>::iterator it2 = it + (end - my_it);
107 
108  // as long as there are indices both in
109  // the end of the entries list and in the
110  // input list
111  while (my_it != end && it2 != entries.end())
112  {
113  if (*my_it < *it2)
114  *it++ = *my_it++;
115  else if (*my_it == *it2)
116  {
117  *it++ = *it2++;
118  ++my_it;
119  }
120  else
121  *it++ = *it2++;
122  }
123  // in case there are indices left in the
124  // input list
125  while (my_it != end)
126  *it++ = *my_it++;
127 
128  // in case there are indices left in the
129  // end of entries
130  while (it2 != entries.end())
131  *it++ = *it2++;
132 
133  // resize and return
134  const size_type new_size = it - entries.begin();
135  Assert(new_size <= stop_size, ExcInternalError());
136  entries.resize(new_size);
137  return;
138  }
139 
140  // unsorted case or case with too few
141  // elements
142  ForwardIterator my_it = begin;
143 
144  // If necessary, increase the size of the
145  // array.
146  if (stop_size > entries.capacity())
147  entries.reserve(stop_size);
148 
149  size_type col = *my_it;
150  std::vector<size_type>::iterator it, it2;
151  // insert the first element as for one
152  // entry only first check the last
153  // element (or if line is still empty)
154  if ((entries.size() == 0) || (entries.back() < col))
155  {
156  entries.push_back(col);
157  it = entries.end() - 1;
158  }
159  else
160  {
161  // do a binary search to find the place
162  // where to insert:
163  it2 = Utilities::lower_bound(entries.begin(), entries.end(), col);
164 
165  // If this entry is a duplicate, continue
166  // immediately Insert at the right place
167  // in the vector. Vector grows
168  // automatically to fit elements. Always
169  // doubles its size.
170  if (*it2 != col)
171  it = entries.insert(it2, col);
172  else
173  it = it2;
174  }
175 
176  ++my_it;
177  // Now try to be smart and insert with
178  // bias in the direction we are
179  // walking. This has the advantage that
180  // for sorted lists, we always search in
181  // the right direction, what should
182  // decrease the work needed in here.
183  for (; my_it != end; ++my_it)
184  {
185  col = *my_it;
186  // need a special insertion command when
187  // we're at the end of the list
188  if (col > entries.back())
189  {
190  entries.push_back(col);
191  it = entries.end() - 1;
192  }
193  // search to the right (preferred search
194  // direction)
195  else if (col > *it)
196  {
197  it2 = Utilities::lower_bound(it++, entries.end(), col);
198  if (*it2 != col)
199  it = entries.insert(it2, col);
200  }
201  // search to the left
202  else if (col < *it)
203  {
204  it2 = Utilities::lower_bound(entries.begin(), it, col);
205  if (*it2 != col)
206  it = entries.insert(it2, col);
207  }
208  // if we're neither larger nor smaller,
209  // then this was a duplicate and we can
210  // just continue.
211  }
212 }
213 
214 
217 {
218  return entries.capacity() * sizeof(size_type) + sizeof(Line);
219 }
220 
221 
224  , have_entries(false)
225  , rowset(0)
226 {}
227 
228 
229 
232  , have_entries(false)
233  , rowset(0)
234 {
235  (void)s;
236  Assert(s.rows == 0 && s.cols == 0,
237  ExcMessage(
238  "This constructor can only be called if the provided argument "
239  "is the sparsity pattern for an empty matrix. This constructor can "
240  "not be used to copy-construct a non-empty sparsity pattern."));
241 }
242 
243 
244 
246  const size_type n,
247  const IndexSet &rowset_)
249  , have_entries(false)
250  , rowset(0)
251 {
252  reinit(m, n, rowset_);
253 }
254 
255 
258  , have_entries(false)
259  , rowset(0)
260 {
261  reinit(rowset_.size(), rowset_.size(), rowset_);
262 }
263 
264 
267  , have_entries(false)
268  , rowset(0)
269 {
270  reinit(n, n);
271 }
272 
273 
274 
277 {
278  (void)s;
279  Assert(s.n_rows() == 0 && s.n_cols() == 0,
280  ExcMessage(
281  "This operator can only be called if the provided argument "
282  "is the sparsity pattern for an empty matrix. This operator can "
283  "not be used to copy a non-empty sparsity pattern."));
284 
285  Assert(n_rows() == 0 && n_cols() == 0,
286  ExcMessage("This operator can only be called if the current object is "
287  "empty."));
288 
289  return *this;
290 }
291 
292 
293 
294 void
296  const size_type n,
297  const IndexSet &rowset_)
298 {
299  resize(m, n);
300  have_entries = false;
301  rowset = rowset_;
302 
303  Assert(rowset.size() == 0 || rowset.size() == m,
304  ExcMessage(
305  "The IndexSet argument to this function needs to either "
306  "be empty (indicating the complete set of rows), or have size "
307  "equal to the desired number of rows as specified by the "
308  "first argument to this function. (Of course, the number "
309  "of indices in this IndexSet may be less than the number "
310  "of rows, but the *size* of the IndexSet must be equal.)"));
311 
312  std::vector<Line> new_lines(rowset.size() == 0 ? n_rows() :
313  rowset.n_elements());
314  lines.swap(new_lines);
315 }
316 
317 
318 
319 void
321 {}
322 
323 
324 
325 bool
327 {
328  return ((rows == 0) && (cols == 0));
329 }
330 
331 
332 
335 {
336  if (!have_entries)
337  return 0;
338 
339  size_type m = 0;
340  for (const auto &line : lines)
341  {
342  m = std::max(m, static_cast<size_type>(line.entries.size()));
343  }
344 
345  return m;
346 }
347 
348 
349 
350 void
352  const size_type & row,
353  const ArrayView<const size_type> &columns,
354  const bool indices_are_sorted)
355 {
356  add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
357 }
358 
359 
360 
361 bool
363 {
366  Assert(
367  rowset.size() == 0 || rowset.is_element(i),
368  ExcMessage(
369  "The row IndexSet does not contain the index i. This sparsity pattern "
370  "object cannot know whether the entry (i, j) exists or not."));
371 
372  // Avoid a segmentation fault in below code if the row index happens to
373  // not be present in the IndexSet rowset:
374  if (!(rowset.size() == 0 || rowset.is_element(i)))
375  return false;
376 
377  if (!have_entries)
378  return false;
379 
380  const size_type rowindex =
381  rowset.size() == 0 ? i : rowset.index_within_set(i);
382 
383  return std::binary_search(lines[rowindex].entries.begin(),
384  lines[rowindex].entries.end(),
385  j);
386 }
387 
388 
389 
390 void
392 {
394 
395  // loop over all elements presently in the sparsity pattern and add the
396  // transpose element. note:
397  //
398  // 1. that the sparsity pattern changes which we work on, but not the present
399  // row
400  //
401  // 2. that the @p{add} function can be called on elements that already exist
402  // without any harm
403  for (size_type row = 0; row < lines.size(); ++row)
404  {
405  const size_type rowindex =
406  rowset.size() == 0 ? row : rowset.nth_index_in_set(row);
407 
408  for (const size_type row_entry : lines[row].entries)
409  // add the transpose entry if this is not the diagonal
410  if (rowindex != row_entry)
411  add(row_entry, rowindex);
412  }
413 }
414 
415 
416 
417 void
419 {
420  AssertIndexRange(row, n_rows());
421  if (!have_entries)
422  return;
423 
424  if (rowset.size() > 0 && !rowset.is_element(row))
425  return;
426 
427  const size_type rowindex =
428  rowset.size() == 0 ? row : rowset.index_within_set(row);
429 
430  AssertIndexRange(rowindex, lines.size());
431  lines[rowindex].entries = std::vector<size_type>();
432 }
433 
434 
435 
438 {
440  view.reinit(rows.n_elements(), this->n_cols());
441  AssertDimension(rows.size(), this->n_rows());
442 
443  const auto end = rows.end();
445  for (auto it = rows.begin(); it != end; ++it, ++view_row)
446  {
447  const size_type rowindex =
448  rowset.size() == 0 ? *it : rowset.index_within_set(*it);
449 
450  view.lines[view_row].entries = lines[rowindex].entries;
451  view.have_entries |= (lines[rowindex].entries.size() > 0);
452  }
453  return view;
454 }
455 
456 
457 
458 template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
459 void
461  const SparsityPatternTypeLeft & sp_A,
462  const SparsityPatternTypeRight &sp_B)
463 {
464  Assert(sp_A.n_rows() == sp_B.n_rows(),
465  ExcDimensionMismatch(sp_A.n_rows(), sp_B.n_rows()));
466 
467  this->reinit(sp_A.n_cols(), sp_B.n_cols());
468  // we will go through all the
469  // rows in the matrix A, and for each column in a row we add the whole
470  // row of matrix B with that row number. This means that we will insert
471  // a lot of entries to each row, which is best handled by the
472  // DynamicSparsityPattern class.
473 
474  std::vector<size_type> new_cols;
475  new_cols.reserve(sp_B.max_entries_per_row());
476 
477  // C_{kl} = A_{ik} B_{il}
478  for (size_type i = 0; i < sp_A.n_rows(); ++i)
479  {
480  // get all column numbers from sp_B in a temporary vector:
481  new_cols.resize(sp_B.row_length(i));
482  {
483  const auto last_il = sp_B.end(i);
484  auto * col_ptr = new_cols.data();
485  for (auto il = sp_B.begin(i); il != last_il; ++il)
486  *col_ptr++ = il->column();
487  }
488  std::sort(new_cols.begin(), new_cols.end());
489 
490  // now for each k, add new_cols to the target sparsity
491  const auto last_ik = sp_A.end(i);
492  for (auto ik = sp_A.begin(i); ik != last_ik; ++ik)
493  this->add_entries(ik->column(), new_cols.begin(), new_cols.end(), true);
494  }
495 }
496 
497 
498 
499 template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
500 void
502  const SparsityPatternTypeLeft & left,
503  const SparsityPatternTypeRight &right)
504 {
505  Assert(left.n_cols() == right.n_rows(),
506  ExcDimensionMismatch(left.n_cols(), right.n_rows()));
507 
508  this->reinit(left.n_rows(), right.n_cols());
509 
510  typename SparsityPatternTypeLeft::iterator it_left = left.begin(),
511  end_left = left.end();
512  for (; it_left != end_left; ++it_left)
513  {
514  const unsigned int j = it_left->column();
515 
516  // We are sitting on entry (i,j) of the left sparsity pattern. We then
517  // need to add all entries (i,k) to the final sparsity pattern where (j,k)
518  // exists in the right sparsity pattern -- i.e., we need to iterate over
519  // row j.
520  typename SparsityPatternTypeRight::iterator it_right = right.begin(j),
521  end_right = right.end(j);
522  for (; it_right != end_right; ++it_right)
523  this->add(it_left->row(), it_right->column());
524  }
525 }
526 
527 
528 
529 void
530 DynamicSparsityPattern::print(std::ostream &out) const
531 {
532  for (size_type row = 0; row < lines.size(); ++row)
533  {
534  out << '[' << (rowset.size() == 0 ? row : rowset.nth_index_in_set(row));
535 
536  for (const auto entry : lines[row].entries)
537  out << ',' << entry;
538 
539  out << ']' << std::endl;
540  }
541 
542  AssertThrow(out.fail() == false, ExcIO());
543 }
544 
545 
546 
547 void
549 {
550  for (size_type row = 0; row < lines.size(); ++row)
551  {
552  const size_type rowindex =
553  rowset.size() == 0 ? row : rowset.nth_index_in_set(row);
554 
555  for (const auto entry : lines[row].entries)
556  // while matrix entries are usually
557  // written (i,j), with i vertical and
558  // j horizontal, gnuplot output is
559  // x-y, that is we have to exchange
560  // the order of output
561  out << entry << " " << -static_cast<signed int>(rowindex) << std::endl;
562  }
563 
564 
565  AssertThrow(out.fail() == false, ExcIO());
566 }
567 
568 
569 
572 {
573  size_type b = 0;
574  for (size_type row = 0; row < lines.size(); ++row)
575  {
576  const size_type rowindex =
577  rowset.size() == 0 ? row : rowset.nth_index_in_set(row);
578 
579  for (const auto entry : lines[row].entries)
580  if (static_cast<size_type>(
581  std::abs(static_cast<int>(rowindex - entry))) > b)
582  b = std::abs(static_cast<signed int>(rowindex - entry));
583  }
584 
585  return b;
586 }
587 
588 
589 
592 {
593  if (!have_entries)
594  return 0;
595 
596  size_type n = 0;
597  for (const auto &line : lines)
598  {
599  n += line.entries.size();
600  }
601 
602  return n;
603 }
604 
605 
606 
607 IndexSet
609 {
610  std::set<types::global_dof_index> cols;
611  for (const auto &line : lines)
612  cols.insert(line.entries.begin(), line.entries.end());
613 
614  IndexSet res(this->n_cols());
615  res.add_indices(cols.begin(), cols.end());
616  return res;
617 }
618 
619 
620 
621 IndexSet
623 {
624  const IndexSet all_rows = complete_index_set(this->n_rows());
625  const IndexSet &locally_stored_rows = rowset.size() == 0 ? all_rows : rowset;
626 
627  std::vector<types::global_dof_index> rows;
628  auto line = lines.begin();
629  AssertDimension(locally_stored_rows.n_elements(), lines.size());
630  for (const auto row : locally_stored_rows)
631  {
632  if (line->entries.size() > 0)
633  rows.push_back(row);
634 
635  ++line;
636  }
637 
638  IndexSet res(this->n_rows());
639  res.add_indices(rows.begin(), rows.end());
640  return res;
641 }
642 
643 
644 
647 {
648  size_type mem = sizeof(DynamicSparsityPattern) +
650  sizeof(rowset);
651 
652  for (const auto &line : lines)
654 
655  return mem;
656 }
657 
658 
659 
663  const DynamicSparsityPattern::size_type col) const
664 {
665  AssertIndexRange(row, n_rows());
666  AssertIndexRange(col, n_cols());
667  Assert(rowset.size() == 0 || rowset.is_element(row), ExcInternalError());
668 
669  const DynamicSparsityPattern::size_type local_row =
670  rowset.size() != 0u ? rowset.index_within_set(row) : row;
671 
672  // now we need to do a binary search. Note that col indices are assumed to
673  // be sorted.
674  const auto &cols = lines[local_row].entries;
675  auto it = Utilities::lower_bound(cols.begin(), cols.end(), col);
676 
677  if ((it != cols.end()) && (*it == col))
678  return (it - cols.begin());
679  else
681 }
682 
683 
684 
685 // explicit instantiations
686 template void
688 template void
690  const size_type *,
691  const bool);
692 #ifndef DEAL_II_VECTOR_ITERATOR_IS_POINTER
693 template void
694 DynamicSparsityPattern::Line::add_entries(std::vector<size_type>::iterator,
695  std::vector<size_type>::iterator,
696  const bool);
697 template void
699  std::vector<size_type>::const_iterator,
700  std::vector<size_type>::const_iterator,
701  const bool);
702 #endif
703 
704 template void
706  const DynamicSparsityPattern &);
707 template void
709  const SparsityPattern &);
710 template void
712  const DynamicSparsityPattern &);
713 template void
715  const SparsityPattern &);
716 
717 template void
719  const SparsityPattern &);
720 template void
722  const SparsityPattern &);
723 template void
725  const DynamicSparsityPattern &);
726 template void
728  const DynamicSparsityPattern &);
729 
iterator begin() const
Definition: array_view.h:594
iterator end() const
Definition: array_view.h:603
DynamicSparsityPattern get_view(const IndexSet &rows) const
types::global_dof_index size_type
void compute_mmult_pattern(const SparsityPatternTypeLeft &left, const SparsityPatternTypeRight &right)
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
size_type column_index(const size_type row, const size_type col) const
DynamicSparsityPattern & operator=(const DynamicSparsityPattern &)
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
void print(std::ostream &out) const
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
void clear_row(const size_type row)
void print_gnuplot(std::ostream &out) const
void compute_Tmmult_pattern(const SparsityPatternTypeLeft &left, const SparsityPatternTypeRight &right)
bool exists(const size_type i, const size_type j) const
void add(const size_type i, const size_type j)
size_type size() const
Definition: index_set.h:1648
size_type index_within_set(const size_type global_index) const
Definition: index_set.h:1863
size_type n_elements() const
Definition: index_set.h:1796
bool is_element(const size_type index) const
Definition: index_set.h:1756
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1844
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1714
virtual void resize(const size_type rows, const size_type cols)
size_type n_rows() const
size_type n_cols() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1586
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1759
#define AssertIndexRange(index, range)
Definition: exceptions.h:1827
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
IndexSet complete_index_set(const IndexSet::size_type N)
Definition: index_set.h:1076
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:999
const types::global_dof_index invalid_size_type
Definition: types.h:222
unsigned int global_dof_index
Definition: types.h:82
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)