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