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