Reference documentation for deal.II version 9.3.3
\(\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
24template <class SparsityPatternBase>
26 : rows(0)
27 , columns(0)
28{}
29
30
31
32template <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
44template <class SparsityPatternBase>
47 : Subscriptor()
48 , rows(0)
49 , columns(0)
50{
51 (void)s;
52 Assert(s.rows == 0 && s.columns == 0,
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
61template <class SparsityPatternBase>
63{
64 // clear all memory
65 try
66 {
67 reinit(0, 0);
68 }
69 catch (...)
70 {}
71}
72
73
74
75template <class SparsityPatternBase>
76void
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)
103 sub_objects[i][j] = p;
104 }
105}
106
107
108template <class SparsityPatternBase>
112{
113 Assert(rows == bsp.rows, ExcDimensionMismatch(rows, bsp.rows));
114 Assert(columns == bsp.columns, ExcDimensionMismatch(columns, bsp.columns));
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
121
122 return *this;
123}
124
125
126
127template <class SparsityPatternBase>
128void
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
166template <class SparsityPatternBase>
167void
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
177template <class SparsityPatternBase>
178bool
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
190template <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
209template <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;
220
221
222
223template <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
237template <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
250template <class SparsityPatternBase>
251void
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
277template <>
278void
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 {
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
307template <class SparsityPatternBase>
308void
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)
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
334template <class SparsityPatternBase>
335void
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
381void
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(),
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
419bool
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
430std::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
448void
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
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)
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)
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
513void
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
528void
529BlockDynamicSparsityPattern::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
543void
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
558namespace 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
size_type block_size(const unsigned int i) const
unsigned int size() const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
#define Assert(cond, exc)
Definition: exceptions.h:1465
void copy_from(const BlockDynamicSparsityPattern &dsp)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
Table< 2, SmartPointer< SparsityPatternType, BlockSparsityPatternBase< SparsityPatternType > > > sub_objects
void reinit(const size_type n_block_rows, const size_type n_block_columns)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
BlockSparsityPattern()=default
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_gnuplot(std::ostream &out) const
void print(std::ostream &out) const
SparsityPatternType & block(const size_type row, const size_type column)
void reinit(const size_type m, const size_type n, const unsigned int max_per_row)
types::global_dof_index size_type
void reinit(const size_type m, const size_type n, const IndexSet &rowset=IndexSet())
BlockSparsityPatternBase & operator=(const BlockSparsityPatternBase &)
void print_svg(std::ostream &out) const
void reinit(const size_type n_block_rows, const size_type n_block_columns)
void reinit(const std::vector< size_type > &row_block_sizes, const std::vector< size_type > &col_block_sizes)
std::enable_if< IsBlockVector< VectorType >::value, void >::type collect_sizes(VectorType &vector)
Definition: operators.h:97
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)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
std::string compress(const std::string &input)
Definition: utilities.cc:392
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:618