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