Reference documentation for deal.II version Git e3a3ec7800 2020-08-07 14:08:19 +0200
\(\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 - 2020 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 
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 
65 {
66  (void)s;
67  Assert(s.empty(),
68  ExcMessage(
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)
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)
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)
121  : SparsityPattern()
122 {
123  Assert(original.rows == original.cols, ExcNotQuadratic());
124  Assert(original.is_compressed(), ExcNotCompressed());
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(),
202  ExcMessage(
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 
216 void
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 
228 void
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
255  store_diagonal_first_in_row = (m == n);
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)),
321  ExcInternalError());
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 
337 void
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]],
357  &colnums[rowstart[rows]],
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),
404  ExcInternalError());
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]),
415  ExcInternalError());
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]),
420  ExcInternalError());
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 
441 void
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.
489 void
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 
549 template <typename number>
550 void
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 
613 void
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 
623 bool
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());
634  Assert(rows == 0, ExcInternalError());
635  Assert(cols == 0, 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 
703 void
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 
731 template <typename ForwardIterator>
732 void
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],
760  ExcNotEnoughSpace(row,
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 
780 bool
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 
816 std::pair<SparsityPatternBase::size_type, SparsityPatternBase::size_type>
818 {
819  Assert(compressed == true, ExcNotCompressed());
820  AssertIndexRange(global_index, 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::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 std::size_t
1040 {
1041  return sizeof(*this) + SparsityPatternBase::memory_consumption();
1042 }
1043 
1044 
1045 
1046 #ifndef DOXYGEN
1047 // explicit instantiations
1048 template void
1049 SparsityPattern::copy_from<float>(const FullMatrix<float> &);
1050 template void
1051 SparsityPattern::copy_from<double>(const FullMatrix<double> &);
1052 
1053 template void
1054 SparsityPattern::add_entries<const SparsityPattern::size_type *>(
1055  const size_type,
1056  const size_type *,
1057  const size_type *,
1058  const bool);
1059 # ifndef DEAL_II_VECTOR_ITERATOR_IS_POINTER
1060 template void
1062  std::vector<SparsityPattern::size_type>::const_iterator>(
1063  const size_type,
1064  std::vector<size_type>::const_iterator,
1065  std::vector<size_type>::const_iterator,
1066  const bool);
1067 # endif
1068 template void
1069 SparsityPattern::add_entries<std::vector<SparsityPattern::size_type>::iterator>(
1070  const size_type,
1071  std::vector<size_type>::iterator,
1072  std::vector<size_type>::iterator,
1073  const bool);
1074 #endif
1075 
Iterator lower_bound(Iterator first, Iterator last, const T &val)
Definition: utilities.h:1075
void block_write(std::ostream &out) const
size_type m() const
const types::global_dof_index invalid_size_type
Definition: types.h:205
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:1568
unsigned int max_row_length
static ::ExceptionBase & ExcMatrixIsCompressed()
iterator begin() const
Contents is actually a matrix.
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:560
std::size_t memory_consumption() const
#define AssertIndexRange(index, range)
Definition: exceptions.h:1636
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
static const char U
size_type bandwidth() const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1521
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:542
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)
size_type n() const
SparsityPatternBase::size_type size_type
#define Assert(cond, exc)
Definition: exceptions.h:1411
const IndexSet & row_index_set() const
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:665
std::size_t memory_consumption() const
void add(const size_type i, const size_type j)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:363
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
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
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
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:362
T min(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcEmptyObject()
void block_read(std::istream &in)
TrilinosWrappers::types::int_type global_index(const Epetra_BlockMap &map, const ::types::global_dof_index i)
static const size_type invalid_entry
iterator begin() const
Definition: array_view.h:551
bool exists(const size_type i, const size_type j) const
unsigned int row_length(const size_type row) const
void copy(const T *begin, const T *end, U *dest)
T max(const T &t, const MPI_Comm &mpi_communicator)
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