Reference documentation for deal.II version GIT 58febcd5cf 2023-09-30 20:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
17 #include <deal.II/base/utilities.h>
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 
52  : SparsityPattern()
53 {
54  (void)s;
55  Assert(s.empty(),
56  ExcMessage(
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)
69  : SparsityPattern()
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)
79  : SparsityPattern()
80 {
81  reinit(m, n, row_lengths);
82 }
83 
84 
85 
87  const unsigned int max_per_row)
88  : SparsityPattern()
89 {
90  reinit(m, m, max_per_row);
91 }
92 
93 
94 
96  const std::vector<unsigned int> &row_lengths)
97  : SparsityPattern()
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)
107  : SparsityPattern()
108 {
109  Assert(original.n_rows() == original.n_cols(), ExcNotQuadratic());
110  Assert(original.is_compressed(), ExcNotCompressed());
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(),
188  ExcMessage(
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 
202 void
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
227  store_diagonal_first_in_row = (m == n);
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.empty() ?
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)),
293  ExcInternalError());
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 
309 void
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 
319 void
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 
331 void
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]],
351  &colnums[rowstart[rows]],
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),
398  ExcInternalError());
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]),
409  ExcInternalError());
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]),
414  ExcInternalError());
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 
435 void
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.
483 void
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 
543 template <typename number>
544 void
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 
608 bool
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 
631 bool
633 {
634  Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
635  AssertIndexRange(i, n_rows());
636  AssertIndexRange(j, n_cols());
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 
649 std::pair<SparsityPattern::size_type, SparsityPattern::size_type>
651 {
652  Assert(compressed == true, ExcNotCompressed());
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());
717  AssertIndexRange(i, n_rows());
718  AssertIndexRange(j, n_cols());
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 
751 void
753 {
754  Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
755  AssertIndexRange(i, n_rows());
756  AssertIndexRange(j, n_cols());
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 
779 template <typename ForwardIterator>
780 void
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],
810  ExcNotEnoughSpace(row,
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 
831 void
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 
841 void
843 {
844  Assert((rowstart != nullptr) && (colnums != nullptr), ExcEmptyObject());
846  // Note that we only require a quadratic matrix here, no special treatment
847  // of diagonals
848  Assert(n_rows() == n_cols(), ExcNotQuadratic());
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());
879  AssertIndexRange(i, n_rows());
880  AssertIndexRange(j, n_cols());
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 
893 void
894 SparsityPattern::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 
914 void
915 SparsityPattern::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 
934 void
935 SparsityPattern::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 
967 void
968 SparsityPattern::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 << ' '
975  << store_diagonal_first_in_row << "][";
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 
991 void
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 
1030 std::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
1041 template void
1042 SparsityPattern::copy_from<float>(const FullMatrix<float> &);
1043 template void
1044 SparsityPattern::copy_from<double>(const FullMatrix<double> &);
1045 
1046 template void
1047 SparsityPattern::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
1053 template 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
1061 template void
1062 SparsityPattern::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:838
iterator begin() const
Definition: array_view.h:591
bool empty() const
Definition: array_view.h:582
iterator end() const
Definition: array_view.h:600
std::size_t size() const
Definition: array_view.h:573
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
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
virtual void add_entries(const ArrayView< const std::pair< size_type, size_type >> &entries)
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
bool empty() const
unsigned int max_row_length
iterator end() const
size_type operator()(const size_type i, const size_type j) const
types::global_dof_index size_type
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:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotCompressed()
static ::ExceptionBase & ExcNotEnoughSpace(int arg1, int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1789
static ::ExceptionBase & ExcEmptyObject()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
static ::ExceptionBase & ExcIO()
static ::ExceptionBase & ExcMatrixIsCompressed()
static ::ExceptionBase & ExcMessage(std::string arg1)
static ::ExceptionBase & ExcNotQuadratic()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
@ 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 * begin(VectorType &V)
VectorType::value_type * end(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:1008
void copy(const T *begin, const T *end, U *dest)
const types::global_dof_index invalid_size_type
Definition: types.h:222