Reference documentation for deal.II version 9.5.0
Searching...
No Matches
diagonal_matrix.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2016 - 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
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#ifndef dealii_diagonal_matrix_h
17#define dealii_diagonal_matrix_h
18
19
20#include <deal.II/base/config.h>
21
22#include <deal.II/lac/vector.h>
24
26
27// forward declarations
28#ifndef DOXYGEN
29template <typename number>
30class Vector;
31namespace LinearAlgebra
32{
33 namespace distributed
34 {
35 template <typename, typename>
36 class Vector;
37 } // namespace distributed
38} // namespace LinearAlgebra
39#endif
40
61template <typename VectorType = Vector<double>>
63{
64public:
65 using value_type = typename VectorType::value_type;
66 using size_type = typename VectorType::size_type;
67
72 DiagonalMatrix() = default;
73
79 explicit DiagonalMatrix(const VectorType &vec);
80
85 void
86 reinit(const VectorType &vec);
87
94 void
96
101 VectorType &
103
107 void
109
113 const VectorType &
114 get_vector() const;
115
121 m() const;
122
128 n() const;
129
140 operator()(const size_type i, const size_type j) const;
141
151 value_type &
152 operator()(const size_type i, const size_type j);
153
165 template <typename number2>
166 void
168 const size_type n_cols,
169 const size_type *col_indices,
170 const number2 * values,
171 const bool elide_zero_values = true,
172 const bool col_indices_are_sorted = false);
173
180 void
181 add(const size_type i, const size_type j, const value_type value);
182
186 void
187 vmult(VectorType &dst, const VectorType &src) const;
188
194 void
195 Tvmult(VectorType &dst, const VectorType &src) const;
196
202 void
203 vmult_add(VectorType &dst, const VectorType &src) const;
204
210 void
211 Tvmult_add(VectorType &dst, const VectorType &src) const;
212
219 apply(const unsigned int index, const value_type src) const;
220
230 void
231 apply_to_subrange(const unsigned int begin_range,
232 const unsigned int end_range,
233 const value_type * src_pointer_to_current_range,
234 value_type * dst_pointer_to_current_range) const;
235
243 void
244 initialize_dof_vector(VectorType &dst) const;
245
249 std::size_t
251
252private:
256 VectorType diagonal;
257};
258
259/* ---------------------------------- Inline functions ------------------- */
260
261#ifndef DOXYGEN
262
263template <typename VectorType>
265 : diagonal(vec)
266{}
267
268
269
270template <typename VectorType>
271void
273{
274 diagonal.reinit(0);
275}
276
277
278
279template <typename VectorType>
280std::size_t
282{
283 return diagonal.memory_consumption();
284}
285
286
287
288template <typename VectorType>
289void
290DiagonalMatrix<VectorType>::reinit(const VectorType &vec)
291{
292 diagonal = vec;
293}
294
295
296
297template <typename VectorType>
298void
300{
301 dst.reinit(diagonal);
302}
303
304
305
306template <typename VectorType>
307void
309{
310 diagonal.compress(operation);
311}
312
313
314
315template <typename VectorType>
316VectorType &
318{
319 return diagonal;
320}
321
322
323
324template <typename VectorType>
325const VectorType &
327{
328 return diagonal;
329}
330
331
332
333template <typename VectorType>
334typename VectorType::size_type
336{
337 return diagonal.size();
338}
339
340
341
342template <typename VectorType>
343typename VectorType::size_type
345{
346 return diagonal.size();
347}
348
349
350
351template <typename VectorType>
352typename VectorType::value_type
354 const size_type j) const
355{
356 Assert(i == j, ExcIndexRange(j, i, i + 1));
357 (void)j;
358 return diagonal(i);
359}
360
361
362
363template <typename VectorType>
364typename VectorType::value_type &
365DiagonalMatrix<VectorType>::operator()(const size_type i, const size_type j)
366{
367 Assert(i == j, ExcIndexRange(j, i, i + 1));
368 (void)j;
369 return diagonal(i);
370}
371
372
373
374template <typename VectorType>
375template <typename number2>
376void
378 const size_type n_cols,
379 const size_type *col_indices,
380 const number2 * values,
381 const bool,
382 const bool)
383{
384 for (size_type i = 0; i < n_cols; ++i)
385 if (col_indices[i] == row)
386 diagonal(row) += values[i];
387}
388
389
390
391template <typename VectorType>
392void
394 const size_type j,
395 const value_type value)
396{
397 if (i == j)
398 diagonal(i) += value;
399}
400
401
402
403namespace internal
404{
405 namespace DiagonalMatrix
406 {
407 template <typename VectorType>
408 void
409 assign_and_scale(VectorType & dst,
410 const VectorType &src,
411 const VectorType &diagonal)
412 {
413 dst = src;
414 dst.scale(diagonal);
415 }
416
417 template <typename Number>
418 void
419 assign_and_scale(
423 &diagonal)
424 {
425 const auto dst_ptr = dst.begin();
426 const auto src_ptr = src.begin();
427 const auto diagonal_ptr = diagonal.begin();
428
430 for (unsigned int i = 0; i < src.locally_owned_size(); ++i)
431 dst_ptr[i] = src_ptr[i] * diagonal_ptr[i];
432 }
433 } // namespace DiagonalMatrix
434} // namespace internal
435
436
437
438template <typename VectorType>
439void
440DiagonalMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
441{
442 internal::DiagonalMatrix::assign_and_scale(dst, src, diagonal);
443}
444
445
446
447template <typename VectorType>
448void
449DiagonalMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
450{
451 vmult(dst, src);
452}
453
454
455
456template <typename VectorType>
457void
459 const VectorType &src) const
460{
461 VectorType tmp(src);
462 tmp.scale(diagonal);
463 dst += tmp;
464}
465
466
467
468template <typename VectorType>
469void
471 const VectorType &src) const
472{
474}
475
476
477
478template <typename VectorType>
479typename VectorType::value_type
480DiagonalMatrix<VectorType>::apply(const unsigned int index,
481 const value_type src) const
482{
483 AssertIndexRange(index, diagonal.locally_owned_elements().n_elements());
484 return diagonal.local_element(index) * src;
485}
486
487
488
489template <typename VectorType>
490void
492 const unsigned int begin_range,
493 const unsigned int end_range,
494 const value_type * src_pointer_to_current_range,
495 value_type * dst_pointer_to_current_range) const
496{
497 AssertIndexRange(begin_range,
498 diagonal.locally_owned_elements().n_elements() + 1);
499 AssertIndexRange(end_range,
500 diagonal.locally_owned_elements().n_elements() + 1);
501
502 const value_type * diagonal_entry = diagonal.begin() + begin_range;
503 const unsigned int length = end_range - begin_range;
504
506 for (unsigned int i = 0; i < length; ++i)
507 dst_pointer_to_current_range[i] =
508 diagonal_entry[i] * src_pointer_to_current_range[i];
509}
510
511
512#endif
513
515
516#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:138
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#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.