Reference documentation for deal.II version 9.6.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
chunk_sparsity_pattern.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2008 - 2022 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
19
20
22
23
28
29
30
32 : Subscriptor()
33 , chunk_size(s.chunk_size)
34 , sparsity_pattern(s.sparsity_pattern)
35{
36 Assert(s.rows == 0 && s.cols == 0,
38 "This constructor can only be called if the provided argument "
39 "is the sparsity pattern for an empty matrix. This constructor can "
40 "not be used to copy-construct a non-empty sparsity pattern."));
41
42 reinit(0, 0, 0, chunk_size);
43}
44
45
46
48 const size_type n,
49 const size_type max_per_row,
50 const size_type chunk_size)
51{
53
54 reinit(m, n, max_per_row, chunk_size);
55}
56
57
58
60 const size_type m,
61 const size_type n,
62 const std::vector<size_type> &row_lengths,
63 const size_type chunk_size)
64{
66
67 reinit(m, n, row_lengths, chunk_size);
68}
69
70
71
73 const size_type max_per_row,
74 const size_type chunk_size)
75{
76 reinit(n, n, max_per_row, chunk_size);
77}
78
79
80
82 const size_type m,
83 const std::vector<size_type> &row_lengths,
84 const size_type chunk_size)
85{
87
88 reinit(m, m, row_lengths, chunk_size);
89}
90
91
92
95{
96 Assert(s.rows == 0 && s.cols == 0,
98 "This operator can only be called if the provided argument "
99 "is the sparsity pattern for an empty matrix. This operator can "
100 "not be used to copy a non-empty sparsity pattern."));
101
102 Assert(rows == 0 && cols == 0,
103 ExcMessage("This operator can only be called if the current object is "
104 "empty."));
105
106 // perform the checks in the underlying object as well
108
109 return *this;
110}
111
112
113
114void
116 const size_type n,
117 const size_type max_per_row,
118 const size_type chunk_size)
119{
121
122 // simply map this function to the other @p{reinit} function
123 const std::vector<size_type> row_lengths(m, max_per_row);
124 reinit(m, n, row_lengths, chunk_size);
125}
126
127
128
129void
131 const size_type n,
132 const ArrayView<const size_type> &row_lengths,
133 const size_type chunk_size)
134{
135 Assert(row_lengths.size() == m, ExcInvalidNumber(m));
137
138 rows = m;
139 cols = n;
140
141 this->chunk_size = chunk_size;
142
143 // pass down to the necessary information to the underlying object. we need
144 // to calculate how many chunks we need: we need to round up (m/chunk_size)
145 // and (n/chunk_size). rounding up in integer arithmetic equals
146 // ((m+chunk_size-1)/chunk_size):
147 const size_type m_chunks = (m + chunk_size - 1) / chunk_size,
148 n_chunks = (n + chunk_size - 1) / chunk_size;
149
150 // compute the maximum number of chunks in each row. the passed array
151 // denotes the number of entries in each row of the big matrix -- in the
152 // worst case, these are all in independent chunks, so we have to calculate
153 // it as follows (as an example: let chunk_size==2, row_lengths={2,2,...},
154 // and entries in row zero at columns {0,2} and for row one at {4,6} -->
155 // we'll need 4 chunks for the first chunk row!) :
156 std::vector<unsigned int> chunk_row_lengths(m_chunks, 0);
157 for (size_type i = 0; i < m; ++i)
158 chunk_row_lengths[i / chunk_size] += row_lengths[i];
159
160 // for the case that the reduced sparsity pattern optimizes the diagonal but
161 // the actual sparsity pattern does not, need to take one more entry in the
162 // row to fit the user-required entry
163 if (m != n && m_chunks == n_chunks)
164 for (unsigned int i = 0; i < m_chunks; ++i)
165 ++chunk_row_lengths[i];
166
167 sparsity_pattern.reinit(m_chunks, n_chunks, chunk_row_lengths);
168}
169
170
171
172void
177
178
179
180template <typename SparsityPatternType>
181void
182ChunkSparsityPattern::copy_from(const SparsityPatternType &dsp,
183 const size_type chunk_size)
184{
186 this->chunk_size = chunk_size;
187 rows = dsp.n_rows();
188 cols = dsp.n_cols();
189
190 // simple case: just use the given sparsity pattern
191 if (chunk_size == 1)
192 {
194 return;
195 }
196
197 // create a temporary compressed sparsity pattern that collects all entries
198 // from the input sparsity pattern and then initialize the underlying small
199 // sparsity pattern
200 const size_type m_chunks = (dsp.n_rows() + chunk_size - 1) / chunk_size,
201 n_chunks = (dsp.n_cols() + chunk_size - 1) / chunk_size;
202 DynamicSparsityPattern temporary_sp(m_chunks, n_chunks);
203
204 for (size_type row = 0; row < dsp.n_rows(); ++row)
205 {
206 const size_type reduced_row = row / chunk_size;
207
208 // TODO: This could be made more efficient if we cached the
209 // previous column and only called add() if the previous and the
210 // current column lead to different chunk columns
211 for (typename SparsityPatternType::iterator col_num = dsp.begin(row);
212 col_num != dsp.end(row);
213 ++col_num)
214 temporary_sp.add(reduced_row, col_num->column() / chunk_size);
215 }
216
217 sparsity_pattern.copy_from(temporary_sp);
218}
219
220
221
222template <typename number>
223void
225 const size_type chunk_size)
226{
228
229 // count number of entries per row, then initialize the underlying sparsity
230 // pattern. remember to also allocate space for the diagonal entry (if that
231 // hasn't happened yet) if m==n since we always allocate that for diagonal
232 // matrices
233 std::vector<size_type> entries_per_row(matrix.m(), 0);
234 for (size_type row = 0; row < matrix.m(); ++row)
235 {
236 for (size_type col = 0; col < matrix.n(); ++col)
237 if (matrix(row, col) != 0)
238 ++entries_per_row[row];
239
240 if ((matrix.m() == matrix.n()) && (matrix(row, row) == 0))
241 ++entries_per_row[row];
242 }
243
244 reinit(matrix.m(), matrix.n(), entries_per_row, chunk_size);
245
246 // then actually fill it
247 for (size_type row = 0; row < matrix.m(); ++row)
248 for (size_type col = 0; col < matrix.n(); ++col)
249 if (matrix(row, col) != 0)
250 add(row, col);
251
252 // finally compress
253 compress();
254}
255
256
257
258void
260 const size_type n,
261 const std::vector<size_type> &row_lengths,
262 const size_type chunk_size)
263{
265
266 reinit(m, n, make_array_view(row_lengths), chunk_size);
267}
268
269
270
271namespace internal
272{
273 namespace
274 {
275 template <typename SparsityPatternType>
276 void
277 copy_sparsity(const SparsityPatternType &src, SparsityPattern &dst)
278 {
279 dst.copy_from(src);
280 }
281
282 void
283 copy_sparsity(const SparsityPattern &src, SparsityPattern &dst)
284 {
285 dst = src;
286 }
287 } // namespace
288} // namespace internal
289
290
291
292template <typename Sparsity>
293void
295 const size_type n,
296 const Sparsity &sparsity_pattern_for_chunks,
297 const size_type chunk_size_in,
298 const bool)
299{
300 Assert(m > (sparsity_pattern_for_chunks.n_rows() - 1) * chunk_size_in &&
301 m <= sparsity_pattern_for_chunks.n_rows() * chunk_size_in,
302 ExcMessage("Number of rows m is not compatible with chunk size "
303 "and number of rows in sparsity pattern for the chunks."));
304 Assert(n > (sparsity_pattern_for_chunks.n_cols() - 1) * chunk_size_in &&
305 n <= sparsity_pattern_for_chunks.n_cols() * chunk_size_in,
307 "Number of columns m is not compatible with chunk size "
308 "and number of columns in sparsity pattern for the chunks."));
309
310 internal::copy_sparsity(sparsity_pattern_for_chunks, sparsity_pattern);
311 chunk_size = chunk_size_in;
312 rows = m;
313 cols = n;
314}
315
316
317
318bool
320{
321 return sparsity_pattern.empty();
322}
323
324
325
331
332
333
334void
342
343
344bool
346{
349
351}
352
353
354
355void
357{
358 // matrix must be square. note that the for some matrix sizes, the current
359 // sparsity pattern may not be square even if the underlying sparsity
360 // pattern is (e.g. a 10x11 matrix with chunk_size 4)
362
364}
365
366
367
370{
372
373 // find out if we did padding and if this row is affected by it
374 if (n_cols() % chunk_size == 0)
376 else
377 // if columns don't align, then just iterate over all chunks and see
378 // what this leads to
379 {
382 end =
384 unsigned int n = 0;
385 for (; p != end; ++p)
386 if (p->column() != sparsity_pattern.n_cols() - 1)
387 n += chunk_size;
388 else
389 n += (n_cols() % chunk_size);
390 return n;
391 }
392}
393
394
395
398{
399 if ((n_rows() % chunk_size == 0) && (n_cols() % chunk_size == 0))
401 else
402 // some of the chunks reach beyond the extent of this matrix. this
403 // requires a somewhat more complicated computations, in particular if the
404 // columns don't align
405 {
406 if ((n_rows() % chunk_size != 0) && (n_cols() % chunk_size == 0))
407 {
408 // columns align with chunks, but not rows
409 size_type n =
414 return n;
415 }
416
417 else
418 {
419 // if columns don't align, then just iterate over all chunks and see
420 // what this leads to. follow the advice in the documentation of the
421 // sparsity pattern iterators to do the loop over individual rows,
422 // rather than all elements
423 size_type n = 0;
424
425 for (size_type row = 0; row < sparsity_pattern.n_rows(); ++row)
426 {
428 for (; p != sparsity_pattern.end(row); ++p)
429 if ((row != sparsity_pattern.n_rows() - 1) &&
430 (p->column() != sparsity_pattern.n_cols() - 1))
431 n += chunk_size * chunk_size;
432 else if ((row == sparsity_pattern.n_rows() - 1) &&
433 (p->column() != sparsity_pattern.n_cols() - 1))
434 // last chunk row, but not last chunk column. only a smaller
435 // number (n_rows % chunk_size) of rows actually exist
436 n += (n_rows() % chunk_size) * chunk_size;
437 else if ((row != sparsity_pattern.n_rows() - 1) &&
438 (p->column() == sparsity_pattern.n_cols() - 1))
439 // last chunk column, but not row
440 n += (n_cols() % chunk_size) * chunk_size;
441 else
442 // bottom right chunk
443 n += (n_cols() % chunk_size) * (n_rows() % chunk_size);
444 }
445
446 return n;
447 }
448 }
449}
450
451
452
453void
454ChunkSparsityPattern::print(std::ostream &out) const
455{
456 Assert((sparsity_pattern.rowstart != nullptr) &&
457 (sparsity_pattern.colnums != nullptr),
459
460 AssertThrow(out.fail() == false, ExcIO());
461
462 for (size_type i = 0; i < sparsity_pattern.rows; ++i)
463 for (size_type d = 0; (d < chunk_size) && (i * chunk_size + d < n_rows());
464 ++d)
465 {
466 out << '[' << i * chunk_size + d;
468 j < sparsity_pattern.rowstart[i + 1];
469 ++j)
471 for (size_type e = 0;
472 ((e < chunk_size) &&
474 ++e)
475 out << ',' << sparsity_pattern.colnums[j] * chunk_size + e;
476 out << ']' << std::endl;
477 }
478
479 AssertThrow(out.fail() == false, ExcIO());
480}
481
482
483
484void
486{
487 Assert((sparsity_pattern.rowstart != nullptr) &&
488 (sparsity_pattern.colnums != nullptr),
490
491 AssertThrow(out.fail() == false, ExcIO());
492
493 // for each entry in the underlying sparsity pattern, repeat everything
494 // chunk_size x chunk_size times
495 for (size_type i = 0; i < sparsity_pattern.rows; ++i)
497 j < sparsity_pattern.rowstart[i + 1];
498 ++j)
500 for (size_type d = 0;
501 ((d < chunk_size) &&
503 ++d)
504 for (size_type e = 0;
505 (e < chunk_size) && (i * chunk_size + e < n_rows());
506 ++e)
507 // while matrix entries are usually written (i,j), with i vertical
508 // and j horizontal, gnuplot output is x-y, that is we have to
509 // exchange the order of output
510 out << sparsity_pattern.colnums[j] * chunk_size + d << " "
511 << -static_cast<signed int>(i * chunk_size + e) << std::endl;
512
513 AssertThrow(out.fail() == false, ExcIO());
514}
515
516
517
520{
521 // calculate the bandwidth from that of the underlying sparsity
522 // pattern. note that even if the bandwidth of that is zero, then the
523 // bandwidth of the chunky pattern is chunk_size-1, if it is 1 then the
524 // chunky pattern has chunk_size+(chunk_size-1), etc
525 //
526 // we'll cut it off at max(n(),m())
528 std::max(n_rows(), n_cols()));
529}
530
531
532
533bool
535{
536 if (chunk_size == 1)
538 else
539 return false;
540}
541
542
543
544void
545ChunkSparsityPattern::block_write(std::ostream &out) const
546{
547 AssertThrow(out.fail() == false, ExcIO());
548
549 // first the simple objects, bracketed in [...]
550 out << '[' << rows << ' ' << cols << ' ' << chunk_size << ' ' << "][";
551 // then the underlying sparsity pattern
553 out << ']';
554
555 AssertThrow(out.fail() == false, ExcIO());
556}
557
558
559
560void
562{
563 AssertThrow(in.fail() == false, ExcIO());
564
565 char c;
566
567 // first read in simple data
568 in >> c;
569 AssertThrow(c == '[', ExcIO());
570 in >> rows >> cols >> chunk_size;
571
572 in >> c;
573 AssertThrow(c == ']', ExcIO());
574 in >> c;
575 AssertThrow(c == '[', ExcIO());
576
577 // then read the underlying sparsity pattern
579
580 in >> c;
581 AssertThrow(c == ']', ExcIO());
582}
583
584
585
586std::size_t
588{
589 return (sizeof(*this) + sparsity_pattern.memory_consumption());
590}
591
592
593
594#ifndef DOXYGEN
595// explicit instantiations
596template void
599 const size_type);
600template void
602 const size_type,
603 const SparsityPattern &,
604 const size_type,
605 const bool);
606template void
608 const size_type,
609 const size_type,
611 const size_type,
612 const bool);
613template void
615 const size_type);
616template void
618 const size_type);
619#endif
620
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
std::size_t size() const
Definition array_view.h:684
void create_from(const size_type m, const size_type n, const Sparsity &sparsity_pattern_for_chunks, const size_type chunk_size, const bool optimize_diagonal=true)
void add(const size_type i, const size_type j)
void block_write(std::ostream &out) const
std::size_t memory_consumption() const
bool exists(const size_type i, const size_type j) const
iterator end() const
void reinit(const size_type m, const size_type n, const size_type max_per_row, const size_type chunk_size)
void print_gnuplot(std::ostream &out) const
void block_read(std::istream &in)
ChunkSparsityPattern & operator=(const ChunkSparsityPattern &)
size_type max_entries_per_row() const
size_type n_cols() const
size_type n_nonzero_elements() const
size_type row_length(const size_type row) const
void print(std::ostream &out) const
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end, const size_type chunk_size)
size_type n_rows() const
void add(const size_type i, const size_type j)
size_type n_rows() const
size_type n_cols() const
void block_write(std::ostream &out) const
void reinit(const size_type m, const size_type n, const ArrayView< const unsigned int > &row_lengths)
bool stores_only_added_elements() const
size_type bandwidth() const
std::size_t n_nonzero_elements() const
bool exists(const size_type i, const size_type j) const
iterator begin() const
std::unique_ptr< size_type[]> colnums
void copy_from(const size_type n_rows, const size_type n_cols, const ForwardIterator begin, const ForwardIterator end)
static constexpr size_type invalid_entry
std::unique_ptr< std::size_t[]> rowstart
void add(const size_type i, const size_type j)
size_type max_entries_per_row() const
iterator end() const
unsigned int row_length(const size_type row) const
std::size_t memory_consumption() const
void block_read(std::istream &in)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
static ::ExceptionBase & ExcIO()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcEmptyObject()
static ::ExceptionBase & ExcNotQuadratic()
static ::ExceptionBase & ExcInvalidIndex(size_type arg1, size_type arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
static ::ExceptionBase & ExcInvalidNumber(size_type arg1)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)