Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
chunk_sparsity_pattern.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2008 - 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/lac/chunk_sparsity_pattern.h>
18 #include <deal.II/lac/dynamic_sparsity_pattern.h>
19 #include <deal.II/lac/full_matrix.h>
20 
21 
22 DEAL_II_NAMESPACE_OPEN
23 
24 
26 {
27  reinit(0, 0, 0, 1);
28 }
29 
30 
31 
33  : Subscriptor()
36 {
37  Assert(s.rows == 0 && s.cols == 0,
38  ExcMessage(
39  "This constructor can only be called if the provided argument "
40  "is the sparsity pattern for an empty matrix. This constructor can "
41  "not be used to copy-construct a non-empty sparsity pattern."));
42 
43  reinit(0, 0, 0, chunk_size);
44 }
45 
46 
47 
49  const size_type n,
50  const size_type max_per_row,
51  const size_type chunk_size)
52 {
53  Assert(chunk_size > 0, ExcInvalidNumber(chunk_size));
54 
55  reinit(m, n, max_per_row, chunk_size);
56 }
57 
58 
59 
61  const size_type m,
62  const size_type n,
63  const std::vector<size_type> &row_lengths,
64  const size_type chunk_size)
65 {
66  Assert(chunk_size > 0, ExcInvalidNumber(chunk_size));
67 
68  reinit(m, n, row_lengths, chunk_size);
69 }
70 
71 
72 
74  const size_type max_per_row,
75  const size_type chunk_size)
76 {
77  reinit(n, n, max_per_row, chunk_size);
78 }
79 
80 
81 
83  const size_type m,
84  const std::vector<size_type> &row_lengths,
85  const size_type chunk_size)
86 {
87  Assert(chunk_size > 0, ExcInvalidNumber(chunk_size));
88 
89  reinit(m, m, row_lengths, chunk_size);
90 }
91 
92 
93 
96 {
97  Assert(s.rows == 0 && s.cols == 0,
98  ExcMessage(
99  "This operator can only be called if the provided argument "
100  "is the sparsity pattern for an empty matrix. This operator can "
101  "not be used to copy a non-empty sparsity pattern."));
102 
103  Assert(rows == 0 && cols == 0,
104  ExcMessage("This operator can only be called if the current object is "
105  "empty."));
106 
107  // perform the checks in the underlying object as well
109 
110  return *this;
111 }
112 
113 
114 
115 void
117  const size_type n,
118  const size_type max_per_row,
119  const size_type chunk_size)
120 {
121  Assert(chunk_size > 0, ExcInvalidNumber(chunk_size));
122 
123  // simply map this function to the other @p{reinit} function
124  const std::vector<size_type> row_lengths(m, max_per_row);
125  reinit(m, n, row_lengths, chunk_size);
126 }
127 
128 
129 
130 void
132  const size_type n,
133  const ArrayView<const size_type> &row_lengths,
134  const size_type chunk_size)
135 {
136  Assert(row_lengths.size() == m, ExcInvalidNumber(m));
137  Assert(chunk_size > 0, ExcInvalidNumber(chunk_size));
138 
139  rows = m;
140  cols = n;
141 
142  this->chunk_size = chunk_size;
143 
144  // pass down to the necessary information to the underlying object. we need
145  // to calculate how many chunks we need: we need to round up (m/chunk_size)
146  // and (n/chunk_size). rounding up in integer arithmetic equals
147  // ((m+chunk_size-1)/chunk_size):
148  const size_type m_chunks = (m + chunk_size - 1) / chunk_size,
149  n_chunks = (n + chunk_size - 1) / chunk_size;
150 
151  // compute the maximum number of chunks in each row. the passed array
152  // denotes the number of entries in each row of the big matrix -- in the
153  // worst case, these are all in independent chunks, so we have to calculate
154  // it as follows (as an example: let chunk_size==2, row_lengths={2,2,...},
155  // and entries in row zero at columns {0,2} and for row one at {4,6} -->
156  // we'll need 4 chunks for the first chunk row!) :
157  std::vector<unsigned int> chunk_row_lengths(m_chunks, 0);
158  for (size_type i = 0; i < m; ++i)
159  chunk_row_lengths[i / chunk_size] += row_lengths[i];
160 
161  // for the case that the reduced sparsity pattern optimizes the diagonal but
162  // the actual sparsity pattern does not, need to take one more entry in the
163  // row to fit the user-required entry
164  if (m != n && m_chunks == n_chunks)
165  for (unsigned int i = 0; i < m_chunks; ++i)
166  ++chunk_row_lengths[i];
167 
168  sparsity_pattern.reinit(m_chunks, n_chunks, chunk_row_lengths);
169 }
170 
171 
172 
173 void
175 {
177 }
178 
179 
180 
181 template <typename SparsityPatternType>
182 void
183 ChunkSparsityPattern::copy_from(const SparsityPatternType &dsp,
184  const size_type chunk_size)
185 {
186  Assert(chunk_size > 0, ExcInvalidNumber(chunk_size));
187  this->chunk_size = chunk_size;
188  rows = dsp.n_rows();
189  cols = dsp.n_cols();
190 
191  // simple case: just use the given sparsity pattern
192  if (chunk_size == 1)
193  {
195  return;
196  }
197 
198  // create a temporary compressed sparsity pattern that collects all entries
199  // from the input sparsity pattern and then initialize the underlying small
200  // sparsity pattern
201  const size_type m_chunks = (dsp.n_rows() + chunk_size - 1) / chunk_size,
202  n_chunks = (dsp.n_cols() + chunk_size - 1) / chunk_size;
203  DynamicSparsityPattern temporary_sp(m_chunks, n_chunks);
204 
205  for (size_type row = 0; row < dsp.n_rows(); ++row)
206  {
207  const size_type reduced_row = row / chunk_size;
208 
209  // TODO: This could be made more efficient if we cached the
210  // previous column and only called add() if the previous and the
211  // current column lead to different chunk columns
212  for (typename SparsityPatternType::iterator col_num = dsp.begin(row);
213  col_num != dsp.end(row);
214  ++col_num)
215  temporary_sp.add(reduced_row, col_num->column() / chunk_size);
216  }
217 
218  sparsity_pattern.copy_from(temporary_sp);
219 }
220 
221 
222 
223 template <typename number>
224 void
226  const size_type chunk_size)
227 {
228  Assert(chunk_size > 0, ExcInvalidNumber(chunk_size));
229 
230  // count number of entries per row, then initialize the underlying sparsity
231  // pattern. remember to also allocate space for the diagonal entry (if that
232  // hasn't happened yet) if m==n since we always allocate that for diagonal
233  // matrices
234  std::vector<size_type> entries_per_row(matrix.m(), 0);
235  for (size_type row = 0; row < matrix.m(); ++row)
236  {
237  for (size_type col = 0; col < matrix.n(); ++col)
238  if (matrix(row, col) != 0)
239  ++entries_per_row[row];
240 
241  if ((matrix.m() == matrix.n()) && (matrix(row, row) == 0))
242  ++entries_per_row[row];
243  }
244 
245  reinit(matrix.m(), matrix.n(), entries_per_row, chunk_size);
246 
247  // then actually fill it
248  for (size_type row = 0; row < matrix.m(); ++row)
249  for (size_type col = 0; col < matrix.n(); ++col)
250  if (matrix(row, col) != 0)
251  add(row, col);
252 
253  // finally compress
254  compress();
255 }
256 
257 
258 
259 void
261  const size_type n,
262  const std::vector<size_type> &row_lengths,
263  const size_type chunk_size)
264 {
265  Assert(chunk_size > 0, ExcInvalidNumber(chunk_size));
266 
267  reinit(m, n, make_array_view(row_lengths), chunk_size);
268 }
269 
270 
271 
272 namespace internal
273 {
274  namespace
275  {
276  template <typename SparsityPatternType>
277  void
278  copy_sparsity(const SparsityPatternType &src, SparsityPattern &dst)
279  {
280  dst.copy_from(src);
281  }
282 
283  void
284  copy_sparsity(const SparsityPattern &src, SparsityPattern &dst)
285  {
286  dst = src;
287  }
288  } // namespace
289 } // namespace internal
290 
291 
292 
293 template <typename Sparsity>
294 void
296  const size_type n,
297  const Sparsity &sparsity_pattern_for_chunks,
298  const size_type chunk_size_in,
299  const bool)
300 {
301  Assert(m > (sparsity_pattern_for_chunks.n_rows() - 1) * chunk_size_in &&
302  m <= sparsity_pattern_for_chunks.n_rows() * chunk_size_in,
303  ExcMessage("Number of rows m is not compatible with chunk size "
304  "and number of rows in sparsity pattern for the chunks."));
305  Assert(n > (sparsity_pattern_for_chunks.n_cols() - 1) * chunk_size_in &&
306  n <= sparsity_pattern_for_chunks.n_cols() * chunk_size_in,
307  ExcMessage(
308  "Number of columns m is not compatible with chunk size "
309  "and number of columns in sparsity pattern for the chunks."));
310 
311  internal::copy_sparsity(sparsity_pattern_for_chunks, sparsity_pattern);
312  chunk_size = chunk_size_in;
313  rows = m;
314  cols = n;
315 }
316 
317 
318 
319 bool
321 {
322  return sparsity_pattern.empty();
323 }
324 
325 
326 
329 {
331 }
332 
333 
334 
335 void
337 {
338  Assert(i < rows, ExcInvalidIndex(i, rows));
339  Assert(j < cols, ExcInvalidIndex(j, cols));
340 
342 }
343 
344 
345 bool
347 {
348  Assert(i < rows, ExcIndexRange(i, 0, rows));
349  Assert(j < cols, ExcIndexRange(j, 0, cols));
350 
351  return sparsity_pattern.exists(i / chunk_size, j / chunk_size);
352 }
353 
354 
355 
356 void
358 {
359  // matrix must be square. note that the for some matrix sizes, the current
360  // sparsity pattern may not be square even if the underlying sparsity
361  // pattern is (e.g. a 10x11 matrix with chunk_size 4)
363 
365 }
366 
367 
368 
371 {
372  Assert(i < rows, ExcIndexRange(i, 0, rows));
373 
374  // find out if we did padding and if this row is affected by it
375  if (n_cols() % chunk_size == 0)
377  else
378  // if columns don't align, then just iterate over all chunks and see
379  // what this leads to
380  {
383  end =
385  unsigned int n = 0;
386  for (; p != end; ++p)
387  if (p->column() != sparsity_pattern.n_cols() - 1)
388  n += chunk_size;
389  else
390  n += (n_cols() % chunk_size);
391  return n;
392  }
393 }
394 
395 
396 
399 {
400  if ((n_rows() % chunk_size == 0) && (n_cols() % chunk_size == 0))
402  else
403  // some of the chunks reach beyond the extent of this matrix. this
404  // requires a somewhat more complicated computations, in particular if the
405  // columns don't align
406  {
407  if ((n_rows() % chunk_size != 0) && (n_cols() % chunk_size == 0))
408  {
409  // columns align with chunks, but not rows
410  size_type n =
412  n -= (sparsity_pattern.n_rows() * chunk_size - n_rows()) *
414  chunk_size;
415  return n;
416  }
417 
418  else
419  {
420  // if columns don't align, then just iterate over all chunks and see
421  // what this leads to. follow the advice in the documentation of the
422  // sparsity pattern iterators to do the loop over individual rows,
423  // rather than all elements
424  size_type n = 0;
425 
426  for (size_type row = 0; row < sparsity_pattern.n_rows(); ++row)
427  {
429  for (; p != sparsity_pattern.end(row); ++p)
430  if ((row != sparsity_pattern.n_rows() - 1) &&
431  (p->column() != sparsity_pattern.n_cols() - 1))
432  n += chunk_size * chunk_size;
433  else if ((row == sparsity_pattern.n_rows() - 1) &&
434  (p->column() != sparsity_pattern.n_cols() - 1))
435  // last chunk row, but not last chunk column. only a smaller
436  // number (n_rows % chunk_size) of rows actually exist
437  n += (n_rows() % chunk_size) * chunk_size;
438  else if ((row != sparsity_pattern.n_rows() - 1) &&
439  (p->column() == sparsity_pattern.n_cols() - 1))
440  // last chunk column, but not row
441  n += (n_cols() % chunk_size) * chunk_size;
442  else
443  // bottom right chunk
444  n += (n_cols() % chunk_size) * (n_rows() % chunk_size);
445  }
446 
447  return n;
448  }
449  }
450 }
451 
452 
453 
454 void
455 ChunkSparsityPattern::print(std::ostream &out) const
456 {
457  Assert((sparsity_pattern.rowstart != nullptr) &&
458  (sparsity_pattern.colnums != nullptr),
459  ExcEmptyObject());
460 
461  AssertThrow(out, ExcIO());
462 
463  for (size_type i = 0; i < sparsity_pattern.rows; ++i)
464  for (size_type d = 0; (d < chunk_size) && (i * chunk_size + d < n_rows());
465  ++d)
466  {
467  out << '[' << i * chunk_size + d;
468  for (size_type j = sparsity_pattern.rowstart[i];
469  j < sparsity_pattern.rowstart[i + 1];
470  ++j)
472  for (size_type e = 0;
473  ((e < chunk_size) &&
474  (sparsity_pattern.colnums[j] * chunk_size + e < n_cols()));
475  ++e)
476  out << ',' << sparsity_pattern.colnums[j] * chunk_size + e;
477  out << ']' << std::endl;
478  }
479 
480  AssertThrow(out, ExcIO());
481 }
482 
483 
484 
485 void
486 ChunkSparsityPattern::print_gnuplot(std::ostream &out) const
487 {
488  Assert((sparsity_pattern.rowstart != nullptr) &&
489  (sparsity_pattern.colnums != nullptr),
490  ExcEmptyObject());
491 
492  AssertThrow(out, ExcIO());
493 
494  // for each entry in the underlying sparsity pattern, repeat everything
495  // chunk_size x chunk_size times
496  for (size_type i = 0; i < sparsity_pattern.rows; ++i)
497  for (size_type j = sparsity_pattern.rowstart[i];
498  j < sparsity_pattern.rowstart[i + 1];
499  ++j)
501  for (size_type d = 0;
502  ((d < chunk_size) &&
503  (sparsity_pattern.colnums[j] * chunk_size + d < n_cols()));
504  ++d)
505  for (size_type e = 0;
506  (e < chunk_size) && (i * chunk_size + e < n_rows());
507  ++e)
508  // while matrix entries are usually written (i,j), with i vertical
509  // and j horizontal, gnuplot output is x-y, that is we have to
510  // exchange the order of output
511  out << sparsity_pattern.colnums[j] * chunk_size + d << " "
512  << -static_cast<signed int>(i * chunk_size + e) << std::endl;
513 
514  AssertThrow(out, ExcIO());
515 }
516 
517 
518 
521 {
522  // calculate the bandwidth from that of the underlying sparsity
523  // pattern. note that even if the bandwidth of that is zero, then the
524  // bandwidth of the chunky pattern is chunk_size-1, if it is 1 then the
525  // chunky pattern has chunk_size+(chunk_size-1), etc
526  //
527  // we'll cut it off at max(n(),m())
528  return std::min(sparsity_pattern.bandwidth() * chunk_size + (chunk_size - 1),
529  std::max(n_rows(), n_cols()));
530 }
531 
532 
533 
534 bool
536 {
537  if (chunk_size == 1)
539  else
540  return false;
541 }
542 
543 
544 
545 void
546 ChunkSparsityPattern::block_write(std::ostream &out) const
547 {
548  AssertThrow(out, ExcIO());
549 
550  // first the simple objects, bracketed in [...]
551  out << '[' << rows << ' ' << cols << ' ' << chunk_size << ' ' << "][";
552  // then the underlying sparsity pattern
554  out << ']';
555 
556  AssertThrow(out, ExcIO());
557 }
558 
559 
560 
561 void
563 {
564  AssertThrow(in, ExcIO());
565 
566  char c;
567 
568  // first read in simple data
569  in >> c;
570  AssertThrow(c == '[', ExcIO());
571  in >> rows >> cols >> chunk_size;
572 
573  in >> c;
574  AssertThrow(c == ']', ExcIO());
575  in >> c;
576  AssertThrow(c == '[', ExcIO());
577 
578  // then read the underlying sparsity pattern
580 
581  in >> c;
582  AssertThrow(c == ']', ExcIO());
583 }
584 
585 
586 
587 std::size_t
589 {
590  return (sizeof(*this) + sparsity_pattern.memory_consumption());
591 }
592 
593 
594 
595 // explicit instantiations
596 template void
597 ChunkSparsityPattern::copy_from<DynamicSparsityPattern>(
598  const DynamicSparsityPattern &,
599  const size_type);
600 template void
601 ChunkSparsityPattern::create_from<SparsityPattern>(const size_type,
602  const size_type,
603  const SparsityPattern &,
604  const size_type,
605  const bool);
606 template void
607 ChunkSparsityPattern::create_from<DynamicSparsityPattern>(
608  const size_type,
609  const size_type,
610  const DynamicSparsityPattern &,
611  const size_type,
612  const bool);
613 template void
614 ChunkSparsityPattern::copy_from<float>(const FullMatrix<float> &,
615  const size_type);
616 template void
617 ChunkSparsityPattern::copy_from<double>(const FullMatrix<double> &,
618  const size_type);
619 
620 DEAL_II_NAMESPACE_CLOSE
void block_write(std::ostream &out) const
size_type m() const
void block_write(std::ostream &out) const
iterator end() const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end, const size_type chunk_size)
std::size_t memory_consumption() const
static ::ExceptionBase & ExcEmptyObject()
iterator begin() const
static ::ExceptionBase & ExcIO()
types::global_dof_index size_type
void add(const size_type i, const size_type j)
iterator end() const
std::size_t memory_consumption() const
void block_read(std::istream &in)
bool exists(const size_type i, const size_type j) const
size_type bandwidth() const
void create_from(const size_type m, const size_type n, const Sparsity &sparsity_pattern_for_chunks, const size_type chunk_size, const bool optimize_diagonal=true)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
void add(const size_type i, const size_type j)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcInvalidNumber(size_type arg1)
std::size_t size() const
Definition: array_view.h:471
SparsityPattern sparsity_pattern
std::unique_ptr< std::size_t[]> rowstart
size_type max_entries_per_row() const
size_type n_cols() 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
#define Assert(cond, exc)
Definition: exceptions.h:1407
void add(const size_type i, const size_type j)
size_type row_length(const size_type row) const
size_type n_rows() const
size_type n_nonzero_elements() const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
std::size_t n_nonzero_elements() const
ChunkSparsityPattern & operator=(const ChunkSparsityPattern &)
bool exists(const size_type i, const size_type j) const
std::unique_ptr< size_type[]> colnums
static ::ExceptionBase & ExcNotQuadratic()
bool stores_only_added_elements() const
void block_read(std::istream &in)
static const size_type invalid_entry
void print(std::ostream &out) const
unsigned int row_length(const size_type row) const
void print_gnuplot(std::ostream &out) const
void reinit(const size_type m, const size_type n, const size_type max_per_row, const size_type chunk_size)
size_type n_rows() const
size_type n_cols() const
size_type max_entries_per_row() const