Reference documentation for deal.II version Git 0943bc0020 2021-10-22 11:23:14 -0400
\(\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\}}\)
block_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 
18 
20 
22 
23 
24 template <class SparsityPatternBase>
26  : rows(0)
27  , columns(0)
28 {}
29 
30 
31 
32 template <class SparsityPatternBase>
34  const size_type n_block_rows,
35  const size_type n_block_columns)
36  : rows(0)
37  , columns(0)
38 {
39  reinit(n_block_rows, n_block_columns);
40 }
41 
42 
43 
44 template <class SparsityPatternBase>
46  const BlockSparsityPatternBase &s)
47  : Subscriptor()
48  , rows(0)
49  , columns(0)
50 {
51  (void)s;
52  Assert(s.rows == 0 && s.columns == 0,
53  ExcMessage(
54  "This constructor can only be called if the provided argument "
55  "is the sparsity pattern for an empty matrix. This constructor can "
56  "not be used to copy-construct a non-empty sparsity pattern."));
57 }
58 
59 
60 
61 template <class SparsityPatternBase>
63 {
64  // clear all memory
65  try
66  {
67  reinit(0, 0);
68  }
69  catch (...)
70  {}
71 }
72 
73 
74 
75 template <class SparsityPatternBase>
76 void
78  const size_type n_block_rows,
79  const size_type n_block_columns)
80 {
81  // delete previous content and
82  // clean the sub_objects array
83  // completely
84  for (size_type i = 0; i < rows; ++i)
85  for (size_type j = 0; j < columns; ++j)
86  {
87  SparsityPatternBase *sp = sub_objects[i][j];
88  sub_objects[i][j] = nullptr;
89  delete sp;
90  }
91  sub_objects.reinit(0, 0);
92 
93  // then set new sizes
94  rows = n_block_rows;
95  columns = n_block_columns;
96  sub_objects.reinit(rows, columns);
97 
98  // allocate new objects
99  for (size_type i = 0; i < rows; ++i)
100  for (size_type j = 0; j < columns; ++j)
101  {
103  sub_objects[i][j] = p;
104  }
105 }
106 
107 
108 template <class SparsityPatternBase>
112 {
113  Assert(rows == bsp.rows, ExcDimensionMismatch(rows, bsp.rows));
115  // copy objects
116  for (size_type i = 0; i < rows; ++i)
117  for (size_type j = 0; j < columns; ++j)
118  *sub_objects[i][j] = *bsp.sub_objects[i][j];
119  // update index objects
120  collect_sizes();
121 
122  return *this;
123 }
124 
125 
126 
127 template <class SparsityPatternBase>
128 void
130 {
131  std::vector<size_type> row_sizes(rows);
132  std::vector<size_type> col_sizes(columns);
133 
134  // first find out the row sizes
135  // from the first block column
136  for (size_type r = 0; r < rows; ++r)
137  row_sizes[r] = sub_objects[r][0]->n_rows();
138  // then check that the following
139  // block columns have the same
140  // sizes
141  for (size_type c = 1; c < columns; ++c)
142  for (size_type r = 0; r < rows; ++r)
143  Assert(row_sizes[r] == sub_objects[r][c]->n_rows(),
144  ExcIncompatibleRowNumbers(r, 0, r, c));
145 
146  // finally initialize the row
147  // indices with this array
148  row_indices.reinit(row_sizes);
149 
150 
151  // then do the same with the columns
152  for (size_type c = 0; c < columns; ++c)
153  col_sizes[c] = sub_objects[0][c]->n_cols();
154  for (size_type r = 1; r < rows; ++r)
155  for (size_type c = 0; c < columns; ++c)
156  Assert(col_sizes[c] == sub_objects[r][c]->n_cols(),
157  ExcIncompatibleRowNumbers(0, c, r, c));
158 
159  // finally initialize the row
160  // indices with this array
161  column_indices.reinit(col_sizes);
162 }
163 
164 
165 
166 template <class SparsityPatternBase>
167 void
169 {
170  for (size_type i = 0; i < rows; ++i)
171  for (size_type j = 0; j < columns; ++j)
172  sub_objects[i][j]->compress();
173 }
174 
175 
176 
177 template <class SparsityPatternBase>
178 bool
180 {
181  for (size_type i = 0; i < rows; ++i)
182  for (size_type j = 0; j < columns; ++j)
183  if (sub_objects[i][j]->empty() == false)
184  return false;
185  return true;
186 }
187 
188 
189 
190 template <class SparsityPatternBase>
193 {
194  size_type max_entries = 0;
195  for (size_type block_row = 0; block_row < rows; ++block_row)
196  {
197  size_type this_row = 0;
198  for (size_type c = 0; c < columns; ++c)
199  this_row += sub_objects[block_row][c]->max_entries_per_row();
200 
201  if (this_row > max_entries)
202  max_entries = this_row;
203  }
204  return max_entries;
205 }
206 
207 
208 
209 template <class SparsityPatternBase>
212 {
213  // only count in first column, since
214  // all rows should be equivalent
215  size_type count = 0;
216  for (size_type r = 0; r < rows; ++r)
217  count += sub_objects[r][0]->n_rows();
218  return count;
219 }
220 
221 
222 
223 template <class SparsityPatternBase>
226 {
227  // only count in first row, since
228  // all rows should be equivalent
229  size_type count = 0;
230  for (size_type c = 0; c < columns; ++c)
231  count += sub_objects[0][c]->n_cols();
232  return count;
233 }
234 
235 
236 
237 template <class SparsityPatternBase>
240 {
241  size_type count = 0;
242  for (size_type i = 0; i < rows; ++i)
243  for (size_type j = 0; j < columns; ++j)
244  count += sub_objects[i][j]->n_nonzero_elements();
245  return count;
246 }
247 
248 
249 
250 template <class SparsityPatternBase>
251 void
253 {
254  size_type k = 0;
255  for (size_type ib = 0; ib < n_block_rows(); ++ib)
256  {
257  for (size_type i = 0; i < block(ib, 0).n_rows(); ++i)
258  {
259  out << '[' << i + k;
260  size_type l = 0;
261  for (size_type jb = 0; jb < n_block_cols(); ++jb)
262  {
263  const SparsityPatternBase &b = block(ib, jb);
264  for (size_type j = 0; j < b.n_cols(); ++j)
265  if (b.exists(i, j))
266  out << ',' << l + j;
267  l += b.n_cols();
268  }
269  out << ']' << std::endl;
270  }
271  k += block(ib, 0).n_rows();
272  }
273 }
274 
275 
276 #ifndef DOXYGEN
277 template <>
278 void
280 {
281  size_type k = 0;
282  for (size_type ib = 0; ib < n_block_rows(); ++ib)
283  {
284  for (size_type i = 0; i < block(ib, 0).n_rows(); ++i)
285  {
286  out << '[' << i + k;
287  size_type l = 0;
288  for (size_type jb = 0; jb < n_block_cols(); ++jb)
289  {
290  const DynamicSparsityPattern &b = block(ib, jb);
291  if (b.row_index_set().size() == 0 ||
292  b.row_index_set().is_element(i))
293  for (size_type j = 0; j < b.n_cols(); ++j)
294  if (b.exists(i, j))
295  out << ',' << l + j;
296  l += b.n_cols();
297  }
298  out << ']' << std::endl;
299  }
300  k += block(ib, 0).n_rows();
301  }
302 }
303 #endif
304 
305 
306 
307 template <class SparsityPatternBase>
308 void
310  std::ostream &out) const
311 {
312  size_type k = 0;
313  for (size_type ib = 0; ib < n_block_rows(); ++ib)
314  {
315  for (size_type i = 0; i < block(ib, 0).n_rows(); ++i)
316  {
317  size_type l = 0;
318  for (size_type jb = 0; jb < n_block_cols(); ++jb)
319  {
320  const SparsityPatternBase &b = block(ib, jb);
321  for (size_type j = 0; j < b.n_cols(); ++j)
322  if (b.exists(i, j))
323  out << l + j << " " << -static_cast<signed int>(i + k)
324  << std::endl;
325  l += b.n_cols();
326  }
327  }
328  k += block(ib, 0).n_rows();
329  }
330 }
331 
332 
333 
334 template <class SparsityPatternBase>
335 void
337  std::ostream &out) const
338 {
339  const unsigned int m = this->n_rows();
340  const unsigned int n = this->n_cols();
341  out
342  << "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\" viewBox=\"0 0 "
343  << n + 2 << " " << m + 2
344  << " \">\n"
345  "<style type=\"text/css\" >\n"
346  " <![CDATA[\n"
347  " rect.pixel {\n"
348  " fill: #ff0000;\n"
349  " }\n"
350  " ]]>\n"
351  " </style>\n\n"
352  " <rect width=\""
353  << n + 2 << "\" height=\"" << m + 2
354  << "\" fill=\"rgb(128, 128, 128)\"/>\n"
355  " <rect x=\"1\" y=\"1\" width=\""
356  << n + 0.1 << "\" height=\"" << m + 0.1
357  << "\" fill=\"rgb(255, 255, 255)\"/>\n\n";
358 
359  for (unsigned int block_i = 0; block_i < n_block_rows(); ++block_i)
360  for (unsigned int block_j = 0; block_j < n_block_rows(); ++block_j)
361  for (const auto &entry : block(block_i, block_j))
362  {
363  out << " <rect class=\"pixel\" x=\""
364  << column_indices.local_to_global(block_j, entry.column()) + 1
365  << "\" y=\""
366  << row_indices.local_to_global(block_i, entry.row()) + 1
367  << "\" width=\".9\" height=\".9\"/>\n";
368  }
369 
370  out << "</svg>" << std::endl;
371 }
372 
373 
374 
376  const size_type n_columns)
377  : BlockSparsityPatternBase<SparsityPattern>(n_rows, n_columns)
378 {}
379 
380 
381 void
383  const BlockIndices & rows,
384  const BlockIndices & cols,
385  const std::vector<std::vector<unsigned int>> &row_lengths)
386 {
387  AssertDimension(row_lengths.size(), cols.size());
388 
389  this->reinit(rows.size(), cols.size());
390  for (size_type j = 0; j < cols.size(); ++j)
391  for (size_type i = 0; i < rows.size(); ++i)
392  {
393  const size_type start = rows.local_to_global(i, 0);
394  const size_type length = rows.block_size(i);
395 
396  if (row_lengths[j].size() == 1)
397  block(i, j).reinit(rows.block_size(i),
398  cols.block_size(j),
399  row_lengths[j][0]);
400  else
401  {
402  Assert(row_lengths[j].begin() + start + length <=
403  row_lengths[j].end(),
404  ExcInternalError());
405  ArrayView<const unsigned int> block_rows(row_lengths[j].data() +
406  start,
407  length);
408  block(i, j).reinit(rows.block_size(i),
409  cols.block_size(j),
410  block_rows);
411  }
412  }
413  this->collect_sizes();
414  Assert(this->row_indices == rows, ExcInternalError());
415  Assert(this->column_indices == cols, ExcInternalError());
416 }
417 
418 
419 bool
421 {
422  for (size_type i = 0; i < rows; ++i)
423  for (size_type j = 0; j < columns; ++j)
424  if (sub_objects[i][j]->is_compressed() == false)
425  return false;
426  return true;
427 }
428 
429 
430 std::size_t
432 {
433  std::size_t mem = 0;
439  for (size_type r = 0; r < rows; ++r)
440  for (size_type c = 0; c < columns; ++c)
442 
443  return mem;
444 }
445 
446 
447 
448 void
450 {
451  // delete old content, set block
452  // sizes anew
453  reinit(dsp.n_block_rows(), dsp.n_block_cols());
454 
455  // copy over blocks
456  for (size_type i = 0; i < n_block_rows(); ++i)
457  for (size_type j = 0; j < n_block_cols(); ++j)
458  block(i, j).copy_from(dsp.block(i, j));
459 
460  // and finally enquire their new
461  // sizes
462  collect_sizes();
463 }
464 
465 
466 
468  const size_type n_rows,
469  const size_type n_columns)
471 {}
472 
473 
474 
476  const std::vector<size_type> &row_indices,
477  const std::vector<size_type> &col_indices)
478  : BlockSparsityPatternBase<DynamicSparsityPattern>(row_indices.size(),
479  col_indices.size())
480 {
481  for (size_type i = 0; i < row_indices.size(); ++i)
482  for (size_type j = 0; j < col_indices.size(); ++j)
483  this->block(i, j).reinit(row_indices[i], col_indices[j]);
484  this->collect_sizes();
485 }
486 
487 
488 
490  const std::vector<IndexSet> &partitioning)
491  : BlockSparsityPatternBase<DynamicSparsityPattern>(partitioning.size(),
492  partitioning.size())
493 {
494  for (size_type i = 0; i < partitioning.size(); ++i)
495  for (size_type j = 0; j < partitioning.size(); ++j)
496  this->block(i, j).reinit(partitioning[i].size(),
497  partitioning[j].size(),
498  partitioning[i]);
499  this->collect_sizes();
500 }
501 
502 
503 
505  const BlockIndices &row_indices,
506  const BlockIndices &col_indices)
507 {
508  reinit(row_indices, col_indices);
509 }
510 
511 
512 
513 void
515  const std::vector<size_type> &row_block_sizes,
516  const std::vector<size_type> &col_block_sizes)
517 {
519  row_block_sizes.size(), col_block_sizes.size());
520  for (size_type i = 0; i < row_block_sizes.size(); ++i)
521  for (size_type j = 0; j < col_block_sizes.size(); ++j)
522  this->block(i, j).reinit(row_block_sizes[i], col_block_sizes[j]);
523  this->collect_sizes();
524 }
525 
526 
527 
528 void
529 BlockDynamicSparsityPattern::reinit(const std::vector<IndexSet> &partitioning)
530 {
532  partitioning.size());
533  for (size_type i = 0; i < partitioning.size(); ++i)
534  for (size_type j = 0; j < partitioning.size(); ++j)
535  this->block(i, j).reinit(partitioning[i].size(),
536  partitioning[j].size(),
537  partitioning[i]);
538  this->collect_sizes();
539 }
540 
541 
542 
543 void
545  const BlockIndices &col_indices)
546 {
548  col_indices.size());
549  for (size_type i = 0; i < row_indices.size(); ++i)
550  for (size_type j = 0; j < col_indices.size(); ++j)
551  this->block(i, j).reinit(row_indices.block_size(i),
552  col_indices.block_size(j));
553  this->collect_sizes();
554 }
555 
556 
557 #ifdef DEAL_II_WITH_TRILINOS
558 namespace TrilinosWrappers
559 {
561  const size_type n_columns)
562  : ::BlockSparsityPatternBase<SparsityPattern>(n_rows, n_columns)
563  {}
564 
565 
566 
568  const std::vector<size_type> &row_indices,
569  const std::vector<size_type> &col_indices)
570  : BlockSparsityPatternBase<SparsityPattern>(row_indices.size(),
571  col_indices.size())
572  {
573  for (size_type i = 0; i < row_indices.size(); ++i)
574  for (size_type j = 0; j < col_indices.size(); ++j)
575  this->block(i, j).reinit(row_indices[i], col_indices[j]);
576  this->collect_sizes();
577  }
578 
579 
580 
582  const std::vector<IndexSet> &parallel_partitioning,
583  const MPI_Comm & communicator)
584  : BlockSparsityPatternBase<SparsityPattern>(parallel_partitioning.size(),
585  parallel_partitioning.size())
586  {
587  for (size_type i = 0; i < parallel_partitioning.size(); ++i)
588  for (size_type j = 0; j < parallel_partitioning.size(); ++j)
589  this->block(i, j).reinit(parallel_partitioning[i],
590  parallel_partitioning[j],
591  communicator);
592  this->collect_sizes();
593  }
594 
595 
596 
598  const std::vector<IndexSet> &row_parallel_partitioning,
599  const std::vector<IndexSet> &col_parallel_partitioning,
600  const std::vector<IndexSet> &writable_rows,
601  const MPI_Comm & communicator)
603  row_parallel_partitioning.size(),
604  col_parallel_partitioning.size())
605  {
606  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
607  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
608  this->block(i, j).reinit(row_parallel_partitioning[i],
609  col_parallel_partitioning[j],
610  writable_rows[i],
611  communicator);
612  this->collect_sizes();
613  }
614 
615 
616 
617  void
618  BlockSparsityPattern::reinit(const std::vector<size_type> &row_block_sizes,
619  const std::vector<size_type> &col_block_sizes)
620  {
622  row_block_sizes.size(), col_block_sizes.size());
623  for (size_type i = 0; i < row_block_sizes.size(); ++i)
624  for (size_type j = 0; j < col_block_sizes.size(); ++j)
625  this->block(i, j).reinit(row_block_sizes[i], col_block_sizes[j]);
626  this->collect_sizes();
627  }
628 
629 
630 
631  void
633  const std::vector<IndexSet> &parallel_partitioning,
634  const MPI_Comm & communicator)
635  {
637  parallel_partitioning.size(), parallel_partitioning.size());
638  for (size_type i = 0; i < parallel_partitioning.size(); ++i)
639  for (size_type j = 0; j < parallel_partitioning.size(); ++j)
640  this->block(i, j).reinit(parallel_partitioning[i],
641  parallel_partitioning[j],
642  communicator);
643  this->collect_sizes();
644  }
645 
646 
647 
648  void
650  const std::vector<IndexSet> &row_parallel_partitioning,
651  const std::vector<IndexSet> &col_parallel_partitioning,
652  const MPI_Comm & communicator)
653  {
655  row_parallel_partitioning.size(), col_parallel_partitioning.size());
656  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
657  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
658  this->block(i, j).reinit(row_parallel_partitioning[i],
659  col_parallel_partitioning[j],
660  communicator);
661  this->collect_sizes();
662  }
663 
664 
665 
666  void
668  const std::vector<IndexSet> &row_parallel_partitioning,
669  const std::vector<IndexSet> &col_parallel_partitioning,
670  const std::vector<IndexSet> &writable_rows,
671  const MPI_Comm & communicator)
672  {
673  AssertDimension(writable_rows.size(), row_parallel_partitioning.size());
675  row_parallel_partitioning.size(), col_parallel_partitioning.size());
676  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
677  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
678  this->block(i, j).reinit(row_parallel_partitioning[i],
679  col_parallel_partitioning[j],
680  writable_rows[i],
681  communicator);
682  this->collect_sizes();
683  }
684 
685 } // namespace TrilinosWrappers
686 
687 #endif
688 
691 #ifdef DEAL_II_WITH_TRILINOS
693 #endif
694 
static ::ExceptionBase & ExcIncompatibleRowNumbers(int arg1, int arg2, int arg3, int arg4)
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
size_type block_size(const unsigned int i) const
SparsityPatternType & block(const size_type row, const size_type column)
size_type size() const
Definition: index_set.h:1639
void reinit(const size_type n_block_rows, const size_type n_block_columns)
void reinit(const size_type n_block_rows, const size_type n_block_columns)
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)
std::size_t memory_consumption() const
void print(std::ostream &out) const
#define Assert(cond, exc)
Definition: exceptions.h:1461
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void reinit(const unsigned int n_blocks, const size_type n_elements_per_block)
const IndexSet & row_index_set() const
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
VectorType::value_type * end(VectorType &V)
BlockSparsityPatternBase & operator=(const BlockSparsityPatternBase &)
size_type local_to_global(const unsigned int block, const size_type index) const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void print_svg(std::ostream &out) const
bool exists(const size_type i, const size_type j) const
Table< 2, SmartPointer< SparsityPatternType, BlockSparsityPatternBase< SparsityPatternType > > > sub_objects
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
VectorType::value_type * begin(VectorType &V)
void print_gnuplot(std::ostream &out) const
bool is_element(const size_type index) const
Definition: index_set.h:1770
void copy_from(const BlockDynamicSparsityPattern &dsp)
bool exists(const size_type i, const size_type j) const
unsigned int size() const
BlockSparsityPattern()=default
std::enable_if< std::is_fundamental< T >::value, std::size_t >::type memory_consumption(const T &t)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
size_type n_cols() const