Reference documentation for deal.II version 9.4.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\}}\)
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
34#ifdef DEAL_II_MSVC
35__declspec(selectany) // Weak extern binding due to multiple link error
36#endif
38
39
41 : max_dim(0)
42 , rows(0)
43 , cols(0)
44 , max_vec_len(0)
45 , max_row_length(0)
46 , rowstart(nullptr)
47 , colnums(nullptr)
48 , compressed(false)
49{}
50
51
52
55 , store_diagonal_first_in_row(false)
56{
57 reinit(0, 0, 0);
58}
59
60
61
64 , store_diagonal_first_in_row(false)
65{
66 (void)s;
67 Assert(s.empty(),
69 "This constructor can only be called if the provided argument "
70 "is the sparsity pattern for an empty matrix. This constructor can "
71 "not be used to copy-construct a non-empty sparsity pattern."));
72
73 reinit(0, 0, 0);
74}
75
76
77
79 const size_type n,
80 const unsigned int max_per_row)
82 , store_diagonal_first_in_row(m == n)
83{
84 reinit(m, n, max_per_row);
85}
86
87
88
90 const size_type n,
91 const std::vector<unsigned int> &row_lengths)
93 , store_diagonal_first_in_row(m == n)
94{
95 reinit(m, n, row_lengths);
96}
97
98
99
101 const unsigned int max_per_row)
103{
104 reinit(m, m, max_per_row);
105}
106
107
108
110 const std::vector<unsigned int> &row_lengths)
112{
113 reinit(m, m, row_lengths);
114}
115
116
117
119 const unsigned int max_per_row,
120 const size_type extra_off_diagonals)
122{
123 Assert(original.rows == original.cols, ExcNotQuadratic());
125
126 reinit(original.rows, original.cols, max_per_row);
127
128 // now copy the entries from the other object
129 for (size_type row = 0; row < original.rows; ++row)
130 {
131 // copy the elements of this row of the other object
132 //
133 // note that the first object actually is the main-diagonal element,
134 // which we need not copy
135 //
136 // we do the copying in two steps: first we note that the elements in
137 // @p{original} are sorted, so we may first copy all the elements up to
138 // the first side-diagonal one which is to be filled in. then we insert
139 // the side-diagonals, finally copy the rest from that element onwards
140 // which is not a side-diagonal any more.
141 const size_type *const original_row_start =
142 &original.colnums[original.rowstart[row]] + 1;
143 // the following requires that @p{original} be compressed since
144 // otherwise there might be invalid_entry's
145 const size_type *const original_row_end =
146 &original.colnums[original.rowstart[row + 1]];
147
148 // find pointers before and after extra off-diagonals. if at top or
149 // bottom of matrix, then set these pointers such that no copying is
150 // necessary (see the @p{copy} commands)
151 const size_type *const original_last_before_side_diagonals =
152 (row > extra_off_diagonals ?
153 Utilities::lower_bound(original_row_start,
154 original_row_end,
155 row - extra_off_diagonals) :
156 original_row_start);
157
158 const size_type *const original_first_after_side_diagonals =
159 (row < rows - extra_off_diagonals - 1 ?
160 std::upper_bound(original_row_start,
161 original_row_end,
162 row + extra_off_diagonals) :
163 original_row_end);
164
165 // find first free slot. the first slot in each row is the diagonal
166 // element
167 size_type *next_free_slot = &colnums[rowstart[row]] + 1;
168
169 // copy elements before side-diagonals
170 next_free_slot = std::copy(original_row_start,
171 original_last_before_side_diagonals,
172 next_free_slot);
173
174 // insert left and right side-diagonals
175 for (size_type i = 1; i <= std::min(row, extra_off_diagonals);
176 ++i, ++next_free_slot)
177 *next_free_slot = row - i;
178 for (size_type i = 1; i <= std::min(extra_off_diagonals, rows - row - 1);
179 ++i, ++next_free_slot)
180 *next_free_slot = row + i;
181
182 // copy rest
183 next_free_slot = std::copy(original_first_after_side_diagonals,
184 original_row_end,
185 next_free_slot);
186
187 // this error may happen if the sum of previous elements per row and
188 // those of the new diagonals exceeds the maximum number of elements per
189 // row given to this constructor
190 Assert(next_free_slot <= &colnums[rowstart[row + 1]],
191 ExcNotEnoughSpace(0, rowstart[row + 1] - rowstart[row]));
192 }
193}
194
195
196
199{
200 (void)s;
201 Assert(s.empty(),
203 "This operator can only be called if the provided argument "
204 "is the sparsity pattern for an empty matrix. This operator can "
205 "not be used to copy a non-empty sparsity pattern."));
206
207 Assert(this->empty(),
208 ExcMessage("This operator can only be called if the current object is "
209 "empty."));
210
211 return *this;
212}
213
214
215
216void
218 const size_type n,
219 const unsigned int max_per_row)
220{
221 // simply map this function to the other @p{reinit} function
222 const std::vector<unsigned int> row_lengths(m, max_per_row);
223 reinit(m, n, row_lengths);
224}
225
226
227
228void
230 const size_type n,
231 const ArrayView<const unsigned int> &row_lengths)
232{
233 AssertDimension(row_lengths.size(), m);
234
235 rows = m;
236 cols = n;
237
238 // delete empty matrices
239 if ((m == 0) || (n == 0))
240 {
241 rowstart.reset();
242 colnums.reset();
243
244 max_vec_len = max_dim = rows = cols = 0;
245 // if dimension is zero: ignore max_per_row
246 max_row_length = 0;
247 compressed = false;
248
249 return;
250 }
251
252 // first, if the matrix is quadratic, we will have to make sure that each
253 // row has at least one entry for the diagonal element. make this more
254 // obvious by having a variable which we can query
256
257 // find out how many entries we need in the @p{colnums} array. if this
258 // number is larger than @p{max_vec_len}, then we will need to reallocate
259 // memory
260 //
261 // note that the number of elements per row is bounded by the number of
262 // columns
263 //
264 std::size_t vec_len = 0;
265 for (size_type i = 0; i < m; ++i)
266 vec_len += std::min(static_cast<size_type>(store_diagonal_first_in_row ?
267 std::max(row_lengths[i], 1U) :
268 row_lengths[i]),
269 n);
270
271 // sometimes, no entries are requested in the matrix (this most often
272 // happens when blocks in a block matrix are simply zero). in that case,
273 // allocate exactly one element, to have a valid pointer to some memory
274 if (vec_len == 0)
275 {
276 vec_len = 1;
277 max_vec_len = vec_len;
278 colnums = std::make_unique<size_type[]>(max_vec_len);
279 }
280
282 (row_lengths.size() == 0 ?
283 0 :
284 std::min(static_cast<size_type>(
285 *std::max_element(row_lengths.begin(), row_lengths.end())),
286 n));
287
288 if (store_diagonal_first_in_row && (max_row_length == 0) && (m != 0))
289 max_row_length = 1;
290
291 // allocate memory for the rowstart values, if necessary. even though we
292 // re-set the pointers again immediately after deleting their old content,
293 // set them to zero in between because the allocation might fail, in which
294 // case we get an exception and the destructor of this object will be called
295 // -- where we look at the non-nullness of the (now invalid) pointer again
296 // and try to delete the memory a second time.
297 if (rows > max_dim)
298 {
299 max_dim = rows;
300 rowstart = std::make_unique<std::size_t[]>(max_dim + 1);
301 }
302
303 // allocate memory for the column numbers if necessary
304 if (vec_len > max_vec_len)
305 {
306 max_vec_len = vec_len;
307 colnums = std::make_unique<size_type[]>(max_vec_len);
308 }
309
310 // set the rowstart array
311 rowstart[0] = 0;
312 for (size_type i = 1; i <= rows; ++i)
313 rowstart[i] =
314 rowstart[i - 1] +
316 std::max(std::min(static_cast<size_type>(row_lengths[i - 1]), n),
317 static_cast<size_type>(1U)) :
318 std::min(static_cast<size_type>(row_lengths[i - 1]), n));
319 Assert((rowstart[rows] == vec_len) ||
320 ((vec_len == 1) && (rowstart[rows] == 0)),
322
323 // preset the column numbers by a value indicating it is not in use
324 std::fill_n(colnums.get(), vec_len, invalid_entry);
325
326 // if diagonal elements are special: let the first entry in each row be the
327 // diagonal value
329 for (size_type i = 0; i < rows; ++i)
330 colnums[rowstart[i]] = i;
331
332 compressed = false;
333}
334
335
336
337void
339{
340 // nothing to do if the object corresponds to an empty matrix
341 if ((rowstart == nullptr) && (colnums == nullptr))
342 {
343 compressed = true;
344 return;
345 }
346
347 // do nothing if already compressed
348 if (compressed)
349 return;
350
351 size_type next_free_entry = 0, next_row_start = 0, row_length = 0;
352
353 // first find out how many non-zero elements there are, in order to allocate
354 // the right amount of memory
355 const std::size_t nonzero_elements =
356 std::count_if(&colnums[rowstart[0]],
358 [](const size_type col) { return col != invalid_entry; });
359 // now allocate the respective memory
360 std::unique_ptr<size_type[]> new_colnums(new size_type[nonzero_elements]);
361
362
363 // reserve temporary storage to store the entries of one row
364 std::vector<size_type> tmp_entries(max_row_length);
365
366 // Traverse all rows
367 for (size_type line = 0; line < rows; ++line)
368 {
369 // copy used entries, break if first unused entry is reached
370 row_length = 0;
371 for (size_type j = rowstart[line]; j < rowstart[line + 1];
372 ++j, ++row_length)
373 if (colnums[j] != invalid_entry)
374 tmp_entries[row_length] = colnums[j];
375 else
376 break;
377 // now @p{rowstart} is the number of entries in this line
378
379 // Sort only beginning at the second entry, if optimized storage of
380 // diagonal entries is on.
381
382 // if this line is empty or has only one entry, don't sort
383 if (row_length > 1)
384 std::sort((store_diagonal_first_in_row) ? tmp_entries.begin() + 1 :
385 tmp_entries.begin(),
386 tmp_entries.begin() + row_length);
387
388 // insert column numbers into the new field
389 for (size_type j = 0; j < row_length; ++j)
390 new_colnums[next_free_entry++] = tmp_entries[j];
391
392 // note new start of this and the next row
393 rowstart[line] = next_row_start;
394 next_row_start = next_free_entry;
395
396 // some internal checks: either the matrix is not quadratic, or if it
397 // is, then the first element of this row must be the diagonal element
398 // (i.e. with column index==line number)
399 // this test only makes sense if we have written to the index
400 // rowstart_line in new_colnums which is the case if row_length is not 0,
401 // so check this first
403 (row_length != 0 && new_colnums[rowstart[line]] == line),
405 // assert that the first entry does not show up in the remaining ones
406 // and that the remaining ones are unique among themselves (this handles
407 // both cases, quadratic and rectangular matrices)
408 //
409 // the only exception here is if the row contains no entries at all
410 Assert((rowstart[line] == next_row_start) ||
411 (std::find(&new_colnums[rowstart[line] + 1],
412 &new_colnums[next_row_start],
413 new_colnums[rowstart[line]]) ==
414 &new_colnums[next_row_start]),
416 Assert((rowstart[line] == next_row_start) ||
417 (std::adjacent_find(&new_colnums[rowstart[line] + 1],
418 &new_colnums[next_row_start]) ==
419 &new_colnums[next_row_start]),
421 }
422
423 // assert that we have used all allocated space, no more and no less
424 Assert(next_free_entry == nonzero_elements, ExcInternalError());
425
426 // set iterator-past-the-end
427 rowstart[rows] = next_row_start;
428
429 // set colnums to the newly allocated array and delete previous content
430 // in the process
431 colnums = std::move(new_colnums);
432
433 // store the size
434 max_vec_len = nonzero_elements;
435
436 compressed = true;
437}
438
439
440
441void
443{
444 // first determine row lengths for each row. if the matrix is quadratic,
445 // then we might have to add an additional entry for the diagonal, if that
446 // is not yet present. as we have to call compress anyway later on, don't
447 // bother to check whether that diagonal entry is in a certain row or not
448 const bool do_diag_optimize = (sp.n_rows() == sp.n_cols());
449 std::vector<unsigned int> row_lengths(sp.n_rows());
450 for (size_type i = 0; i < sp.n_rows(); ++i)
451 {
452 row_lengths[i] = sp.row_length(i);
453 if (do_diag_optimize && !sp.exists(i, i))
454 ++row_lengths[i];
455 }
456 reinit(sp.n_rows(), sp.n_cols(), row_lengths);
457
458 // now enter all the elements into the matrix, if there are any. note that
459 // if the matrix is quadratic, then we already have the diagonal element
460 // preallocated
461 if (n_rows() != 0 && n_cols() != 0)
462 for (size_type row = 0; row < sp.n_rows(); ++row)
463 {
464 size_type *cols = &colnums[rowstart[row]] + (do_diag_optimize ? 1 : 0);
465 typename SparsityPattern::iterator col_num = sp.begin(row),
466 end_row = sp.end(row);
467
468 for (; col_num != end_row; ++col_num)
469 {
470 const size_type col = col_num->column();
471 if ((col != row) || !do_diag_optimize)
472 *cols++ = col;
473 }
474 }
475
476 // do not need to compress the sparsity pattern since we already have
477 // allocated the right amount of data, and the SparsityPattern data is
478 // sorted, too.
479 compressed = true;
480}
481
482
483
484// Use a special implementation for DynamicSparsityPattern where we can use
485// the column_number method to gain faster access to the
486// entries. DynamicSparsityPattern::iterator can show quadratic complexity in
487// case many rows are empty and the begin() method needs to jump to the next
488// free row. Otherwise, the code is exactly the same as above.
489void
491{
492 const bool do_diag_optimize = (dsp.n_rows() == dsp.n_cols());
493 const auto &row_index_set = dsp.row_index_set();
494
495 std::vector<unsigned int> row_lengths(dsp.n_rows());
496
497 if (row_index_set.size() == 0)
498 {
499 for (size_type i = 0; i < dsp.n_rows(); ++i)
500 {
501 row_lengths[i] = dsp.row_length(i);
502 if (do_diag_optimize && !dsp.exists(i, i))
503 ++row_lengths[i];
504 }
505 }
506 else
507 {
508 for (size_type i = 0; i < dsp.n_rows(); ++i)
509 {
510 if (row_index_set.is_element(i))
511 {
512 row_lengths[i] = dsp.row_length(i);
513 if (do_diag_optimize && !dsp.exists(i, i))
514 ++row_lengths[i];
515 }
516 else
517 {
518 // If the row i is not stored in the DynamicSparsityPattern we
519 // nevertheless need to allocate 1 entry per row for the
520 // "diagonal optimization". (We store a pointer to the next row
521 // in place of the repeated index i for the diagonal element.)
522 row_lengths[i] = do_diag_optimize ? 1 : 0;
523 }
524 }
525 }
526 reinit(dsp.n_rows(), dsp.n_cols(), row_lengths);
527
528 if (n_rows() != 0 && n_cols() != 0)
529 for (size_type row = 0; row < dsp.n_rows(); ++row)
530 {
531 size_type *cols = &colnums[rowstart[row]] + (do_diag_optimize ? 1 : 0);
532 const unsigned int row_length = dsp.row_length(row);
533 for (unsigned int index = 0; index < row_length; ++index)
534 {
535 const size_type col = dsp.column_number(row, index);
536 if ((col != row) || !do_diag_optimize)
537 *cols++ = col;
538 }
539 }
540
541 // do not need to compress the sparsity pattern since we already have
542 // allocated the right amount of data, and the SparsityPatternType data is
543 // sorted, too.
544 compressed = true;
545}
546
547
548
549template <typename number>
550void
552{
553 // first init with the number of entries per row. if this matrix is square
554 // then we also have to allocate memory for the diagonal entry, unless we
555 // have already counted it
556 const bool matrix_is_square = (matrix.m() == matrix.n());
557
558 std::vector<unsigned int> entries_per_row(matrix.m(), 0);
559 for (size_type row = 0; row < matrix.m(); ++row)
560 {
561 for (size_type col = 0; col < matrix.n(); ++col)
562 if (matrix(row, col) != 0)
563 ++entries_per_row[row];
564 if (matrix_is_square && (matrix(row, row) == 0))
565 ++entries_per_row[row];
566 }
567
568 reinit(matrix.m(), matrix.n(), entries_per_row);
569
570 // now set entries. if we enter entries row by row, then we'll get
571 // quadratic complexity in the number of entries per row. this is
572 // not usually a problem (we don't usually create dense matrices),
573 // but there are cases where it matters -- so we may as well be
574 // gentler and hand over a whole row of entries at a time
575 std::vector<size_type> column_indices;
576 column_indices.reserve(
577 entries_per_row.size() > 0 ?
578 *std::max_element(entries_per_row.begin(), entries_per_row.end()) :
579 0);
580 for (size_type row = 0; row < matrix.m(); ++row)
581 {
582 column_indices.resize(entries_per_row[row]);
583
584 size_type current_index = 0;
585 for (size_type col = 0; col < matrix.n(); ++col)
586 if (matrix(row, col) != 0)
587 {
588 column_indices[current_index] = col;
589 ++current_index;
590 }
591 else
592 // the (row,col) entry is zero; check if we need to add it
593 // anyway because it's the diagonal entry of a square
594 // matrix
595 if (matrix_is_square && (col == row))
596 {
597 column_indices[current_index] = row;
598 ++current_index;
599 }
600
601 // check that we really added the correct number of indices
602 Assert(current_index == entries_per_row[row], ExcInternalError());
603
604 // now bulk add all of these entries
605 add_entries(row, column_indices.begin(), column_indices.end(), true);
606 }
607
608 // finally compress
609 compress();
610}
611
612
613void
615 const size_type n,
616 const std::vector<unsigned int> &row_lengths)
617{
618 reinit(m, n, make_array_view(row_lengths));
619}
620
621
622
623bool
625{
626 // let's try to be on the safe side of life by using multiple possibilities
627 // in the check for emptiness... (sorry for this kludge -- emptying matrices
628 // and freeing memory was not present in the original implementation and I
629 // don't know at how many places I missed something in adding it, so I try
630 // to be cautious. wb)
631 if ((rowstart == nullptr) || (rows == 0) || (cols == 0))
632 {
633 Assert(rowstart == nullptr, ExcInternalError());
636 Assert(colnums == nullptr, ExcInternalError());
638
639 return true;
640 }
641 return false;
642}
643
644
645
648{
649 // if compress() has not yet been called, we can get the maximum number of
650 // elements per row using the stored value
651 if (!compressed)
652 return max_row_length;
653
654 // if compress() was called, we use a better algorithm which gives us a
655 // sharp bound
656 size_type m = 0;
657 for (size_type i = 1; i <= rows; ++i)
658 m = std::max(m, static_cast<size_type>(rowstart[i] - rowstart[i - 1]));
659
660 return m;
661}
662
663
664
667{
668 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
672
673 // let's see whether there is something in this line
674 if (rowstart[i] == rowstart[i + 1])
675 return invalid_entry;
676
677 // If special storage of diagonals was requested, we can get the diagonal
678 // element faster by this query.
679 if (store_diagonal_first_in_row && (i == j))
680 return rowstart[i];
681
682 // all other entries are sorted, so we can use a binary search algorithm
683 //
684 // note that the entries are only sorted upon compression, so this would
685 // fail for non-compressed sparsity patterns; however, that is why the
686 // Assertion is at the top of this function, so it may not be called for
687 // noncompressed structures.
688 const size_type *sorted_region_start =
690 &colnums[rowstart[i]]);
691 const size_type *const p =
692 Utilities::lower_bound<const size_type *>(sorted_region_start,
693 &colnums[rowstart[i + 1]],
694 j);
695 if ((p != &colnums[rowstart[i + 1]]) && (*p == j))
696 return (p - colnums.get());
697 else
698 return invalid_entry;
699}
700
701
702
703void
705{
706 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
710
711 for (std::size_t k = rowstart[i]; k < rowstart[i + 1]; ++k)
712 {
713 // entry already exists
714 if (colnums[k] == j)
715 return;
716 // empty entry found, put new entry here
717 if (colnums[k] == invalid_entry)
718 {
719 colnums[k] = j;
720 return;
721 }
722 }
723
724 // if we came thus far, something went wrong: there was not enough space in
725 // this line
726 Assert(false, ExcNotEnoughSpace(i, rowstart[i + 1] - rowstart[i]));
727}
728
729
730
731template <typename ForwardIterator>
732void
734 ForwardIterator begin,
735 ForwardIterator end,
736 const bool indices_are_sorted)
737{
738 if (indices_are_sorted == true)
739 {
740 if (begin != end)
741 {
742 ForwardIterator it = begin;
743 bool has_larger_entries = false;
744 // skip diagonal
745 std::size_t k = rowstart[row] + store_diagonal_first_in_row;
746 for (; k < rowstart[row + 1]; ++k)
747 if (colnums[k] == invalid_entry)
748 break;
749 else if (colnums[k] >= *it)
750 {
751 has_larger_entries = true;
752 break;
753 }
754 if (has_larger_entries == false)
755 for (; it != end; ++it)
756 {
757 if (store_diagonal_first_in_row && *it == row)
758 continue;
759 Assert(k <= rowstart[row + 1],
761 rowstart[row + 1] - rowstart[row]));
762 colnums[k++] = *it;
763 }
764 else
765 // cannot just append the new range at the end, forward to the
766 // other function
767 for (ForwardIterator p = begin; p != end; ++p)
768 add(row, *p);
769 }
770 }
771 else
772 {
773 // forward to the other function.
774 for (ForwardIterator it = begin; it != end; ++it)
775 add(row, *it);
776 }
777}
778
779
780bool
782{
783 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
786
787 for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k)
788 {
789 // entry already exists
790 if (colnums[k] == j)
791 return true;
792 }
793 return false;
794}
795
796
797
800{
801 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
804
805 for (size_type k = rowstart[i]; k < rowstart[i + 1]; ++k)
806 {
807 // entry exists
808 if (colnums[k] == j)
809 return k - rowstart[i];
810 }
812}
813
814
815
816std::pair<SparsityPatternBase::size_type, SparsityPatternBase::size_type>
818{
821
822 // first find the row in which the entry is located. for this note that the
823 // rowstart array indexes the global indices at which each row starts. since
824 // it is sorted, and since there is an element for the one-past-last row, we
825 // can simply use a bisection search on it
826 const size_type row =
827 (std::upper_bound(rowstart.get(), rowstart.get() + rows, global_index) -
828 rowstart.get() - 1);
829
830 // now, the column index is simple since that is what the colnums array
831 // stores:
832 const size_type col = colnums[global_index];
833
834 // so return the respective pair
835 return std::make_pair(row, col);
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 < 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
874void
875SparsityPatternBase::print(std::ostream &out) const
876{
877 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
878
879 AssertThrow(out.fail() == false, ExcIO());
880
881 for (size_type i = 0; i < rows; ++i)
882 {
883 out << '[' << i;
884 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
885 if (colnums[j] != invalid_entry)
886 out << ',' << colnums[j];
887 out << ']' << std::endl;
888 }
889
890 AssertThrow(out.fail() == false, ExcIO());
891}
892
893
894
895void
897{
898 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
899
900 AssertThrow(out.fail() == false, ExcIO());
901
902 for (size_type i = 0; i < rows; ++i)
903 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
904 if (colnums[j] != invalid_entry)
905 // while matrix entries are usually written (i,j), with i vertical and
906 // j horizontal, gnuplot output is x-y, that is we have to exchange
907 // the order of output
908 out << colnums[j] << " " << -static_cast<signed int>(i) << std::endl;
909
910 AssertThrow(out.fail() == false, ExcIO());
911}
912
913void
914SparsityPatternBase::print_svg(std::ostream &out) const
915{
916 const unsigned int m = this->n_rows();
917 const unsigned int n = this->n_cols();
918 out
919 << "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\" viewBox=\"0 0 "
920 << n + 2 << " " << m + 2
921 << " \">\n"
922 "<style type=\"text/css\" >\n"
923 " <![CDATA[\n"
924 " rect.pixel {\n"
925 " fill: #ff0000;\n"
926 " }\n"
927 " ]]>\n"
928 " </style>\n\n"
929 " <rect width=\""
930 << n + 2 << "\" height=\"" << m + 2
931 << "\" fill=\"rgb(128, 128, 128)\"/>\n"
932 " <rect x=\"1\" y=\"1\" width=\""
933 << n + 0.1 << "\" height=\"" << m + 0.1
934 << "\" fill=\"rgb(255, 255, 255)\"/>\n\n";
935
936 for (const auto &entry : *this)
937 {
938 out << " <rect class=\"pixel\" x=\"" << entry.column() + 1 << "\" y=\""
939 << entry.row() + 1 << "\" width=\".9\" height=\".9\"/>\n";
940 }
941 out << "</svg>" << std::endl;
942}
943
944
945
948{
949 Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
950 size_type b = 0;
951 for (size_type i = 0; i < rows; ++i)
952 for (size_type j = rowstart[i]; j < rowstart[i + 1]; ++j)
953 if (colnums[j] != invalid_entry)
954 {
955 if (static_cast<size_type>(
956 std::abs(static_cast<int>(i - colnums[j]))) > b)
957 b = std::abs(static_cast<signed int>(i - colnums[j]));
958 }
959 else
960 // leave if at the end of the entries of this line
961 break;
962 return b;
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 << ' ' << rows << ' ' << cols << ' ' << max_vec_len
973 << ' ' << 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
1037std::size_t
1039{
1040 return sizeof(*this) + SparsityPatternBase::memory_consumption();
1041}
1042
1043
1044
1045#ifndef DOXYGEN
1046// explicit instantiations
1047template void
1048SparsityPattern::copy_from<float>(const FullMatrix<float> &);
1049template void
1050SparsityPattern::copy_from<double>(const FullMatrix<double> &);
1051
1052template void
1053SparsityPattern::add_entries<const SparsityPattern::size_type *>(
1054 const size_type,
1055 const size_type *,
1056 const size_type *,
1057 const bool);
1058# ifndef DEAL_II_VECTOR_ITERATOR_IS_POINTER
1059template void
1061 std::vector<SparsityPattern::size_type>::const_iterator>(
1062 const size_type,
1063 std::vector<size_type>::const_iterator,
1064 std::vector<size_type>::const_iterator,
1065 const bool);
1066# endif
1067template void
1068SparsityPattern::add_entries<std::vector<SparsityPattern::size_type>::iterator>(
1069 const size_type,
1070 std::vector<size_type>::iterator,
1071 std::vector<size_type>::iterator,
1072 const bool);
1073#endif
1074
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:699
iterator begin() const
Definition: array_view.h:585
iterator end() const
Definition: array_view.h:594
std::size_t size() const
Definition: array_view.h:576
std::size_t n_nonzero_elements() const
unsigned int max_row_length
void print_svg(std::ostream &out) const
void add(const size_type i, const size_type j)
void reinit(const size_type m, const size_type n, const unsigned int max_per_row)
unsigned int row_length(const size_type row) const
size_type n_rows() const
std::unique_ptr< size_type[]> colnums
std::unique_ptr< std::size_t[]> rowstart
void print(std::ostream &out) const
bool exists(const size_type i, const size_type j) const
iterator end() const
void print_gnuplot(std::ostream &out) const
static const size_type invalid_entry
size_type bandwidth() const
bool is_compressed() const
size_type max_entries_per_row() const
std::pair< size_type, size_type > matrix_position(const std::size_t global_index) const
std::size_t memory_consumption() const
iterator begin() const
size_type row_position(const size_type i, const size_type j) const
size_type n_cols() const
void block_write(std::ostream &out) const
void add_entries(const size_type row, ForwardIterator begin, ForwardIterator end, const bool indices_are_sorted=false)
SparsityPattern & operator=(const SparsityPattern &)
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
virtual void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths) override
size_type operator()(const size_type i, const size_type j) const
std::size_t memory_consumption() const
void block_read(std::istream &in)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcEmptyObject()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotEnoughSpace(int arg1, int arg2)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcMatrixIsCompressed()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
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
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1153
void copy(const T *begin, const T *end, U *dest)
const types::global_dof_index invalid_size_type
Definition: types.h:210
::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 > &)