Reference documentation for deal.II version GIT relicensing-384-g4769914377 2024-04-13 09:00:02+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\}}\)
Loading...
Searching...
No Matches
block_sparsity_pattern.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2000 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
17
19
21
22
23template <typename SparsityPatternType>
28
29
30
31template <typename SparsityPatternType>
39
40
41
42template <typename SparsityPatternType>
46{
47 (void)s;
48 Assert(s.n_block_rows() == 0 && s.n_block_cols() == 0,
50 "This constructor can only be called if the provided argument "
51 "is the sparsity pattern for an empty matrix. This constructor can "
52 "not be used to copy-construct a non-empty sparsity pattern."));
53}
54
55
56
57template <typename SparsityPatternType>
58void
60 const size_type new_block_rows,
61 const size_type new_block_columns)
62{
63 sub_objects.reinit(0, 0);
64
65 block_rows = new_block_rows;
66 block_columns = new_block_columns;
67
68 sub_objects.reinit(block_rows, block_columns);
69 for (size_type i = 0; i < n_block_rows(); ++i)
70 for (size_type j = 0; j < n_block_cols(); ++j)
71 sub_objects[i][j] = std::make_unique<SparsityPatternType>();
72}
73
74
75template <typename SparsityPatternType>
79{
80 AssertDimension(n_block_rows(), bsp.n_block_rows());
81 AssertDimension(n_block_cols(), bsp.n_block_cols());
82 // copy objects
83 for (size_type i = 0; i < n_block_rows(); ++i)
84 for (size_type j = 0; j < n_block_cols(); ++j)
85 *sub_objects[i][j] = *bsp.sub_objects[i][j];
86 // update index objects
87 collect_sizes();
88
89 return *this;
90}
91
92
93
94template <typename SparsityPatternType>
97{
98 // only count in first column, since
99 // all rows should be equivalent
100 size_type count = 0;
101 for (size_type r = 0; r < n_block_rows(); ++r)
102 count += sub_objects[r][0]->n_rows();
103 return count;
104}
105
106
107
108template <typename SparsityPatternType>
111{
112 // only count in first row, since
113 // all rows should be equivalent
114 size_type count = 0;
115 for (size_type c = 0; c < n_block_cols(); ++c)
116 count += sub_objects[0][c]->n_cols();
117 return count;
119
120
121
122template <typename SparsityPatternType>
123void
125{
126 SparsityPatternBase::resize(compute_n_rows(), compute_n_cols());
127
128 std::vector<size_type> row_sizes(n_block_rows());
129 std::vector<size_type> col_sizes(n_block_cols());
130
131 // first find out the row sizes
132 // from the first block column
133 for (size_type r = 0; r < n_block_rows(); ++r)
134 row_sizes[r] = sub_objects[r][0]->n_rows();
135 // then check that the following
136 // block columns have the same
137 // sizes
138 for (size_type c = 1; c < n_block_cols(); ++c)
139 for (size_type r = 0; r < n_block_rows(); ++r)
140 Assert(row_sizes[r] == sub_objects[r][c]->n_rows(),
141 ExcIncompatibleRowNumbers(r, 0, r, c));
143 // finally initialize the row
144 // indices with this array
145 row_indices.reinit(row_sizes);
146
147
148 // then do the same with the columns
149 for (size_type c = 0; c < n_block_cols(); ++c)
150 col_sizes[c] = sub_objects[0][c]->n_cols();
151 for (size_type r = 1; r < n_block_rows(); ++r)
152 for (size_type c = 0; c < n_block_cols(); ++c)
153 Assert(col_sizes[c] == sub_objects[r][c]->n_cols(),
154 ExcIncompatibleRowNumbers(0, c, r, c));
155
156 // finally initialize the row
157 // indices with this array
158 column_indices.reinit(col_sizes);
159
160 // Resize scratch arrays
161 block_column_indices.resize(n_block_cols());
162 counter_within_block.resize(n_block_cols());
163}
164
165
166
167template <typename SparsityPatternType>
168void
170{
171 for (size_type i = 0; i < n_block_rows(); ++i)
172 for (size_type j = 0; j < n_block_cols(); ++j)
173 sub_objects[i][j]->compress();
174}
175
176
177
178template <typename SparsityPatternType>
179bool
181{
182 for (size_type i = 0; i < n_block_rows(); ++i)
183 for (size_type j = 0; j < n_block_cols(); ++j)
184 if (sub_objects[i][j]->empty() == false)
185 return false;
186 return true;
187}
188
189
190
191template <typename SparsityPatternType>
194{
195 size_type max_entries = 0;
196 for (size_type block_row = 0; block_row < n_block_rows(); ++block_row)
197 {
198 size_type this_row = 0;
199 for (size_type c = 0; c < n_block_cols(); ++c)
200 this_row += sub_objects[block_row][c]->max_entries_per_row();
201
202 if (this_row > max_entries)
203 max_entries = this_row;
204 }
205 return max_entries;
207
208
209
210template <typename SparsityPatternType>
213{
214 size_type count = 0;
215 for (size_type i = 0; i < n_block_rows(); ++i)
216 for (size_type j = 0; j < n_block_cols(); ++j)
217 count += sub_objects[i][j]->n_nonzero_elements();
218 return count;
219}
220
221
222
223template <typename SparsityPatternType>
224void
226{
227 size_type k = 0;
228 for (size_type ib = 0; ib < n_block_rows(); ++ib)
229 {
230 for (size_type i = 0; i < block(ib, 0).n_rows(); ++i)
231 {
232 out << '[' << i + k;
233 size_type l = 0;
234 for (size_type jb = 0; jb < n_block_cols(); ++jb)
235 {
236 const SparsityPatternType &b = block(ib, jb);
237 for (size_type j = 0; j < b.n_cols(); ++j)
238 if (b.exists(i, j))
239 out << ',' << l + j;
240 l += b.n_cols();
241 }
242 out << ']' << std::endl;
243 }
244 k += block(ib, 0).n_rows();
245 }
246}
247
248
249#ifndef DOXYGEN
250template <>
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 DynamicSparsityPattern &b = block(ib, jb);
264 if (b.row_index_set().size() == 0 ||
265 b.row_index_set().is_element(i))
266 for (size_type j = 0; j < b.n_cols(); ++j)
267 if (b.exists(i, j))
268 out << ',' << l + j;
269 l += b.n_cols();
270 }
271 out << ']' << std::endl;
272 }
273 k += block(ib, 0).n_rows();
274 }
275}
276#endif
277
278
279
280template <typename SparsityPatternType>
281void
283 std::ostream &out) const
284{
285 size_type k = 0;
286 for (size_type ib = 0; ib < n_block_rows(); ++ib)
287 {
288 for (size_type i = 0; i < block(ib, 0).n_rows(); ++i)
289 {
290 size_type l = 0;
291 for (size_type jb = 0; jb < n_block_cols(); ++jb)
292 {
293 const SparsityPatternType &b = block(ib, jb);
294 for (size_type j = 0; j < b.n_cols(); ++j)
295 if (b.exists(i, j))
296 out << l + j << " " << -static_cast<signed int>(i + k)
297 << std::endl;
298 l += b.n_cols();
299 }
300 }
301 k += block(ib, 0).n_rows();
302 }
304
305
306
307template <typename SparsityPatternType>
308void
310 std::ostream &out) const
311{
312 const unsigned int m = this->n_rows();
313 const unsigned int n = this->n_cols();
314 out
315 << "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\" viewBox=\"0 0 "
316 << n + 2 << " " << m + 2
317 << " \">\n"
318 "<style type=\"text/css\" >\n"
319 " <![CDATA[\n"
320 " rect.pixel {\n"
321 " fill: #ff0000;\n"
322 " }\n"
323 " ]]>\n"
324 " </style>\n\n"
325 " <rect width=\""
326 << n + 2 << "\" height=\"" << m + 2
327 << "\" fill=\"rgb(128, 128, 128)\"/>\n"
328 " <rect x=\"1\" y=\"1\" width=\""
329 << n + 0.1 << "\" height=\"" << m + 0.1
330 << "\" fill=\"rgb(255, 255, 255)\"/>\n\n";
331
332 for (unsigned int block_i = 0; block_i < n_block_rows(); ++block_i)
333 for (unsigned int block_j = 0; block_j < n_block_cols(); ++block_j)
334 for (const auto &entry : block(block_i, block_j))
335 {
336 out << " <rect class=\"pixel\" x=\""
337 << column_indices.local_to_global(block_j, entry.column()) + 1
338 << "\" y=\""
339 << row_indices.local_to_global(block_i, entry.row()) + 1
340 << "\" width=\".9\" height=\".9\"/>\n";
341 }
342
343 out << "</svg>" << std::endl;
344}
345
346
347
348template <typename SparsityPatternType>
349std::size_t
351{
352 std::size_t mem = 0;
353 mem += (MemoryConsumption::memory_consumption(n_block_rows()) +
358 for (size_type r = 0; r < n_block_rows(); ++r)
359 for (size_type c = 0; c < n_block_cols(); ++c)
360 mem += MemoryConsumption::memory_consumption(*sub_objects[r][c]);
361
362 return mem;
363}
364
365
366
368 const size_type n_columns)
369 : BlockSparsityPatternBase<SparsityPattern>(n_rows, n_columns)
370{}
371
372
373void
375 const BlockIndices &rows,
376 const BlockIndices &cols,
377 const std::vector<std::vector<unsigned int>> &row_lengths)
378{
379 AssertDimension(row_lengths.size(), cols.size());
380
381 this->reinit(rows.size(), cols.size());
382 for (size_type j = 0; j < cols.size(); ++j)
383 for (size_type i = 0; i < rows.size(); ++i)
384 {
385 const size_type start = rows.local_to_global(i, 0);
386 const size_type length = rows.block_size(i);
387
388 if (row_lengths[j].size() == 1)
389 block(i, j).reinit(rows.block_size(i),
390 cols.block_size(j),
391 row_lengths[j][0]);
392 else
393 {
394 Assert(row_lengths[j].begin() + start + length <=
395 row_lengths[j].end(),
397 ArrayView<const unsigned int> block_rows(row_lengths[j].data() +
398 start,
399 length);
400 block(i, j).reinit(rows.block_size(i),
401 cols.block_size(j),
402 block_rows);
403 }
404 }
406 Assert(this->row_indices == rows, ExcInternalError());
407 Assert(this->column_indices == cols, ExcInternalError());
408}
409
410
411bool
413{
414 for (size_type i = 0; i < n_block_rows(); ++i)
415 for (size_type j = 0; j < n_block_cols(); ++j)
416 if (sub_objects[i][j]->is_compressed() == false)
417 return false;
418 return true;
419}
420
421
422
423void
425{
426 // delete old content, set block
427 // sizes anew
428 reinit(dsp.n_block_rows(), dsp.n_block_cols());
429
430 // copy over blocks
431 for (size_type i = 0; i < n_block_rows(); ++i)
432 for (size_type j = 0; j < n_block_cols(); ++j)
433 block(i, j).copy_from(dsp.block(i, j));
434
435 // and finally enquire their new
436 // sizes
438}
439
440
441
447
448
449
451 const std::vector<size_type> &row_indices,
452 const std::vector<size_type> &col_indices)
454 col_indices.size())
455{
456 for (size_type i = 0; i < row_indices.size(); ++i)
457 for (size_type j = 0; j < col_indices.size(); ++j)
458 this->block(i, j).reinit(row_indices[i], col_indices[j]);
459 this->collect_sizes();
460}
461
462
463
465 const std::vector<IndexSet> &partitioning)
467 partitioning.size())
468{
469 for (size_type i = 0; i < partitioning.size(); ++i)
470 for (size_type j = 0; j < partitioning.size(); ++j)
471 this->block(i, j).reinit(partitioning[i].size(),
472 partitioning[j].size(),
473 partitioning[i]);
474 this->collect_sizes();
475}
476
477
478
480 const BlockIndices &row_indices,
481 const BlockIndices &col_indices)
482{
483 reinit(row_indices, col_indices);
484}
485
486
487
488void
490 const std::vector<size_type> &row_block_sizes,
491 const std::vector<size_type> &col_block_sizes)
492{
494 row_block_sizes.size(), col_block_sizes.size());
495 for (size_type i = 0; i < row_block_sizes.size(); ++i)
496 for (size_type j = 0; j < col_block_sizes.size(); ++j)
497 this->block(i, j).reinit(row_block_sizes[i], col_block_sizes[j]);
498 this->collect_sizes();
499}
500
501
502
503void
504BlockDynamicSparsityPattern::reinit(const std::vector<IndexSet> &partitioning)
505{
507 partitioning.size());
508 for (size_type i = 0; i < partitioning.size(); ++i)
509 for (size_type j = 0; j < partitioning.size(); ++j)
510 this->block(i, j).reinit(partitioning[i].size(),
511 partitioning[j].size(),
512 partitioning[i]);
513 this->collect_sizes();
514}
515
516
517
518void
520 const BlockIndices &col_indices)
521{
523 col_indices.size());
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.block_size(i),
527 col_indices.block_size(j));
528 this->collect_sizes();
529}
530
531
532#ifdef DEAL_II_WITH_TRILINOS
533namespace TrilinosWrappers
534{
536 const size_type n_columns)
537 : ::BlockSparsityPatternBase<SparsityPattern>(n_rows, n_columns)
538 {}
539
540
541
543 const std::vector<size_type> &row_indices,
544 const std::vector<size_type> &col_indices)
545 : BlockSparsityPatternBase<SparsityPattern>(row_indices.size(),
546 col_indices.size())
547 {
548 for (size_type i = 0; i < row_indices.size(); ++i)
549 for (size_type j = 0; j < col_indices.size(); ++j)
550 this->block(i, j).reinit(row_indices[i], col_indices[j]);
551 this->collect_sizes();
552 }
553
554
555
557 const std::vector<IndexSet> &parallel_partitioning,
558 const MPI_Comm communicator)
559 : BlockSparsityPatternBase<SparsityPattern>(parallel_partitioning.size(),
560 parallel_partitioning.size())
561 {
562 for (size_type i = 0; i < parallel_partitioning.size(); ++i)
563 for (size_type j = 0; j < parallel_partitioning.size(); ++j)
564 this->block(i, j).reinit(parallel_partitioning[i],
565 parallel_partitioning[j],
566 communicator);
567 this->collect_sizes();
568 }
569
570
571
573 const std::vector<IndexSet> &row_parallel_partitioning,
574 const std::vector<IndexSet> &col_parallel_partitioning,
575 const std::vector<IndexSet> &writable_rows,
576 const MPI_Comm communicator)
578 row_parallel_partitioning.size(),
579 col_parallel_partitioning.size())
580 {
581 for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
582 for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
583 this->block(i, j).reinit(row_parallel_partitioning[i],
584 col_parallel_partitioning[j],
585 writable_rows[i],
586 communicator);
587 this->collect_sizes();
588 }
589
590
591
592 void
593 BlockSparsityPattern::reinit(const std::vector<size_type> &row_block_sizes,
594 const std::vector<size_type> &col_block_sizes)
595 {
597 row_block_sizes.size(), col_block_sizes.size());
598 for (size_type i = 0; i < row_block_sizes.size(); ++i)
599 for (size_type j = 0; j < col_block_sizes.size(); ++j)
600 this->block(i, j).reinit(row_block_sizes[i], col_block_sizes[j]);
601 this->collect_sizes();
602 }
603
604
605
606 void
608 const std::vector<IndexSet> &parallel_partitioning,
609 const MPI_Comm communicator)
610 {
612 parallel_partitioning.size(), parallel_partitioning.size());
613 for (size_type i = 0; i < parallel_partitioning.size(); ++i)
614 for (size_type j = 0; j < parallel_partitioning.size(); ++j)
615 this->block(i, j).reinit(parallel_partitioning[i],
616 parallel_partitioning[j],
617 communicator);
618 this->collect_sizes();
619 }
620
621
622
623 void
625 const std::vector<IndexSet> &row_parallel_partitioning,
626 const std::vector<IndexSet> &col_parallel_partitioning,
627 const MPI_Comm communicator)
628 {
630 row_parallel_partitioning.size(), col_parallel_partitioning.size());
631 for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
632 for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
633 this->block(i, j).reinit(row_parallel_partitioning[i],
634 col_parallel_partitioning[j],
635 communicator);
636 this->collect_sizes();
637 }
638
639
640
641 void
643 const std::vector<IndexSet> &row_parallel_partitioning,
644 const std::vector<IndexSet> &col_parallel_partitioning,
645 const std::vector<IndexSet> &writable_rows,
646 const MPI_Comm communicator)
647 {
648 AssertDimension(writable_rows.size(), row_parallel_partitioning.size());
650 row_parallel_partitioning.size(), col_parallel_partitioning.size());
651 for (size_type i = 0; i < row_parallel_partitioning.size(); ++i)
652 for (size_type j = 0; j < col_parallel_partitioning.size(); ++j)
653 this->block(i, j).reinit(row_parallel_partitioning[i],
654 col_parallel_partitioning[j],
655 writable_rows[i],
656 communicator);
657 this->collect_sizes();
658 }
659
660} // namespace TrilinosWrappers
661
662#endif
663
666#ifdef DEAL_II_WITH_TRILINOS
668#endif
669
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::enable_if_t< std::is_fundamental_v< T >, std::size_t > memory_consumption(const T &t)
const InputIterator OutputIterator out
Definition parallel.h:167