Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
block_sparsity_pattern.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2000 - 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/base/memory_consumption.h>
18 
19 #include <deal.II/lac/block_sparsity_pattern.h>
20 
21 DEAL_II_NAMESPACE_OPEN
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 template <>
277 void
279 {
280  size_type k = 0;
281  for (size_type ib = 0; ib < n_block_rows(); ++ib)
282  {
283  for (size_type i = 0; i < block(ib, 0).n_rows(); ++i)
284  {
285  out << '[' << i + k;
286  size_type l = 0;
287  for (size_type jb = 0; jb < n_block_cols(); ++jb)
288  {
289  const DynamicSparsityPattern &b = block(ib, jb);
290  if (b.row_index_set().size() == 0 ||
291  b.row_index_set().is_element(i))
292  for (size_type j = 0; j < b.n_cols(); ++j)
293  if (b.exists(i, j))
294  out << ',' << l + j;
295  l += b.n_cols();
296  }
297  out << ']' << std::endl;
298  }
299  k += block(ib, 0).n_rows();
300  }
301 }
302 
303 
304 template <class SparsityPatternBase>
305 void
307  std::ostream &out) const
308 {
309  size_type k = 0;
310  for (size_type ib = 0; ib < n_block_rows(); ++ib)
311  {
312  for (size_type i = 0; i < block(ib, 0).n_rows(); ++i)
313  {
314  size_type l = 0;
315  for (size_type jb = 0; jb < n_block_cols(); ++jb)
316  {
317  const SparsityPatternBase &b = block(ib, jb);
318  for (size_type j = 0; j < b.n_cols(); ++j)
319  if (b.exists(i, j))
320  out << l + j << " " << -static_cast<signed int>(i + k)
321  << std::endl;
322  l += b.n_cols();
323  }
324  }
325  k += block(ib, 0).n_rows();
326  }
327 }
328 
329 
330 
332  const size_type n_columns)
333  : BlockSparsityPatternBase<SparsityPattern>(n_rows, n_columns)
334 {}
335 
336 
337 void
339  const BlockIndices & rows,
340  const BlockIndices & cols,
341  const std::vector<std::vector<unsigned int>> &row_lengths)
342 {
343  AssertDimension(row_lengths.size(), cols.size());
344 
345  this->reinit(rows.size(), cols.size());
346  for (size_type j = 0; j < cols.size(); ++j)
347  for (size_type i = 0; i < rows.size(); ++i)
348  {
349  const size_type start = rows.local_to_global(i, 0);
350  const size_type length = rows.block_size(i);
351 
352  if (row_lengths[j].size() == 1)
353  block(i, j).reinit(rows.block_size(i),
354  cols.block_size(j),
355  row_lengths[j][0]);
356  else
357  {
358  Assert(row_lengths[j].begin() + start + length <=
359  row_lengths[j].end(),
360  ExcInternalError());
361  ArrayView<const unsigned int> block_rows(row_lengths[j].data() +
362  start,
363  length);
364  block(i, j).reinit(rows.block_size(i),
365  cols.block_size(j),
366  block_rows);
367  }
368  }
369  this->collect_sizes();
370  Assert(this->row_indices == rows, ExcInternalError());
371  Assert(this->column_indices == cols, ExcInternalError());
372 }
373 
374 
375 bool
377 {
378  for (size_type i = 0; i < rows; ++i)
379  for (size_type j = 0; j < columns; ++j)
380  if (sub_objects[i][j]->is_compressed() == false)
381  return false;
382  return true;
383 }
384 
385 
386 std::size_t
388 {
389  std::size_t mem = 0;
395  for (size_type r = 0; r < rows; ++r)
396  for (size_type c = 0; c < columns; ++c)
398 
399  return mem;
400 }
401 
402 
403 
404 void
406 {
407  // delete old content, set block
408  // sizes anew
409  reinit(dsp.n_block_rows(), dsp.n_block_cols());
410 
411  // copy over blocks
412  for (size_type i = 0; i < n_block_rows(); ++i)
413  for (size_type j = 0; j < n_block_cols(); ++j)
414  block(i, j).copy_from(dsp.block(i, j));
415 
416  // and finally enquire their new
417  // sizes
418  collect_sizes();
419 }
420 
421 
422 
424  const size_type n_rows,
425  const size_type n_columns)
427 {}
428 
429 
430 
432  const std::vector<size_type> &row_indices,
433  const std::vector<size_type> &col_indices)
434  : BlockSparsityPatternBase<DynamicSparsityPattern>(row_indices.size(),
435  col_indices.size())
436 {
437  for (size_type i = 0; i < row_indices.size(); ++i)
438  for (size_type j = 0; j < col_indices.size(); ++j)
439  this->block(i, j).reinit(row_indices[i], col_indices[j]);
440  this->collect_sizes();
441 }
442 
443 
445  const std::vector<IndexSet> &partitioning)
446  : BlockSparsityPatternBase<DynamicSparsityPattern>(partitioning.size(),
447  partitioning.size())
448 {
449  for (size_type i = 0; i < partitioning.size(); ++i)
450  for (size_type j = 0; j < partitioning.size(); ++j)
451  this->block(i, j).reinit(partitioning[i].size(),
452  partitioning[j].size(),
453  partitioning[i]);
454  this->collect_sizes();
455 }
456 
457 
459  const BlockIndices &row_indices,
460  const BlockIndices &col_indices)
461 {
462  reinit(row_indices, col_indices);
463 }
464 
465 
466 void
468  const std::vector<size_type> &row_block_sizes,
469  const std::vector<size_type> &col_block_sizes)
470 {
472  row_block_sizes.size(), col_block_sizes.size());
473  for (size_type i = 0; i < row_block_sizes.size(); ++i)
474  for (size_type j = 0; j < col_block_sizes.size(); ++j)
475  this->block(i, j).reinit(row_block_sizes[i], col_block_sizes[j]);
476  this->collect_sizes();
477 }
478 
479 void
480 BlockDynamicSparsityPattern::reinit(const std::vector<IndexSet> &partitioning)
481 {
483  partitioning.size());
484  for (size_type i = 0; i < partitioning.size(); ++i)
485  for (size_type j = 0; j < partitioning.size(); ++j)
486  this->block(i, j).reinit(partitioning[i].size(),
487  partitioning[j].size(),
488  partitioning[i]);
489  this->collect_sizes();
490 }
491 
492 void
494  const BlockIndices &col_indices)
495 {
497  col_indices.size());
498  for (size_type i = 0; i < row_indices.size(); ++i)
499  for (size_type j = 0; j < col_indices.size(); ++j)
500  this->block(i, j).reinit(row_indices.block_size(i),
501  col_indices.block_size(j));
502  this->collect_sizes();
503 }
504 
505 
506 #ifdef DEAL_II_WITH_TRILINOS
507 namespace TrilinosWrappers
508 {
510  const size_type n_columns)
511  : ::BlockSparsityPatternBase<SparsityPattern>(n_rows, n_columns)
512  {}
513 
514 
515 
517  const std::vector<size_type> &row_indices,
518  const std::vector<size_type> &col_indices)
519  : BlockSparsityPatternBase<SparsityPattern>(row_indices.size(),
520  col_indices.size())
521  {
522  for (size_type i = 0; i < row_indices.size(); ++i)
523  for (size_type j = 0; j < col_indices.size(); ++j)
524  this->block(i, j).reinit(row_indices[i], col_indices[j]);
525  this->collect_sizes();
526  }
527 
528 
529 
531  const std::vector<Epetra_Map> &parallel_partitioning)
532  : BlockSparsityPatternBase<SparsityPattern>(parallel_partitioning.size(),
533  parallel_partitioning.size())
534  {
535  for (size_type i = 0; i < parallel_partitioning.size(); ++i)
536  for (size_type j = 0; j < parallel_partitioning.size(); ++j)
537  this->block(i, j).reinit(parallel_partitioning[i],
538  parallel_partitioning[j]);
539  this->collect_sizes();
540  }
541 
542 
543 
545  const std::vector<IndexSet> &parallel_partitioning,
546  const MPI_Comm & communicator)
547  : BlockSparsityPatternBase<SparsityPattern>(parallel_partitioning.size(),
548  parallel_partitioning.size())
549  {
550  for (size_type i = 0; i < parallel_partitioning.size(); ++i)
551  for (size_type j = 0; j < parallel_partitioning.size(); ++j)
552  this->block(i, j).reinit(parallel_partitioning[i],
553  parallel_partitioning[j],
554  communicator);
555  this->collect_sizes();
556  }
557 
558 
559 
561  const std::vector<IndexSet> &row_parallel_partitioning,
562  const std::vector<IndexSet> &col_parallel_partitioning,
563  const std::vector<IndexSet> &writable_rows,
564  const MPI_Comm & communicator)
566  row_parallel_partitioning.size(),
567  col_parallel_partitioning.size())
568  {
569  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
570  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
571  this->block(i, j).reinit(row_parallel_partitioning[i],
572  col_parallel_partitioning[j],
573  writable_rows[i],
574  communicator);
575  this->collect_sizes();
576  }
577 
578 
579 
580  void
581  BlockSparsityPattern::reinit(const std::vector<size_type> &row_block_sizes,
582  const std::vector<size_type> &col_block_sizes)
583  {
585  row_block_sizes.size(), col_block_sizes.size());
586  for (size_type i = 0; i < row_block_sizes.size(); ++i)
587  for (size_type j = 0; j < col_block_sizes.size(); ++j)
588  this->block(i, j).reinit(row_block_sizes[i], col_block_sizes[j]);
589  this->collect_sizes();
590  }
591 
592 
593 
594  void
596  const std::vector<Epetra_Map> &parallel_partitioning)
597  {
599  parallel_partitioning.size(), parallel_partitioning.size());
600  for (size_type i = 0; i < parallel_partitioning.size(); ++i)
601  for (size_type j = 0; j < parallel_partitioning.size(); ++j)
602  this->block(i, j).reinit(parallel_partitioning[i],
603  parallel_partitioning[j]);
604  this->collect_sizes();
605  }
606 
607 
608 
609  void
611  const std::vector<IndexSet> &parallel_partitioning,
612  const MPI_Comm & communicator)
613  {
615  parallel_partitioning.size(), parallel_partitioning.size());
616  for (size_type i = 0; i < parallel_partitioning.size(); ++i)
617  for (size_type j = 0; j < parallel_partitioning.size(); ++j)
618  this->block(i, j).reinit(parallel_partitioning[i],
619  parallel_partitioning[j],
620  communicator);
621  this->collect_sizes();
622  }
623 
624 
625 
626  void
628  const std::vector<IndexSet> &row_parallel_partitioning,
629  const std::vector<IndexSet> &col_parallel_partitioning,
630  const MPI_Comm & communicator)
631  {
633  row_parallel_partitioning.size(), col_parallel_partitioning.size());
634  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
635  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
636  this->block(i, j).reinit(row_parallel_partitioning[i],
637  col_parallel_partitioning[j],
638  communicator);
639  this->collect_sizes();
640  }
641 
642 
643 
644  void
646  const std::vector<IndexSet> &row_parallel_partitioning,
647  const std::vector<IndexSet> &col_parallel_partitioning,
648  const std::vector<IndexSet> &writable_rows,
649  const MPI_Comm & communicator)
650  {
651  AssertDimension(writable_rows.size(), row_parallel_partitioning.size());
653  row_parallel_partitioning.size(), col_parallel_partitioning.size());
654  for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
655  for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
656  this->block(i, j).reinit(row_parallel_partitioning[i],
657  col_parallel_partitioning[j],
658  writable_rows[i],
659  communicator);
660  this->collect_sizes();
661  }
662 
663 } // namespace TrilinosWrappers
664 
665 #endif
666 
669 #ifdef DEAL_II_WITH_TRILINOS
671 #endif
672 
673 DEAL_II_NAMESPACE_CLOSE
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:1567
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:1600
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:1407
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
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)
bool exists(const size_type i, const size_type j) const
Table< 2, SmartPointer< SparsityPatternType, BlockSparsityPatternBase< SparsityPatternType > > > sub_objects
void print_gnuplot(std::ostream &out) const
bool is_element(const size_type index) const
Definition: index_set.h:1732
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)
static ::ExceptionBase & ExcInternalError()
size_type n_cols() const