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