deal.II version GIT relicensing-3540-g7552a02177 2025-06-20 13:50:00+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
sparsity_pattern.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 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
15
17
22
23#include <algorithm>
24#include <cmath>
25#include <functional>
26#include <iomanip>
27#include <iostream>
28#include <memory>
29#include <numeric>
30
32
34
35
36
39 , store_diagonal_first_in_row(false)
40 , max_dim(0)
41 , max_vec_len(0)
42 , max_row_length(0)
43 , compressed(false)
44{
45 reinit(0, 0, 0);
46}
47
48
49
52{
53 Assert(s.empty(),
55 "This constructor can only be called if the provided argument "
56 "is the sparsity pattern for an empty matrix. This constructor can "
57 "not be used to copy-construct a non-empty sparsity pattern."));
58
59 reinit(0, 0, 0);
60}
61
62
63
65 const size_type n,
66 const unsigned int max_per_row)
68{
69 reinit(m, n, max_per_row);
70}
71
72
73
75 const size_type n,
76 const std::vector<unsigned int> &row_lengths)
78{
79 reinit(m, n, row_lengths);
80}
81
82
83
85 const unsigned int max_per_row)
87{
88 reinit(m, m, max_per_row);
89}
90
91
92
94 const std::vector<unsigned int> &row_lengths)
96{
97 reinit(m, m, row_lengths);
98}
99
100
101
103 const unsigned int max_per_row,
104 const size_type extra_off_diagonals)
106{
107 Assert(original.n_rows() == original.n_cols(), ExcNotQuadratic());
109
110 reinit(original.n_rows(), original.n_cols(), max_per_row);
111
112 // now copy the entries from the other object
113 for (size_type row = 0; row < original.rows; ++row)
114 {
115 // copy the elements of this row of the other object
116 //
117 // note that the first object actually is the main-diagonal element,
118 // which we need not copy
119 //
120 // we do the copying in two steps: first we note that the elements in
121 // @p{original} are sorted, so we may first copy all the elements up to
122 // the first side-diagonal one which is to be filled in. then we insert
123 // the side-diagonals, finally copy the rest from that element onwards
124 // which is not a side-diagonal any more.
125 const size_type *const original_row_start =
126 &original.colnums[original.rowstart[row]] + 1;
127 // the following requires that @p{original} be compressed since
128 // otherwise there might be invalid_entry's
129 const size_type *const original_row_end =
130 &original.colnums[original.rowstart[row + 1]];
131
132 // find pointers before and after extra off-diagonals. if at top or
133 // bottom of matrix, then set these pointers such that no copying is
134 // necessary (see the @p{copy} commands)
135 const size_type *const original_last_before_side_diagonals =
136 (row > extra_off_diagonals ?
137 Utilities::lower_bound(original_row_start,
138 original_row_end,
139 row - extra_off_diagonals) :
140 original_row_start);
141
142 const size_type *const original_first_after_side_diagonals =
143 (row < rows - extra_off_diagonals - 1 ?
144 std::upper_bound(original_row_start,
145 original_row_end,
146 row + extra_off_diagonals) :
147 original_row_end);
148
149 // find first free slot. the first slot in each row is the diagonal
150 // element
151 size_type *next_free_slot = &colnums[rowstart[row]] + 1;
152
153 // copy elements before side-diagonals
154 next_free_slot = std::copy(original_row_start,
155 original_last_before_side_diagonals,
156 next_free_slot);
157
158 // insert left and right side-diagonals
159 for (size_type i = 1; i <= std::min(row, extra_off_diagonals);
160 ++i, ++next_free_slot)
161 *next_free_slot = row - i;
162 for (size_type i = 1; i <= std::min(extra_off_diagonals, rows - row - 1);
163 ++i, ++next_free_slot)
164 *next_free_slot = row + i;
165
166 // copy rest
167 next_free_slot = std::copy(original_first_after_side_diagonals,
168 original_row_end,
169 next_free_slot);
170
171 // this error may happen if the sum of previous elements per row and
172 // those of the new diagonals exceeds the maximum number of elements per
173 // row given to this constructor
174 Assert(next_free_slot <= &colnums[rowstart[row + 1]],
175 ExcNotEnoughSpace(0, rowstart[row + 1] - rowstart[row]));
176 }
177}
178
179
180
183{
184 Assert(s.empty(),
186 "This operator can only be called if the provided argument "
187 "is the sparsity pattern for an empty matrix. This operator can "
188 "not be used to copy a non-empty sparsity pattern."));
189
190 Assert(this->empty(),
191 ExcMessage("This operator can only be called if the current object is "
192 "empty."));
193
194 return *this;
195}
196
197
198
199void
201 const size_type n,
202 const ArrayView<const unsigned int> &row_lengths)
203{
204 AssertDimension(row_lengths.size(), m);
205 resize(m, n);
206
207 // delete empty matrices
208 if ((m == 0) || (n == 0))
209 {
210 rowstart.reset();
211 colnums.reset();
212
213 max_vec_len = max_dim = 0;
214 // if dimension is zero: ignore max_per_row
215 max_row_length = 0;
216 compressed = false;
217
218 return;
219 }
220
221 // first, if the matrix is quadratic, we will have to make sure that each
222 // row has at least one entry for the diagonal element. make this more
223 // obvious by having a variable which we can query
225
226 // find out how many entries we need in the @p{colnums} array. if this
227 // number is larger than @p{max_vec_len}, then we will need to reallocate
228 // memory
229 //
230 // note that the number of elements per row is bounded by the number of
231 // columns
232 //
233 std::size_t vec_len = 0;
234 for (size_type i = 0; i < m; ++i)
235 vec_len += std::min(static_cast<size_type>(store_diagonal_first_in_row ?
236 std::max(row_lengths[i], 1U) :
237 row_lengths[i]),
238 n);
239
240 // sometimes, no entries are requested in the matrix (this most often
241 // happens when blocks in a block matrix are simply zero). in that case,
242 // allocate exactly one element, to have a valid pointer to some memory
243 if (vec_len == 0)
244 {
245 vec_len = 1;
246 max_vec_len = vec_len;
247 colnums = std::make_unique<size_type[]>(max_vec_len);
248 }
249
251 (row_lengths.empty() ?
252 0 :
253 std::min(static_cast<size_type>(
254 *std::max_element(row_lengths.begin(), row_lengths.end())),
255 n));
256
257 if (store_diagonal_first_in_row && (max_row_length == 0) && (m != 0))
258 max_row_length = 1;
259
260 // allocate memory for the rowstart values, if necessary. even though we
261 // re-set the pointers again immediately after deleting their old content,
262 // set them to zero in between because the allocation might fail, in which
263 // case we get an exception and the destructor of this object will be called
264 // -- where we look at the non-nullness of the (now invalid) pointer again
265 // and try to delete the memory a second time.
266 if (rows > max_dim)
267 {
268 max_dim = rows;
269 rowstart = std::make_unique<std::size_t[]>(max_dim + 1);
270 }
271
272 // allocate memory for the column numbers if necessary
273 if (vec_len > max_vec_len)
274 {
275 max_vec_len = vec_len;
276 colnums = std::make_unique<size_type[]>(max_vec_len);
277 }
278
279 // set the rowstart array
280 rowstart[0] = 0;
281 for (size_type i = 1; i <= rows; ++i)
282 rowstart[i] =
283 rowstart[i - 1] +
285 std::max(std::min(static_cast<size_type>(row_lengths[i - 1]), n),
286 static_cast<size_type>(1U)) :
287 std::min(static_cast<size_type>(row_lengths[i - 1]), n));
288 Assert((rowstart[rows] == vec_len) ||
289 ((vec_len == 1) && (rowstart[rows] == 0)),
291
292 // preset the column numbers by a value indicating it is not in use
293 std::fill_n(colnums.get(), vec_len, invalid_entry);
294
295 // if diagonal elements are special: let the first entry in each row be the
296 // diagonal value
298 for (size_type i = 0; i < n_rows(); ++i)
299 colnums[rowstart[i]] = i;
300
301 compressed = false;
302}
303
304
305
306void
308 const size_type n,
309 const std::vector<unsigned int> &row_lengths)
310{
311 reinit(m, n, make_array_view(row_lengths));
312}
313
314
315
316void
318 const size_type n,
319 const unsigned int max_per_row)
320{
321 // simply map this function to the other @p{reinit} function
322 const std::vector<unsigned int> row_lengths(m, max_per_row);
323 reinit(m, n, row_lengths);
324}
325
326
327
328void
330{
331 // nothing to do if the object corresponds to an empty matrix
332 if ((rowstart == nullptr) && (colnums == nullptr))
333 {
334 compressed = true;
335 return;
336 }
337
338 // do nothing if already compressed
339 if (compressed)
340 return;
341
342 std::size_t next_free_entry = 0, next_row_start = 0, row_length = 0;
343
344 // first find out how many non-zero elements there are, in order to allocate
345 // the right amount of memory
346 const std::size_t nonzero_elements =
347 std::count_if(&colnums[rowstart[0]],
349 [](const size_type col) { return col != invalid_entry; });
350 // now allocate the respective memory
351 std::unique_ptr<size_type[]> new_colnums(new size_type[nonzero_elements]);
352
353
354 // reserve temporary storage to store the entries of one row
355 std::vector<size_type> tmp_entries(max_row_length);
356
357 // Traverse all rows
358 for (size_type line = 0; line < n_rows(); ++line)
359 {
360 // copy used entries, break if first unused entry is reached
361 row_length = 0;
362 for (std::size_t j = rowstart[line]; j < rowstart[line + 1];
363 ++j, ++row_length)
364 if (colnums[j] != invalid_entry)
365 tmp_entries[row_length] = colnums[j];
366 else
367 break;
368 // now @p{rowstart} is the number of entries in this line
369
370 // Sort only beginning at the second entry, if optimized storage of
371 // diagonal entries is on.
372
373 // if this line is empty or has only one entry, don't sort
374 if (row_length > 1)
375 std::sort((store_diagonal_first_in_row) ? tmp_entries.begin() + 1 :
376 tmp_entries.begin(),
377 tmp_entries.begin() + row_length);
378
379 // insert column numbers into the new field
380 for (size_type j = 0; j < row_length; ++j)
381 new_colnums[next_free_entry++] = tmp_entries[j];
382
383 // note new start of this and the next row
384 rowstart[line] = next_row_start;
385 next_row_start = next_free_entry;
386
387 // some internal checks: either the matrix is not quadratic, or if it
388 // is, then the first element of this row must be the diagonal element
389 // (i.e. with column index==line number)
390 // this test only makes sense if we have written to the index
391 // rowstart_line in new_colnums which is the case if row_length is not 0,
392 // so check this first
394 (row_length != 0 && new_colnums[rowstart[line]] == line),
396 // assert that the first entry does not show up in the remaining ones
397 // and that the remaining ones are unique among themselves (this handles
398 // both cases, quadratic and rectangular matrices)
399 //
400 // the only exception here is if the row contains no entries at all
401 Assert((rowstart[line] == next_row_start) ||
402 (std::find(&new_colnums[rowstart[line] + 1],
403 &new_colnums[next_row_start],
404 new_colnums[rowstart[line]]) ==
405 &new_colnums[next_row_start]),
407 Assert((rowstart[line] == next_row_start) ||
408 (std::adjacent_find(&new_colnums[rowstart[line] + 1],
409 &new_colnums[next_row_start]) ==
410 &new_colnums[next_row_start]),
412 }
413
414 // assert that we have used all allocated space, no more and no less
415 Assert(next_free_entry == nonzero_elements, ExcInternalError());
416
417 // set iterator-past-the-end
418 rowstart[rows] = next_row_start;
419
420 // set colnums to the newly allocated array and delete previous content
421 // in the process
422 colnums = std::move(new_colnums);
423
424 // store the size
425 max_vec_len = nonzero_elements;
426
427 compressed = true;
428}
429
430
431
432void
434{
435 // first determine row lengths for each row. if the matrix is quadratic,
436 // then we might have to add an additional entry for the diagonal, if that
437 // is not yet present. as we have to call compress anyway later on, don't
438 // bother to check whether that diagonal entry is in a certain row or not
439 const bool do_diag_optimize = (sp.n_rows() == sp.n_cols());
440 std::vector<unsigned int> row_lengths(sp.n_rows());
441 for (size_type i = 0; i < sp.n_rows(); ++i)
442 {
443 row_lengths[i] = sp.row_length(i);
444 if (do_diag_optimize && !sp.exists(i, i))
445 ++row_lengths[i];
446 }
447 reinit(sp.n_rows(), sp.n_cols(), row_lengths);
448
449 // now enter all the elements into the matrix, if there are any. note that
450 // if the matrix is quadratic, then we already have the diagonal element
451 // preallocated
452 if (n_rows() != 0 && n_cols() != 0)
453 for (size_type row = 0; row < sp.n_rows(); ++row)
454 {
455 size_type *cols = &colnums[rowstart[row]] + (do_diag_optimize ? 1 : 0);
456 typename SparsityPattern::iterator col_num = sp.begin(row),
457 end_row = sp.end(row);
458
459 for (; col_num != end_row; ++col_num)
460 {
461 const size_type col = col_num->column();
462 if ((col != row) || !do_diag_optimize)
463 *cols++ = col;
464 }
465 }
466
467 // do not need to compress the sparsity pattern since we already have
468 // allocated the right amount of data, and the SparsityPattern data is
469 // sorted, too.
470 compressed = true;
471}
472
473
474
475// Use a special implementation for DynamicSparsityPattern where we can use
476// the column_number method to gain faster access to the
477// entries. DynamicSparsityPattern::iterator can show quadratic complexity in
478// case many rows are empty and the begin() method needs to jump to the next
479// free row. Otherwise, the code is exactly the same as above.
480void
482{
483 const bool do_diag_optimize = (dsp.n_rows() == dsp.n_cols());
484 const auto &row_index_set = dsp.row_index_set();
485
486 std::vector<unsigned int> row_lengths(dsp.n_rows());
487
488 if (row_index_set.size() == 0)
489 {
490 for (size_type i = 0; i < dsp.n_rows(); ++i)
491 {
492 row_lengths[i] = dsp.row_length(i);
493 if (do_diag_optimize && !dsp.exists(i, i))
494 ++row_lengths[i];
495 }
496 }
497 else
498 {
499 for (size_type i = 0; i < dsp.n_rows(); ++i)
500 {
501 if (row_index_set.is_element(i))
502 {
503 row_lengths[i] = dsp.row_length(i);
504 if (do_diag_optimize && !dsp.exists(i, i))
505 ++row_lengths[i];
506 }
507 else
508 {
509 // If the row i is not stored in the DynamicSparsityPattern we
510 // nevertheless need to allocate 1 entry per row for the
511 // "diagonal optimization". (We store a pointer to the next row
512 // in place of the repeated index i for the diagonal element.)
513 row_lengths[i] = do_diag_optimize ? 1 : 0;
514 }
515 }
516 }
517 reinit(dsp.n_rows(), dsp.n_cols(), row_lengths);
518
519 if (n_rows() != 0 && n_cols() != 0)
520 for (size_type row = 0; row < dsp.n_rows(); ++row)
521 {
522 size_type *cols = &colnums[rowstart[row]] + (do_diag_optimize ? 1 : 0);
523 const unsigned int row_length = dsp.row_length(row);
524 for (unsigned int index = 0; index < row_length; ++index)
525 {
526 const size_type col = dsp.column_number(row, index);
527 if ((col != row) || !do_diag_optimize)
528 *cols++ = col;
529 }
530 }
531
532 // do not need to compress the sparsity pattern since we already have
533 // allocated the right amount of data, and the SparsityPatternType data is
534 // sorted, too.
535 compressed = true;
536}
537
538
539
540template <typename number>
541void
543{
544 // first init with the number of entries per row. if this matrix is square
545 // then we also have to allocate memory for the diagonal entry, unless we
546 // have already counted it
547 const bool matrix_is_square = (matrix.m() == matrix.n());
548
549 std::vector<unsigned int> entries_per_row(matrix.m(), 0);
550 for (size_type row = 0; row < matrix.m(); ++row)
551 {
552 for (size_type col = 0; col < matrix.n(); ++col)
553 if (matrix(row, col) != 0)
554 ++entries_per_row[row];
555 if (matrix_is_square && (matrix(row, row) == 0))
556 ++entries_per_row[row];
557 }
558
559 reinit(matrix.m(), matrix.n(), entries_per_row);
560
561 // now set entries. if we enter entries row by row, then we'll get
562 // quadratic complexity in the number of entries per row. this is
563 // not usually a problem (we don't usually create dense matrices),
564 // but there are cases where it matters -- so we may as well be
565 // gentler and hand over a whole row of entries at a time
566 std::vector<size_type> column_indices;
567 column_indices.reserve(
568 entries_per_row.size() > 0 ?
569 *std::max_element(entries_per_row.begin(), entries_per_row.end()) :
570 0);
571 for (size_type row = 0; row < matrix.m(); ++row)
572 {
573 column_indices.resize(entries_per_row[row]);
574
575 size_type current_index = 0;
576 for (size_type col = 0; col < matrix.n(); ++col)
577 if (matrix(row, col) != 0)
578 {
579 column_indices[current_index] = col;
580 ++current_index;
581 }
582 else
583 // the (row,col) entry is zero; check if we need to add it
584 // anyway because it's the diagonal entry of a square
585 // matrix
586 if (matrix_is_square && (col == row))
587 {
588 column_indices[current_index] = row;
589 ++current_index;
590 }
591
592 // check that we really added the correct number of indices
593 Assert(current_index == entries_per_row[row], ExcInternalError());
594
595 // now bulk add all of these entries
596 add_entries(row, column_indices.begin(), column_indices.end(), true);
597 }
598
599 // finally compress
600 compress();
601}
602
603
604
605bool
607{
608 // let's try to be on the safe side of life by using multiple possibilities
609 // in the check for emptiness... (sorry for this kludge -- emptying matrices
610 // and freeing memory was not present in the original implementation and I
611 // don't know at how many places I missed something in adding it, so I try
612 // to be cautious. wb)
613 if ((rowstart == nullptr) || (n_rows() == 0) || (n_cols() == 0))
614 {
615 Assert(rowstart == nullptr, ExcInternalError());
616 Assert(n_rows() == 0, ExcInternalError());
617 Assert(n_cols() == 0, ExcInternalError());
618 Assert(colnums == nullptr, ExcInternalError());
620
621 return true;
622 }
623 return false;
624}
625
626
627
630{
631 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
635
636 // let's see whether there is something in this line
637 if (rowstart[i] == rowstart[i + 1])
638 return invalid_entry;
639
640 // If special storage of diagonals was requested, we can get the diagonal
641 // element faster by this query.
642 if (store_diagonal_first_in_row && (i == j))
643 return rowstart[i];
644
645 // all other entries are sorted, so we can use a binary search algorithm
646 //
647 // note that the entries are only sorted upon compression, so this would
648 // fail for non-compressed sparsity patterns; however, that is why the
649 // Assertion is at the top of this function, so it may not be called for
650 // noncompressed structures.
651 const size_type *sorted_region_start =
653 &colnums[rowstart[i]]);
654 const size_type *const p =
655 Utilities::lower_bound<const size_type *>(sorted_region_start,
656 &colnums[rowstart[i + 1]],
657 j);
658 if ((p != &colnums[rowstart[i + 1]]) && (*p == j))
659 return (p - colnums.get());
660 else
661 return invalid_entry;
662}
663
664
665
666bool
668{
669 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
672
673 for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k)
674 {
675 // entry already exists
676 if (colnums[k] == j)
677 return true;
678 }
679 return false;
680}
681
682
683
684std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
685SparsityPattern::matrix_position(const std::size_t global_index) const
686{
688 AssertIndexRange(global_index, n_nonzero_elements());
689
690 // first find the row in which the entry is located. for this note that the
691 // rowstart array indexes the global indices at which each row starts. since
692 // it is sorted, and since there is an element for the one-past-last row, we
693 // can simply use a bisection search on it
694 const size_type row =
695 (std::upper_bound(rowstart.get(), rowstart.get() + rows, global_index) -
696 rowstart.get() - 1);
697
698 // now, the column index is simple since that is what the colnums array
699 // stores:
700 const size_type col = colnums[global_index];
701
702 // so return the respective pair
703 return std::make_pair(row, col);
704}
705
706
707
710{
711 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
714
715 for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k)
716 {
717 // entry exists
718 if (colnums[k] == j)
719 return k - rowstart[i];
720 }
722}
723
724
725
728{
729 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
730 size_type b = 0;
731 for (size_type i = 0; i < n_rows(); ++i)
732 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
733 if (colnums[j] != invalid_entry)
734 {
735 if (static_cast<size_type>(
736 std::abs(static_cast<int>(i - colnums[j]))) > b)
737 b = std::abs(static_cast<signed int>(i - colnums[j]));
738 }
739 else
740 // leave if at the end of the entries of this line
741 break;
742 return b;
743}
744
745
746
749{
750 // if compress() has not yet been called, we can get the maximum number of
751 // elements per row using the stored value
752 if (!compressed)
753 return max_row_length;
754
755 // if compress() was called, we use a better algorithm which gives us a
756 // sharp bound
757 size_type m = 0;
758 for (size_type i = 1; i <= rows; ++i)
759 m = std::max(m, static_cast<size_type>(rowstart[i] - rowstart[i - 1]));
760
761 return m;
762}
763
764
765
766void
768{
769 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
773
774 for (std::size_t k = rowstart[i]; k < rowstart[i + 1]; ++k)
775 {
776 // entry already exists
777 if (colnums[k] == j)
778 return;
779 // empty entry found, put new entry here
780 if (colnums[k] == invalid_entry)
781 {
782 colnums[k] = j;
783 return;
784 }
785 }
786
787 // if we came thus far, something went wrong: there was not enough space in
788 // this line
789 Assert(false, ExcNotEnoughSpace(i, rowstart[i + 1] - rowstart[i]));
790}
791
792
793
794template <typename ForwardIterator>
795void
797 ForwardIterator begin,
798 ForwardIterator end,
799 const bool indices_are_sorted)
800{
801 AssertIndexRange(row, n_rows());
802 if (indices_are_sorted == true)
803 {
804 if (begin != end)
805 {
806 ForwardIterator it = begin;
807 bool has_larger_entries = false;
808 // skip diagonal
809 std::size_t k = rowstart[row] + store_diagonal_first_in_row;
810 for (; k < rowstart[row + 1]; ++k)
811 if (colnums[k] == invalid_entry)
812 break;
813 else if (colnums[k] >= *it)
814 {
815 has_larger_entries = true;
816 break;
817 }
818 if (has_larger_entries == false)
819 for (; it != end; ++it)
820 {
821 AssertIndexRange(*it, n_cols());
822 if (store_diagonal_first_in_row && *it == row)
823 continue;
824 Assert(k <= rowstart[row + 1],
826 rowstart[row + 1] - rowstart[row]));
827 colnums[k++] = *it;
828 }
829 else
830 // cannot just append the new range at the end, forward to the
831 // other function
832 for (ForwardIterator p = begin; p != end; ++p)
833 add(row, *p);
834 }
835 }
836 else
837 {
838 // forward to the other function.
839 for (ForwardIterator it = begin; it != end; ++it)
840 add(row, *it);
841 }
842}
843
844
845
846void
848 const ArrayView<const size_type> &columns,
849 const bool indices_are_sorted)
850{
851 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
852}
853
854
855
856void
858{
859 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
861 // Note that we only require a quadratic matrix here, no special treatment
862 // of diagonals
864
865 // loop over all elements presently in the sparsity pattern and add the
866 // transpose element. note:
867 //
868 // 1. that the sparsity pattern changes which we work on, but not the
869 // present row
870 //
871 // 2. that the @p{add} function can be called on elements that already exist
872 // without any harm
873 for (size_type row = 0; row < n_rows(); ++row)
874 for (size_type k = rowstart[row]; k < rowstart[row + 1]; ++k)
875 {
876 // check whether we are at the end of the entries of this row. if so,
877 // go to next row
878 if (colnums[k] == invalid_entry)
879 break;
880
881 // otherwise add the transpose entry if this is not the diagonal (that
882 // would not harm, only take time to check up)
883 if (colnums[k] != row)
884 add(colnums[k], row);
885 }
886}
887
888
889
890bool
892{
894 return false;
895
896 // it isn't quite necessary to compare *all* member variables. by only
897 // comparing the essential ones, we can say that two sparsity patterns are
898 // equal even if one is compressed and the other is not (in which case some
899 // of the member variables are not yet set correctly)
900 if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed)
901 return false;
902
903 if (rows > 0)
904 {
905 for (size_type i = 0; i < rows + 1; ++i)
906 if (rowstart[i] != sp2.rowstart[i])
907 return false;
908
909 for (size_type i = 0; i < rowstart[rows]; ++i)
910 if (colnums[i] != sp2.colnums[i])
911 return false;
912 }
913
914 return true;
915}
916
917
918
919void
920SparsityPattern::print(std::ostream &out) const
921{
922 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
923
924 AssertThrow(out.fail() == false, ExcIO());
925
926 for (size_type i = 0; i < n_rows(); ++i)
927 {
928 out << '[' << i;
929 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
930 if (colnums[j] != invalid_entry)
931 out << ',' << colnums[j];
932 out << ']' << std::endl;
933 }
934
935 AssertThrow(out.fail() == false, ExcIO());
936}
937
938
939
940void
941SparsityPattern::print_gnuplot(std::ostream &out) const
942{
943 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
944
945 AssertThrow(out.fail() == false, ExcIO());
946
947 for (size_type i = 0; i < n_rows(); ++i)
948 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
949 if (colnums[j] != invalid_entry)
950 // while matrix entries are usually written (i,j), with i vertical and
951 // j horizontal, gnuplot output is x-y, that is we have to exchange
952 // the order of output
953 out << colnums[j] << " " << -static_cast<signed int>(i) << std::endl;
954
955 AssertThrow(out.fail() == false, ExcIO());
956}
957
958
959
960void
961SparsityPattern::print_svg(std::ostream &out) const
962{
963 const unsigned int m = this->n_rows();
964 const unsigned int n = this->n_cols();
965 out
966 << "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\" viewBox=\"0 0 "
967 << n + 2 << " " << m + 2
968 << " \">\n"
969 "<style type=\"text/css\" >\n"
970 " <![CDATA[\n"
971 " rect.pixel {\n"
972 " fill: #ff0000;\n"
973 " }\n"
974 " ]]>\n"
975 " </style>\n\n"
976 " <rect width=\""
977 << n + 2 << "\" height=\"" << m + 2
978 << "\" fill=\"rgb(128, 128, 128)\"/>\n"
979 " <rect x=\"1\" y=\"1\" width=\""
980 << n + 0.1 << "\" height=\"" << m + 0.1
981 << "\" fill=\"rgb(255, 255, 255)\"/>\n\n";
982
983 for (const auto &entry : *this)
984 {
985 out << " <rect class=\"pixel\" x=\"" << entry.column() + 1 << "\" y=\""
986 << entry.row() + 1 << "\" width=\".9\" height=\".9\"/>\n";
987 }
988 out << "</svg>" << std::endl;
989}
990
991
992
993void
994SparsityPattern::block_write(std::ostream &out) const
995{
996 AssertThrow(out.fail() == false, ExcIO());
997
998 // first the simple objects, bracketed in [...]
999 out << '[' << max_dim << ' ' << n_rows() << ' ' << n_cols() << ' '
1000 << max_vec_len << ' ' << max_row_length << ' ' << compressed << ' '
1001 << store_diagonal_first_in_row << "][";
1002 // then write out real data
1003 out.write(reinterpret_cast<const char *>(rowstart.get()),
1004 reinterpret_cast<const char *>(rowstart.get() + max_dim + 1) -
1005 reinterpret_cast<const char *>(rowstart.get()));
1006 out << "][";
1007 out.write(reinterpret_cast<const char *>(colnums.get()),
1008 reinterpret_cast<const char *>(colnums.get() + max_vec_len) -
1009 reinterpret_cast<const char *>(colnums.get()));
1010 out << ']';
1011
1012 AssertThrow(out.fail() == false, ExcIO());
1013}
1014
1015
1016
1017void
1019{
1020 AssertThrow(in.fail() == false, ExcIO());
1021
1022 char c;
1023
1024 // first read in simple data
1025 in >> c;
1026 AssertThrow(c == '[', ExcIO());
1027 in >> max_dim >> rows >> cols >> max_vec_len >> max_row_length >>
1029
1030 in >> c;
1031 AssertThrow(c == ']', ExcIO());
1032 in >> c;
1033 AssertThrow(c == '[', ExcIO());
1034
1035 // reallocate space
1036 rowstart = std::make_unique<std::size_t[]>(max_dim + 1);
1037 colnums = std::make_unique<size_type[]>(max_vec_len);
1038
1039 // then read data
1040 in.read(reinterpret_cast<char *>(rowstart.get()),
1041 reinterpret_cast<char *>(rowstart.get() + max_dim + 1) -
1042 reinterpret_cast<char *>(rowstart.get()));
1043 in >> c;
1044 AssertThrow(c == ']', ExcIO());
1045 in >> c;
1046 AssertThrow(c == '[', ExcIO());
1047 in.read(reinterpret_cast<char *>(colnums.get()),
1048 reinterpret_cast<char *>(colnums.get() + max_vec_len) -
1049 reinterpret_cast<char *>(colnums.get()));
1050 in >> c;
1051 AssertThrow(c == ']', ExcIO());
1052}
1053
1054
1055
1056std::size_t
1058{
1059 return (max_dim * sizeof(size_type) + sizeof(*this) +
1060 max_vec_len * sizeof(size_type));
1061}
1062
1063
1064
1065#ifndef DOXYGEN
1066// explicit instantiations
1067template void
1068SparsityPattern::copy_from<float>(const FullMatrix<float> &);
1069template void
1070SparsityPattern::copy_from<double>(const FullMatrix<double> &);
1071
1072template void
1073SparsityPattern::add_entries<const SparsityPattern::size_type *>(
1074 const size_type,
1075 const size_type *,
1076 const size_type *,
1077 const bool);
1078# ifndef DEAL_II_VECTOR_ITERATOR_IS_POINTER
1079template void
1081 std::vector<SparsityPattern::size_type>::const_iterator>(
1082 const size_type,
1083 std::vector<size_type>::const_iterator,
1084 std::vector<size_type>::const_iterator,
1085 const bool);
1086# endif
1087template void
1088SparsityPattern::add_entries<std::vector<SparsityPattern::size_type>::iterator>(
1089 const size_type,
1090 std::vector<size_type>::iterator,
1091 std::vector<size_type>::iterator,
1092 const bool);
1093#endif
1094
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
iterator begin() const
Definition array_view.h:707
bool empty() const
Definition array_view.h:698
iterator end() const
Definition array_view.h:716
std::size_t size() const
Definition array_view.h:689
const IndexSet & row_index_set() const
size_type row_length(const size_type row) const
size_type column_number(const size_type row, const size_type index) const
bool exists(const size_type i, const size_type j) const
virtual void resize(const size_type rows, const size_type cols)
size_type n_rows() const
size_type n_cols() const
std::pair< size_type, size_type > matrix_position(const std::size_t global_index) const
void block_write(std::ostream &out) const
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
void print_svg(std::ostream &out) const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
size_type bandwidth() const
void print_gnuplot(std::ostream &out) const
virtual void add_row_entries(const size_type &row, const ArrayView< const size_type > &columns, const bool indices_are_sorted=false) override
bool is_compressed() const
SparsityPattern & operator=(const SparsityPattern &)
std::size_t n_nonzero_elements() const
bool exists(const size_type i, const size_type j) const
iterator begin() const
std::unique_ptr< size_type[]> colnums
std::size_t max_vec_len
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
static constexpr size_type invalid_entry
size_type row_position(const size_type i, const size_type j) const
std::unique_ptr< std::size_t[]> rowstart
void print(std::ostream &out) const
void add(const size_type i, const size_type j)
size_type max_entries_per_row() const
unsigned int max_row_length
iterator end() const
size_type operator()(const size_type i, const size_type j) const
unsigned int row_length(const size_type row) const
std::size_t memory_consumption() const
bool operator==(const SparsityPattern &) const
void block_read(std::istream &in)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
static ::ExceptionBase & ExcNotEnoughSpace(int arg1, int arg2)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcMatrixIsCompressed()
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition utilities.h:1033
constexpr types::global_dof_index invalid_size_type
Definition types.h:250
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)