Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
dynamic_sparsity_pattern.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2008 - 2023 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
18
21
22#include <algorithm>
23#include <cmath>
24#include <functional>
25#include <numeric>
26#include <set>
27
29
30
31
32template <typename ForwardIterator>
33void
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()),
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,
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
257 : DynamicSparsityPattern(rowset_.size(), rowset_.size(), rowset_)
258{}
259
260
263 , have_entries(false)
264 , rowset(0)
265{
266 reinit(n, n);
267}
268
269
270
273{
274 (void)s;
275 Assert(s.n_rows() == 0 && s.n_cols() == 0,
277 "This operator can only be called if the provided argument "
278 "is the sparsity pattern for an empty matrix. This operator can "
279 "not be used to copy a non-empty sparsity pattern."));
280
281 Assert(n_rows() == 0 && n_cols() == 0,
282 ExcMessage("This operator can only be called if the current object is "
283 "empty."));
284
285 return *this;
286}
287
288
289
290void
292 const size_type n,
293 const IndexSet &rowset_)
294{
295 resize(m, n);
296 have_entries = false;
297 rowset = rowset_;
298
299 Assert(rowset.size() == 0 || rowset.size() == m,
301 "The IndexSet argument to this function needs to either "
302 "be empty (indicating the complete set of rows), or have size "
303 "equal to the desired number of rows as specified by the "
304 "first argument to this function. (Of course, the number "
305 "of indices in this IndexSet may be less than the number "
306 "of rows, but the *size* of the IndexSet must be equal.)"));
307
308 std::vector<Line> new_lines(rowset.size() == 0 ? n_rows() :
310 lines.swap(new_lines);
311}
312
313
314
315void
317{}
318
319
320
321bool
323{
324 return ((rows == 0) && (cols == 0));
325}
326
327
328
331{
332 if (!have_entries)
333 return 0;
334
335 size_type m = 0;
336 for (const auto &line : lines)
337 {
338 m = std::max(m, static_cast<size_type>(line.entries.size()));
339 }
340
341 return m;
342}
343
344
345
346void
348 const size_type & row,
349 const ArrayView<const size_type> &columns,
350 const bool indices_are_sorted)
351{
352 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
353}
354
355
356
357bool
359{
362 Assert(
363 rowset.size() == 0 || rowset.is_element(i),
365 "The row IndexSet does not contain the index i. This sparsity pattern "
366 "object cannot know whether the entry (i, j) exists or not."));
367
368 // Avoid a segmentation fault in below code if the row index happens to
369 // not be present in the IndexSet rowset:
370 if (!(rowset.size() == 0 || rowset.is_element(i)))
371 return false;
372
373 if (!have_entries)
374 return false;
375
376 const size_type rowindex =
377 rowset.size() == 0 ? i : rowset.index_within_set(i);
378
379 return std::binary_search(lines[rowindex].entries.begin(),
380 lines[rowindex].entries.end(),
381 j);
382}
383
384
385
386void
388{
390
391 // loop over all elements presently in the sparsity pattern and add the
392 // transpose element. note:
393 //
394 // 1. that the sparsity pattern changes which we work on, but not the present
395 // row
396 //
397 // 2. that the @p{add} function can be called on elements that already exist
398 // without any harm
399 for (size_type row = 0; row < lines.size(); ++row)
400 {
401 const size_type rowindex =
402 rowset.size() == 0 ? row : rowset.nth_index_in_set(row);
403
404 for (const size_type row_entry : lines[row].entries)
405 // add the transpose entry if this is not the diagonal
406 if (rowindex != row_entry)
407 add(row_entry, rowindex);
408 }
409}
410
411
412
413void
415{
416 AssertIndexRange(row, n_rows());
417 if (!have_entries)
418 return;
419
420 if (rowset.size() > 0 && !rowset.is_element(row))
421 return;
422
423 const size_type rowindex =
424 rowset.size() == 0 ? row : rowset.index_within_set(row);
425
426 AssertIndexRange(rowindex, lines.size());
427 lines[rowindex].entries = std::vector<size_type>();
428}
429
430
431
434{
436 view.reinit(rows.n_elements(), this->n_cols());
437 AssertDimension(rows.size(), this->n_rows());
438
439 const auto end = rows.end();
441 for (auto it = rows.begin(); it != end; ++it, ++view_row)
442 {
443 const size_type rowindex =
444 rowset.size() == 0 ? *it : rowset.index_within_set(*it);
445
446 view.lines[view_row].entries = lines[rowindex].entries;
447 view.have_entries |= (lines[rowindex].entries.size() > 0);
448 }
449 return view;
450}
451
452
453
454template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
455void
457 const SparsityPatternTypeLeft & sp_A,
458 const SparsityPatternTypeRight &sp_B)
459{
460 Assert(sp_A.n_rows() == sp_B.n_rows(),
461 ExcDimensionMismatch(sp_A.n_rows(), sp_B.n_rows()));
462
463 this->reinit(sp_A.n_cols(), sp_B.n_cols());
464 // we will go through all the
465 // rows in the matrix A, and for each column in a row we add the whole
466 // row of matrix B with that row number. This means that we will insert
467 // a lot of entries to each row, which is best handled by the
468 // DynamicSparsityPattern class.
469
470 std::vector<size_type> new_cols;
471 new_cols.reserve(sp_B.max_entries_per_row());
472
473 // C_{kl} = A_{ik} B_{il}
474 for (size_type i = 0; i < sp_A.n_rows(); ++i)
475 {
476 // get all column numbers from sp_B in a temporary vector:
477 new_cols.resize(sp_B.row_length(i));
478 {
479 const auto last_il = sp_B.end(i);
480 auto * col_ptr = new_cols.data();
481 for (auto il = sp_B.begin(i); il != last_il; ++il)
482 *col_ptr++ = il->column();
483 }
484 std::sort(new_cols.begin(), new_cols.end());
485
486 // now for each k, add new_cols to the target sparsity
487 const auto last_ik = sp_A.end(i);
488 for (auto ik = sp_A.begin(i); ik != last_ik; ++ik)
489 this->add_entries(ik->column(), new_cols.begin(), new_cols.end(), true);
490 }
491}
492
493
494
495template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
496void
498 const SparsityPatternTypeLeft & left,
499 const SparsityPatternTypeRight &right)
500{
501 Assert(left.n_cols() == right.n_rows(),
502 ExcDimensionMismatch(left.n_cols(), right.n_rows()));
503
504 this->reinit(left.n_rows(), right.n_cols());
505
506 typename SparsityPatternTypeLeft::iterator it_left = left.begin(),
507 end_left = left.end();
508 for (; it_left != end_left; ++it_left)
509 {
510 const unsigned int j = it_left->column();
511
512 // We are sitting on entry (i,j) of the left sparsity pattern. We then
513 // need to add all entries (i,k) to the final sparsity pattern where (j,k)
514 // exists in the right sparsity pattern -- i.e., we need to iterate over
515 // row j.
516 typename SparsityPatternTypeRight::iterator it_right = right.begin(j),
517 end_right = right.end(j);
518 for (; it_right != end_right; ++it_right)
519 this->add(it_left->row(), it_right->column());
520 }
521}
522
523
524
525void
526DynamicSparsityPattern::print(std::ostream &out) const
527{
528 for (size_type row = 0; row < lines.size(); ++row)
529 {
530 out << '[' << (rowset.size() == 0 ? row : rowset.nth_index_in_set(row));
531
532 for (const auto entry : lines[row].entries)
533 out << ',' << entry;
534
535 out << ']' << std::endl;
536 }
537
538 AssertThrow(out.fail() == false, ExcIO());
539}
540
541
542
543void
545{
546 for (size_type row = 0; row < lines.size(); ++row)
547 {
548 const size_type rowindex =
549 rowset.size() == 0 ? row : rowset.nth_index_in_set(row);
550
551 for (const auto entry : lines[row].entries)
552 // while matrix entries are usually
553 // written (i,j), with i vertical and
554 // j horizontal, gnuplot output is
555 // x-y, that is we have to exchange
556 // the order of output
557 out << entry << " " << -static_cast<signed int>(rowindex) << std::endl;
558 }
559
560
561 AssertThrow(out.fail() == false, ExcIO());
562}
563
564
565
568{
569 size_type b = 0;
570 for (size_type row = 0; row < lines.size(); ++row)
571 {
572 const size_type rowindex =
573 rowset.size() == 0 ? row : rowset.nth_index_in_set(row);
574
575 for (const auto entry : lines[row].entries)
576 if (static_cast<size_type>(
577 std::abs(static_cast<int>(rowindex - entry))) > b)
578 b = std::abs(static_cast<signed int>(rowindex - entry));
579 }
580
581 return b;
582}
583
584
585
588{
589 if (!have_entries)
590 return 0;
591
592 size_type n = 0;
593 for (const auto &line : lines)
594 {
595 n += line.entries.size();
596 }
597
598 return n;
599}
600
601
602
605{
606 std::set<types::global_dof_index> cols;
607 for (const auto &line : lines)
608 cols.insert(line.entries.begin(), line.entries.end());
609
610 IndexSet res(this->n_cols());
611 res.add_indices(cols.begin(), cols.end());
612 return res;
613}
614
615
616
619{
620 const IndexSet all_rows = complete_index_set(this->n_rows());
621 const IndexSet &locally_stored_rows = rowset.size() == 0 ? all_rows : rowset;
622
623 std::vector<types::global_dof_index> rows;
624 auto line = lines.begin();
625 AssertDimension(locally_stored_rows.n_elements(), lines.size());
626 for (const auto row : locally_stored_rows)
627 {
628 if (line->entries.size() > 0)
629 rows.push_back(row);
630
631 ++line;
632 }
633
634 IndexSet res(this->n_rows());
635 res.add_indices(rows.begin(), rows.end());
636 return res;
637}
638
639
640
643{
644 size_type mem = sizeof(DynamicSparsityPattern) +
646 sizeof(rowset);
647
648 for (const auto &line : lines)
650
651 return mem;
652}
653
654
655
659 const DynamicSparsityPattern::size_type col) const
660{
661 AssertIndexRange(row, n_rows());
662 AssertIndexRange(col, n_cols());
664
665 const DynamicSparsityPattern::size_type local_row =
666 rowset.size() != 0u ? rowset.index_within_set(row) : row;
667
668 // now we need to do a binary search. Note that col indices are assumed to
669 // be sorted.
670 const auto &cols = lines[local_row].entries;
671 auto it = Utilities::lower_bound(cols.begin(), cols.end(), col);
672
673 if ((it != cols.end()) && (*it == col))
674 return (it - cols.begin());
675 else
677}
678
679
680
681// explicit instantiations
682template void
683DynamicSparsityPattern::Line::add_entries(size_type *, size_type *, const bool);
684template void
686 const size_type *,
687 const bool);
688#ifndef DEAL_II_VECTOR_ITERATOR_IS_POINTER
689template void
690DynamicSparsityPattern::Line::add_entries(std::vector<size_type>::iterator,
691 std::vector<size_type>::iterator,
692 const bool);
693template void
695 std::vector<size_type>::const_iterator,
696 std::vector<size_type>::const_iterator,
697 const bool);
698#endif
699
700template void
702 const DynamicSparsityPattern &);
703template void
705 const SparsityPattern &);
706template void
708 const DynamicSparsityPattern &);
709template void
711 const SparsityPattern &);
712
713template void
715 const SparsityPattern &);
716template void
718 const SparsityPattern &);
719template void
721 const DynamicSparsityPattern &);
722template void
724 const DynamicSparsityPattern &);
725
iterator begin() const
Definition array_view.h:594
iterator end() const
Definition array_view.h:603
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_unique_and_sorted=false)
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 &)
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:1661
size_type index_within_set(const size_type global_index) const
Definition index_set.h:1883
size_type n_elements() const
Definition index_set.h:1816
bool is_element(const size_type index) const
Definition index_set.h:1776
size_type nth_index_in_set(const size_type local_index) const
Definition index_set.h:1864
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition index_set.h:1727
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
IndexSet complete_index_set(const IndexSet::size_type N)
Definition index_set.h:1089
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1016
const types::global_dof_index invalid_size_type
Definition types.h:222
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
void add_entries(ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted)