Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
24template <class SparsityPatternType>
26 : block_rows(0)
27 , block_columns(0)
28{}
29
30
31
32template <class SparsityPatternType>
34 const size_type block_rows,
35 const size_type block_columns)
37{
39}
40
41
42
43template <class SparsityPatternType>
47{
48 (void)s;
49 Assert(s.n_block_rows() == 0 && s.n_block_cols() == 0,
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
58template <class SparsityPatternType>
59void
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
76template <class 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
95template <class 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
109template <class 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;
120
121
122
123template <class SparsityPatternType>
124void
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));
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
168template <class SparsityPatternType>
169void
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
179template <class SparsityPatternType>
180bool
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
192template <class 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;
208
209
210
211template <class 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
224template <class SparsityPatternType>
225void
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
251template <>
252void
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
281template <class SparsityPatternType>
282void
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 }
305
306
307
308template <class SparsityPatternType>
309void
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
349template <class SparsityPatternType>
350std::size_t
352{
353 std::size_t mem = 0;
354 mem += (MemoryConsumption::memory_consumption(n_block_rows()) +
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
374void
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(),
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 }
407 Assert(this->row_indices == rows, ExcInternalError());
408 Assert(this->column_indices == cols, ExcInternalError());
409}
410
411
412bool
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
424void
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
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)
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)
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
489void
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
504void
505BlockDynamicSparsityPattern::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
519void
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
534namespace 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
SparsityPattern & block(const size_type row, const size_type column)
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)
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:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::enable_if_t< std::is_fundamental< T >::value, std::size_t > memory_consumption(const T &t)