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