Reference documentation for deal.II version 9.6.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
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 size_type 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 (size_type 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
630bool
632{
633 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
636
637 for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k)
638 {
639 // entry already exists
640 if (colnums[k] == j)
641 return true;
642 }
643 return false;
644}
645
646
647
648std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
649SparsityPattern::matrix_position(const std::size_t global_index) const
650{
652 AssertIndexRange(global_index, n_nonzero_elements());
653
654 // first find the row in which the entry is located. for this note that the
655 // rowstart array indexes the global indices at which each row starts. since
656 // it is sorted, and since there is an element for the one-past-last row, we
657 // can simply use a bisection search on it
658 const size_type row =
659 (std::upper_bound(rowstart.get(), rowstart.get() + rows, global_index) -
660 rowstart.get() - 1);
661
662 // now, the column index is simple since that is what the colnums array
663 // stores:
664 const size_type col = colnums[global_index];
665
666 // so return the respective pair
667 return std::make_pair(row, col);
668}
669
670
671
674{
675 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
676 size_type b = 0;
677 for (size_type i = 0; i < n_rows(); ++i)
678 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
679 if (colnums[j] != invalid_entry)
680 {
681 if (static_cast<size_type>(
682 std::abs(static_cast<int>(i - colnums[j]))) > b)
683 b = std::abs(static_cast<signed int>(i - colnums[j]));
684 }
685 else
686 // leave if at the end of the entries of this line
687 break;
688 return b;
689}
690
691
692
695{
696 // if compress() has not yet been called, we can get the maximum number of
697 // elements per row using the stored value
698 if (!compressed)
699 return max_row_length;
700
701 // if compress() was called, we use a better algorithm which gives us a
702 // sharp bound
703 size_type m = 0;
704 for (size_type i = 1; i <= rows; ++i)
705 m = std::max(m, static_cast<size_type>(rowstart[i] - rowstart[i - 1]));
706
707 return m;
708}
709
710
711
714{
715 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
719
720 // let's see whether there is something in this line
721 if (rowstart[i] == rowstart[i + 1])
722 return invalid_entry;
723
724 // If special storage of diagonals was requested, we can get the diagonal
725 // element faster by this query.
726 if (store_diagonal_first_in_row && (i == j))
727 return rowstart[i];
728
729 // all other entries are sorted, so we can use a binary search algorithm
730 //
731 // note that the entries are only sorted upon compression, so this would
732 // fail for non-compressed sparsity patterns; however, that is why the
733 // Assertion is at the top of this function, so it may not be called for
734 // noncompressed structures.
735 const size_type *sorted_region_start =
737 &colnums[rowstart[i]]);
738 const size_type *const p =
740 &colnums[rowstart[i + 1]],
741 j);
742 if ((p != &colnums[rowstart[i + 1]]) && (*p == j))
743 return (p - colnums.get());
744 else
745 return invalid_entry;
746}
747
748
749
750void
752{
753 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
757
758 for (std::size_t k = rowstart[i]; k < rowstart[i + 1]; ++k)
759 {
760 // entry already exists
761 if (colnums[k] == j)
762 return;
763 // empty entry found, put new entry here
764 if (colnums[k] == invalid_entry)
765 {
766 colnums[k] = j;
767 return;
768 }
769 }
770
771 // if we came thus far, something went wrong: there was not enough space in
772 // this line
773 Assert(false, ExcNotEnoughSpace(i, rowstart[i + 1] - rowstart[i]));
774}
775
776
777
778template <typename ForwardIterator>
779void
781 ForwardIterator begin,
782 ForwardIterator end,
783 const bool indices_are_sorted)
784{
785 AssertIndexRange(row, n_rows());
786 if (indices_are_sorted == true)
787 {
788 if (begin != end)
789 {
790 ForwardIterator it = begin;
791 bool has_larger_entries = false;
792 // skip diagonal
793 std::size_t k = rowstart[row] + store_diagonal_first_in_row;
794 for (; k < rowstart[row + 1]; ++k)
795 if (colnums[k] == invalid_entry)
796 break;
797 else if (colnums[k] >= *it)
798 {
799 has_larger_entries = true;
800 break;
801 }
802 if (has_larger_entries == false)
803 for (; it != end; ++it)
804 {
805 AssertIndexRange(*it, n_cols());
806 if (store_diagonal_first_in_row && *it == row)
807 continue;
808 Assert(k <= rowstart[row + 1],
810 rowstart[row + 1] - rowstart[row]));
811 colnums[k++] = *it;
812 }
813 else
814 // cannot just append the new range at the end, forward to the
815 // other function
816 for (ForwardIterator p = begin; p != end; ++p)
817 add(row, *p);
818 }
819 }
820 else
821 {
822 // forward to the other function.
823 for (ForwardIterator it = begin; it != end; ++it)
824 add(row, *it);
825 }
826}
827
828
829
830void
832 const ArrayView<const size_type> &columns,
833 const bool indices_are_sorted)
834{
835 add_entries(row, columns.begin(), columns.end(), indices_are_sorted);
836}
837
838
839
840void
842{
843 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
845 // Note that we only require a quadratic matrix here, no special treatment
846 // of diagonals
848
849 // loop over all elements presently in the sparsity pattern and add the
850 // transpose element. note:
851 //
852 // 1. that the sparsity pattern changes which we work on, but not the
853 // present row
854 //
855 // 2. that the @p{add} function can be called on elements that already exist
856 // without any harm
857 for (size_type row = 0; row < n_rows(); ++row)
858 for (size_type k = rowstart[row]; k < rowstart[row + 1]; ++k)
859 {
860 // check whether we are at the end of the entries of this row. if so,
861 // go to next row
862 if (colnums[k] == invalid_entry)
863 break;
864
865 // otherwise add the transpose entry if this is not the diagonal (that
866 // would not harm, only take time to check up)
867 if (colnums[k] != row)
868 add(colnums[k], row);
869 }
870}
871
872
873
876{
877 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
880
881 for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k)
882 {
883 // entry exists
884 if (colnums[k] == j)
885 return k - rowstart[i];
886 }
888}
889
890
891
892void
893SparsityPattern::print(std::ostream &out) const
894{
895 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
896
897 AssertThrow(out.fail() == false, ExcIO());
898
899 for (size_type i = 0; i < n_rows(); ++i)
900 {
901 out << '[' << i;
902 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
903 if (colnums[j] != invalid_entry)
904 out << ',' << colnums[j];
905 out << ']' << std::endl;
906 }
907
908 AssertThrow(out.fail() == false, ExcIO());
909}
910
911
912
913void
914SparsityPattern::print_gnuplot(std::ostream &out) const
915{
916 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
917
918 AssertThrow(out.fail() == false, ExcIO());
919
920 for (size_type i = 0; i < n_rows(); ++i)
921 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
922 if (colnums[j] != invalid_entry)
923 // while matrix entries are usually written (i,j), with i vertical and
924 // j horizontal, gnuplot output is x-y, that is we have to exchange
925 // the order of output
926 out << colnums[j] << " " << -static_cast<signed int>(i) << std::endl;
927
928 AssertThrow(out.fail() == false, ExcIO());
929}
930
931
932
933void
934SparsityPattern::print_svg(std::ostream &out) const
935{
936 const unsigned int m = this->n_rows();
937 const unsigned int n = this->n_cols();
938 out
939 << "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\" viewBox=\"0 0 "
940 << n + 2 << " " << m + 2
941 << " \">\n"
942 "<style type=\"text/css\" >\n"
943 " <![CDATA[\n"
944 " rect.pixel {\n"
945 " fill: #ff0000;\n"
946 " }\n"
947 " ]]>\n"
948 " </style>\n\n"
949 " <rect width=\""
950 << n + 2 << "\" height=\"" << m + 2
951 << "\" fill=\"rgb(128, 128, 128)\"/>\n"
952 " <rect x=\"1\" y=\"1\" width=\""
953 << n + 0.1 << "\" height=\"" << m + 0.1
954 << "\" fill=\"rgb(255, 255, 255)\"/>\n\n";
955
956 for (const auto &entry : *this)
957 {
958 out << " <rect class=\"pixel\" x=\"" << entry.column() + 1 << "\" y=\""
959 << entry.row() + 1 << "\" width=\".9\" height=\".9\"/>\n";
960 }
961 out << "</svg>" << std::endl;
962}
963
964
965
966void
967SparsityPattern::block_write(std::ostream &out) const
968{
969 AssertThrow(out.fail() == false, ExcIO());
970
971 // first the simple objects, bracketed in [...]
972 out << '[' << max_dim << ' ' << n_rows() << ' ' << n_cols() << ' '
973 << max_vec_len << ' ' << max_row_length << ' ' << compressed << ' '
975 // then write out real data
976 out.write(reinterpret_cast<const char *>(rowstart.get()),
977 reinterpret_cast<const char *>(rowstart.get() + max_dim + 1) -
978 reinterpret_cast<const char *>(rowstart.get()));
979 out << "][";
980 out.write(reinterpret_cast<const char *>(colnums.get()),
981 reinterpret_cast<const char *>(colnums.get() + max_vec_len) -
982 reinterpret_cast<const char *>(colnums.get()));
983 out << ']';
984
985 AssertThrow(out.fail() == false, ExcIO());
986}
987
988
989
990void
992{
993 AssertThrow(in.fail() == false, ExcIO());
994
995 char c;
996
997 // first read in simple data
998 in >> c;
999 AssertThrow(c == '[', ExcIO());
1000 in >> max_dim >> rows >> cols >> max_vec_len >> max_row_length >>
1002
1003 in >> c;
1004 AssertThrow(c == ']', ExcIO());
1005 in >> c;
1006 AssertThrow(c == '[', ExcIO());
1007
1008 // reallocate space
1009 rowstart = std::make_unique<std::size_t[]>(max_dim + 1);
1010 colnums = std::make_unique<size_type[]>(max_vec_len);
1011
1012 // then read data
1013 in.read(reinterpret_cast<char *>(rowstart.get()),
1014 reinterpret_cast<char *>(rowstart.get() + max_dim + 1) -
1015 reinterpret_cast<char *>(rowstart.get()));
1016 in >> c;
1017 AssertThrow(c == ']', ExcIO());
1018 in >> c;
1019 AssertThrow(c == '[', ExcIO());
1020 in.read(reinterpret_cast<char *>(colnums.get()),
1021 reinterpret_cast<char *>(colnums.get() + max_vec_len) -
1022 reinterpret_cast<char *>(colnums.get()));
1023 in >> c;
1024 AssertThrow(c == ']', ExcIO());
1025}
1026
1027
1028
1029std::size_t
1031{
1032 return (max_dim * sizeof(size_type) + sizeof(*this) +
1033 max_vec_len * sizeof(size_type));
1034}
1035
1036
1037
1038#ifndef DOXYGEN
1039// explicit instantiations
1040template void
1042template void
1044
1045template void
1047 const size_type,
1048 const size_type *,
1049 const size_type *,
1050 const bool);
1051# ifndef DEAL_II_VECTOR_ITERATOR_IS_POINTER
1052template void
1054 std::vector<SparsityPattern::size_type>::const_iterator>(
1055 const size_type,
1056 std::vector<size_type>::const_iterator,
1057 std::vector<size_type>::const_iterator,
1058 const bool);
1059# endif
1060template void
1061SparsityPattern::add_entries<std::vector<SparsityPattern::size_type>::iterator>(
1062 const size_type,
1063 std::vector<size_type>::iterator,
1064 std::vector<size_type>::iterator,
1065 const bool);
1066#endif
1067
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:949
iterator begin() const
Definition array_view.h:702
bool empty() const
Definition array_view.h:693
iterator end() const
Definition array_view.h:711
std::size_t size() const
Definition array_view.h:684
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
void block_read(std::istream &in)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
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:1032
const types::global_dof_index invalid_size_type
Definition types.h:233
::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 > &)