Reference documentation for deal.II version GIT relicensing-422-gb369f187d8 2024-04-17 18:10: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
diagonal_matrix.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2016 - 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#ifndef dealii_diagonal_matrix_h
16#define dealii_diagonal_matrix_h
17
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/lac/vector.h>
23
25
26// forward declarations
27#ifndef DOXYGEN
28template <typename number>
29class Vector;
30namespace LinearAlgebra
31{
32 namespace distributed
33 {
34 template <typename, typename>
35 class Vector;
36 } // namespace distributed
37} // namespace LinearAlgebra
38#endif
39
60template <typename VectorType = Vector<double>>
62{
63public:
64 using value_type = typename VectorType::value_type;
65 using size_type = typename VectorType::size_type;
66
71 DiagonalMatrix() = default;
72
78 explicit DiagonalMatrix(const VectorType &vec);
79
84 void
85 reinit(const VectorType &vec);
86
93 void
95
100 VectorType &
102
106 void
108
112 const VectorType &
113 get_vector() const;
114
120 m() const;
121
127 n() const;
128
139 operator()(const size_type i, const size_type j) const;
140
150 value_type &
151 operator()(const size_type i, const size_type j);
152
164 template <typename number2>
165 void
166 add(const size_type row,
167 const size_type n_cols,
168 const size_type *col_indices,
169 const number2 *values,
170 const bool elide_zero_values = true,
171 const bool col_indices_are_sorted = false);
172
179 void
180 add(const size_type i, const size_type j, const value_type value);
181
185 void
186 vmult(VectorType &dst, const VectorType &src) const;
187
193 void
194 Tvmult(VectorType &dst, const VectorType &src) const;
195
201 void
202 vmult_add(VectorType &dst, const VectorType &src) const;
203
209 void
210 Tvmult_add(VectorType &dst, const VectorType &src) const;
211
218 apply(const unsigned int index, const value_type src) const;
219
229 void
230 apply_to_subrange(const unsigned int begin_range,
231 const unsigned int end_range,
232 const value_type *src_pointer_to_current_range,
233 value_type *dst_pointer_to_current_range) const;
234
242 void
243 initialize_dof_vector(VectorType &dst) const;
244
248 std::size_t
250
251private:
255 VectorType diagonal;
256};
257
258/* ---------------------------------- Inline functions ------------------- */
259
260#ifndef DOXYGEN
261
262template <typename VectorType>
264 : diagonal(vec)
265{}
266
267
268
269template <typename VectorType>
270void
272{
273 diagonal.reinit(0);
274}
275
276
277
278template <typename VectorType>
279std::size_t
281{
282 return diagonal.memory_consumption();
283}
284
285
286
287template <typename VectorType>
288void
289DiagonalMatrix<VectorType>::reinit(const VectorType &vec)
290{
291 diagonal = vec;
292}
293
294
295
296template <typename VectorType>
297void
299{
300 dst.reinit(diagonal);
301}
302
303
304
305template <typename VectorType>
306void
308{
309 diagonal.compress(operation);
310}
311
312
313
314template <typename VectorType>
315VectorType &
317{
318 return diagonal;
319}
320
321
322
323template <typename VectorType>
324const VectorType &
326{
327 return diagonal;
328}
329
330
331
332template <typename VectorType>
333typename VectorType::size_type
335{
336 return diagonal.size();
337}
338
339
340
341template <typename VectorType>
342typename VectorType::size_type
344{
345 return diagonal.size();
346}
347
348
349
350template <typename VectorType>
351typename VectorType::value_type
353 const size_type j) const
354{
355 Assert(i == j, ExcIndexRange(j, i, i + 1));
356 (void)j;
357 return diagonal(i);
358}
359
360
361
362template <typename VectorType>
363typename VectorType::value_type &
364DiagonalMatrix<VectorType>::operator()(const size_type i, const size_type j)
365{
366 Assert(i == j, ExcIndexRange(j, i, i + 1));
367 (void)j;
368 return diagonal(i);
369}
370
371
372
373template <typename VectorType>
374template <typename number2>
375void
376DiagonalMatrix<VectorType>::add(const size_type row,
377 const size_type n_cols,
378 const size_type *col_indices,
379 const number2 *values,
380 const bool,
381 const bool)
382{
383 for (size_type i = 0; i < n_cols; ++i)
384 if (col_indices[i] == row)
385 diagonal(row) += values[i];
386}
387
388
389
390template <typename VectorType>
391void
392DiagonalMatrix<VectorType>::add(const size_type i,
393 const size_type j,
394 const value_type value)
395{
396 if (i == j)
397 diagonal(i) += value;
398}
399
400
401
402namespace internal
403{
404 namespace DiagonalMatrix
405 {
406 template <typename VectorType>
407 void
408 assign_and_scale(VectorType &dst,
409 const VectorType &src,
410 const VectorType &diagonal)
411 {
412 dst = src;
413 dst.scale(diagonal);
414 }
415
416 template <typename Number>
417 void
418 assign_and_scale(
422 &diagonal)
423 {
424 auto *const dst_ptr = dst.begin();
425 const auto *const src_ptr = src.begin();
426 const auto *const diagonal_ptr = diagonal.begin();
427
429 for (unsigned int i = 0; i < src.locally_owned_size(); ++i)
430 dst_ptr[i] = src_ptr[i] * diagonal_ptr[i];
431 }
432 } // namespace DiagonalMatrix
433} // namespace internal
434
435
436
437template <typename VectorType>
438void
439DiagonalMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
440{
441 internal::DiagonalMatrix::assign_and_scale(dst, src, diagonal);
442}
443
444
445
446template <typename VectorType>
447void
448DiagonalMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
449{
450 vmult(dst, src);
451}
452
453
454
455template <typename VectorType>
456void
458 const VectorType &src) const
459{
460 VectorType tmp(src);
461 tmp.scale(diagonal);
462 dst += tmp;
463}
464
465
466
467template <typename VectorType>
468void
470 const VectorType &src) const
471{
472 vmult_add(dst, src);
473}
474
475
476
477template <typename VectorType>
478typename VectorType::value_type
479DiagonalMatrix<VectorType>::apply(const unsigned int index,
480 const value_type src) const
481{
482 AssertIndexRange(index, diagonal.locally_owned_elements().n_elements());
483 return diagonal.local_element(index) * src;
484}
485
486
487
488template <typename VectorType>
489void
491 const unsigned int begin_range,
492 const unsigned int end_range,
493 const value_type *src_pointer_to_current_range,
494 value_type *dst_pointer_to_current_range) const
495{
496 AssertIndexRange(begin_range,
497 diagonal.locally_owned_elements().n_elements() + 1);
498 AssertIndexRange(end_range,
499 diagonal.locally_owned_elements().n_elements() + 1);
500
501 const value_type *diagonal_entry = diagonal.begin() + begin_range;
502 const unsigned int length = end_range - begin_range;
503
505 for (unsigned int i = 0; i < length; ++i)
506 dst_pointer_to_current_range[i] =
507 diagonal_entry[i] * src_pointer_to_current_range[i];
508}
509
510
511#endif
512
514
515#endif
void initialize_dof_vector(VectorType &dst) const
size_type m() const
void Tvmult_add(VectorType &dst, const VectorType &src) const
value_type apply(const unsigned int index, const value_type src) const
void apply_to_subrange(const unsigned int begin_range, const unsigned int end_range, const value_type *src_pointer_to_current_range, value_type *dst_pointer_to_current_range) const
VectorType diagonal
value_type & operator()(const size_type i, const size_type j)
VectorType & get_vector()
std::size_t memory_consumption() const
void add(const size_type i, const size_type j, const value_type value)
DiagonalMatrix(const VectorType &vec)
value_type operator()(const size_type i, const size_type j) const
typename VectorType::size_type size_type
void compress(VectorOperation::values operation)
const VectorType & get_vector() const
void vmult(VectorType &dst, const VectorType &src) const
size_type n() const
void reinit(const VectorType &vec)
void add(const size_type row, const size_type n_cols, const size_type *col_indices, const number2 *values, const bool elide_zero_values=true, const bool col_indices_are_sorted=false)
typename VectorType::value_type value_type
void vmult_add(VectorType &dst, const VectorType &src) const
DiagonalMatrix()=default
void Tvmult(VectorType &dst, const VectorType &src) const
size_type locally_owned_size() const
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:143
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
@ diagonal
Matrix is diagonal.