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