Reference documentation for deal.II version Git bed997f895 2020-09-22 11:49:20 -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 template <class SparsityPatternBase>
307 void
309  std::ostream &out) const
310 {
311  size_type k = 0;
312  for (size_type ib = 0; ib < n_block_rows(); ++ib)
313  {
314  for (size_type i = 0; i < block(ib, 0).n_rows(); ++i)
315  {
316  size_type l = 0;
317  for (size_type jb = 0; jb < n_block_cols(); ++jb)
318  {
319  const SparsityPatternBase &b = block(ib, jb);
320  for (size_type j = 0; j < b.n_cols(); ++j)
321  if (b.exists(i, j))
322  out << l + j << " " << -static_cast<signed int>(i + k)
323  << std::endl;
324  l += b.n_cols();
325  }
326  }
327  k += block(ib, 0).n_rows();
328  }
329 }
330 
331 
332 
334  const size_type n_columns)
335  : BlockSparsityPatternBase<SparsityPattern>(n_rows, n_columns)
336 {}
337 
338 
339 void
341  const BlockIndices & rows,
342  const BlockIndices & cols,
343  const std::vector<std::vector<unsigned int>> &row_lengths)
344 {
345  AssertDimension(row_lengths.size(), cols.size());
346 
347  this->reinit(rows.size(), cols.size());
348  for (size_type j = 0; j < cols.size(); ++j)
349  for (size_type i = 0; i < rows.size(); ++i)
350  {
351  const size_type start = rows.local_to_global(i, 0);
352  const size_type length = rows.block_size(i);
353 
354  if (row_lengths[j].size() == 1)
355  block(i, j).reinit(rows.block_size(i),
356  cols.block_size(j),
357  row_lengths[j][0]);
358  else
359  {
360  Assert(row_lengths[j].begin() + start + length <=
361  row_lengths[j].end(),
362  ExcInternalError());
363  ArrayView<const unsigned int> block_rows(row_lengths[j].data() +
364  start,
365  length);
366  block(i, j).reinit(rows.block_size(i),
367  cols.block_size(j),
368  block_rows);
369  }
370  }
371  this->collect_sizes();
372  Assert(this->row_indices == rows, ExcInternalError());
373  Assert(this->column_indices == cols, ExcInternalError());
374 }
375 
376 
377 bool
379 {
380  for (size_type i = 0; i < rows; ++i)
381  for (size_type j = 0; j < columns; ++j)
382  if (sub_objects[i][j]->is_compressed() == false)
383  return false;
384  return true;
385 }
386 
387 
388 std::size_t
390 {
391  std::size_t mem = 0;
397  for (size_type r = 0; r < rows; ++r)
398  for (size_type c = 0; c < columns; ++c)
400 
401  return mem;
402 }
403 
404 
405 
406 void
408 {
409  // delete old content, set block
410  // sizes anew
411  reinit(dsp.n_block_rows(), dsp.n_block_cols());
412 
413  // copy over blocks
414  for (size_type i = 0; i < n_block_rows(); ++i)
415  for (size_type j = 0; j < n_block_cols(); ++j)
416  block(i, j).copy_from(dsp.block(i, j));
417 
418  // and finally enquire their new
419  // sizes
420  collect_sizes();
421 }
422 
423 
424 
426  const size_type n_rows,
427  const size_type n_columns)
429 {}
430 
431 
432 
434  const std::vector<size_type> &row_indices,
435  const std::vector<size_type> &col_indices)
436  : BlockSparsityPatternBase<DynamicSparsityPattern>(row_indices.size(),
437  col_indices.size())
438 {
439  for (size_type i = 0; i < row_indices.size(); ++i)
440  for (size_type j = 0; j < col_indices.size(); ++j)
441  this->block(i, j).reinit(row_indices[i], col_indices[j]);
442  this->collect_sizes();
443 }
444 
445 
447  const std::vector<IndexSet> &partitioning)
448  : BlockSparsityPatternBase<DynamicSparsityPattern>(partitioning.size(),
449  partitioning.size())
450 {
451  for (size_type i = 0; i < partitioning.size(); ++i)
452  for (size_type j = 0; j < partitioning.size(); ++j)
453  this->block(i, j).reinit(partitioning[i].size(),
454  partitioning[j].size(),
455  partitioning[i]);
456  this->collect_sizes();
457 }
458 
459 
461  const BlockIndices &row_indices,
462  const BlockIndices &col_indices)
463 {
464  reinit(row_indices, col_indices);
465 }
466 
467 
468 void
470  const std::vector<size_type> &row_block_sizes,
471  const std::vector<size_type> &col_block_sizes)
472 {
474  row_block_sizes.size(), col_block_sizes.size());
475  for (size_type i = 0; i < row_block_sizes.size(); ++i)
476  for (size_type j = 0; j < col_block_sizes.size(); ++j)
477  this->block(i, j).reinit(row_block_sizes[i], col_block_sizes[j]);
478  this->collect_sizes();
479 }
480 
481 void
482 BlockDynamicSparsityPattern::reinit(const std::vector<IndexSet> &partitioning)
483 {
485  partitioning.size());
486  for (size_type i = 0; i < partitioning.size(); ++i)
487  for (size_type j = 0; j < partitioning.size(); ++j)
488  this->block(i, j).reinit(partitioning[i].size(),
489  partitioning[j].size(),
490  partitioning[i]);
491  this->collect_sizes();
492 }
493 
494 void
496  const BlockIndices &col_indices)
497 {
499  col_indices.size());
500  for (size_type i = 0; i < row_indices.size(); ++i)
501  for (size_type j = 0; j < col_indices.size(); ++j)
502  this->block(i, j).reinit(row_indices.block_size(i),
503  col_indices.block_size(j));
504  this->collect_sizes();
505 }
506 
507 
508 #ifdef DEAL_II_WITH_TRILINOS
509 namespace TrilinosWrappers
510 {
512  const size_type n_columns)
513  : ::BlockSparsityPatternBase<SparsityPattern>(n_rows, n_columns)
514  {}
515 
516 
517 
519  const std::vector<size_type> &row_indices,
520  const std::vector<size_type> &col_indices)
521  : BlockSparsityPatternBase<SparsityPattern>(row_indices.size(),
522  col_indices.size())
523  {
524  for (size_type i = 0; i < row_indices.size(); ++i)
525  for (size_type j = 0; j < col_indices.size(); ++j)
526  this->block(i, j).reinit(row_indices[i], col_indices[j]);
527  this->collect_sizes();
528  }
529 
530 
531 
533  const std::vector<IndexSet> &parallel_partitioning,
534  const MPI_Comm & communicator)
535  : BlockSparsityPatternBase<SparsityPattern>(parallel_partitioning.size(),
536  parallel_partitioning.size())
537  {
538  for (size_type i = 0; i < parallel_partitioning.size(); ++i)
539  for (size_type j = 0; j < parallel_partitioning.size(); ++j)
540  this->block(i, j).reinit(parallel_partitioning[i],
541  parallel_partitioning[j],
542  communicator);
543  this->collect_sizes();
544  }
545 
546 
547 
549  const std::vector<IndexSet> &row_parallel_partitioning,
550  const std::vector<IndexSet> &col_parallel_partitioning,
551  const std::vector<IndexSet> &writable_rows,
552  const MPI_Comm & communicator)
554  row_parallel_partitioning.size(),
555  col_parallel_partitioning.size())
556  {
557  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
558  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
559  this->block(i, j).reinit(row_parallel_partitioning[i],
560  col_parallel_partitioning[j],
561  writable_rows[i],
562  communicator);
563  this->collect_sizes();
564  }
565 
566 
567 
568  void
569  BlockSparsityPattern::reinit(const std::vector<size_type> &row_block_sizes,
570  const std::vector<size_type> &col_block_sizes)
571  {
573  row_block_sizes.size(), col_block_sizes.size());
574  for (size_type i = 0; i < row_block_sizes.size(); ++i)
575  for (size_type j = 0; j < col_block_sizes.size(); ++j)
576  this->block(i, j).reinit(row_block_sizes[i], col_block_sizes[j]);
577  this->collect_sizes();
578  }
579 
580 
581 
582  void
584  const std::vector<IndexSet> &parallel_partitioning,
585  const MPI_Comm & communicator)
586  {
588  parallel_partitioning.size(), parallel_partitioning.size());
589  for (size_type i = 0; i < parallel_partitioning.size(); ++i)
590  for (size_type j = 0; j < parallel_partitioning.size(); ++j)
591  this->block(i, j).reinit(parallel_partitioning[i],
592  parallel_partitioning[j],
593  communicator);
594  this->collect_sizes();
595  }
596 
597 
598 
599  void
601  const std::vector<IndexSet> &row_parallel_partitioning,
602  const std::vector<IndexSet> &col_parallel_partitioning,
603  const MPI_Comm & communicator)
604  {
606  row_parallel_partitioning.size(), col_parallel_partitioning.size());
607  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
608  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
609  this->block(i, j).reinit(row_parallel_partitioning[i],
610  col_parallel_partitioning[j],
611  communicator);
612  this->collect_sizes();
613  }
614 
615 
616 
617  void
619  const std::vector<IndexSet> &row_parallel_partitioning,
620  const std::vector<IndexSet> &col_parallel_partitioning,
621  const std::vector<IndexSet> &writable_rows,
622  const MPI_Comm & communicator)
623  {
624  AssertDimension(writable_rows.size(), row_parallel_partitioning.size());
626  row_parallel_partitioning.size(), col_parallel_partitioning.size());
627  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
628  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
629  this->block(i, j).reinit(row_parallel_partitioning[i],
630  col_parallel_partitioning[j],
631  writable_rows[i],
632  communicator);
633  this->collect_sizes();
634  }
635 
636 } // namespace TrilinosWrappers
637 
638 #endif
639 
642 #ifdef DEAL_II_WITH_TRILINOS
644 #endif
645 
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:1568
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:1632
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
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1411
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:363
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)
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:362
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:1763
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