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