deal.II version GIT relicensing-3507-gc50d9d6ce3 2025-06-12 13:20:01+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>
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 return diagonal(i);
357}
358
359
360
361template <typename VectorType>
362typename VectorType::value_type &
363DiagonalMatrix<VectorType>::operator()(const size_type i, const size_type j)
364{
365 Assert(i == j, ExcIndexRange(j, i, i + 1));
366 return diagonal(i);
367}
368
369
370
371template <typename VectorType>
372template <typename number2>
373void
374DiagonalMatrix<VectorType>::add(const size_type row,
375 const size_type n_cols,
376 const size_type *col_indices,
377 const number2 *values,
378 const bool,
379 const bool)
380{
381 for (size_type i = 0; i < n_cols; ++i)
382 if (col_indices[i] == row)
383 diagonal(row) += values[i];
384}
385
386
387
388template <typename VectorType>
389void
390DiagonalMatrix<VectorType>::add(const size_type i,
391 const size_type j,
392 const value_type value)
393{
394 if (i == j)
395 diagonal(i) += value;
396}
397
398
399
400namespace internal
401{
402 namespace DiagonalMatrix
403 {
404 template <typename VectorType>
405 void
406 assign_and_scale(VectorType &dst,
407 const VectorType &src,
408 const VectorType &diagonal)
409 {
410 dst = src;
411 dst.scale(diagonal);
412 }
413
414 template <typename Number>
415 void
416 assign_and_scale(
420 &diagonal)
421 {
422 auto *const dst_ptr = dst.begin();
423 const auto *const src_ptr = src.begin();
424 const auto *const diagonal_ptr = diagonal.begin();
425
427 for (unsigned int i = 0; i < src.locally_owned_size(); ++i)
428 dst_ptr[i] = src_ptr[i] * diagonal_ptr[i];
429 }
430 } // namespace DiagonalMatrix
431} // namespace internal
432
433
434
435template <typename VectorType>
436void
437DiagonalMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
438{
439 internal::DiagonalMatrix::assign_and_scale(dst, src, diagonal);
440}
441
442
443
444template <typename VectorType>
445void
446DiagonalMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
447{
448 vmult(dst, src);
449}
450
451
452
453template <typename VectorType>
454void
456 const VectorType &src) const
457{
458 VectorType tmp(src);
459 tmp.scale(diagonal);
460 dst += tmp;
461}
462
463
464
465template <typename VectorType>
466void
468 const VectorType &src) const
469{
470 vmult_add(dst, src);
471}
472
473
474
475template <typename VectorType>
476typename VectorType::value_type
477DiagonalMatrix<VectorType>::apply(const unsigned int index,
478 const value_type src) const
479{
480 AssertIndexRange(index, diagonal.locally_owned_elements().n_elements());
481 return diagonal.local_element(index) * src;
482}
483
484
485
486template <typename VectorType>
487void
489 const unsigned int begin_range,
490 const unsigned int end_range,
491 const value_type *src_pointer_to_current_range,
492 value_type *dst_pointer_to_current_range) const
493{
494 AssertIndexRange(begin_range,
495 diagonal.locally_owned_elements().n_elements() + 1);
496 AssertIndexRange(end_range,
497 diagonal.locally_owned_elements().n_elements() + 1);
498
499 const value_type *diagonal_entry = diagonal.begin() + begin_range;
500 const unsigned int length = end_range - begin_range;
501
503 for (unsigned int i = 0; i < length; ++i)
504 dst_pointer_to_current_range[i] =
505 diagonal_entry[i] * src_pointer_to_current_range[i];
506}
507
508
509#endif
510
512
513#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:209
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#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.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType