Reference documentation for deal.II version GIT 29f9da0a34 2023-12-07 10:00: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\}}\)
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
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 #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
29 template <typename number>
30 class Vector;
31 namespace LinearAlgebra
32 {
33  namespace distributed
34  {
35  template <typename, typename>
36  class Vector;
37  } // namespace distributed
38 } // namespace LinearAlgebra
39 #endif
40 
61 template <typename VectorType = Vector<double>>
63 {
64 public:
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
108  clear();
109 
113  const VectorType &
114  get_vector() const;
115 
120  size_type
121  m() const;
122 
127  size_type
128  n() const;
129 
139  value_type
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
167  add(const size_type row,
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 
218  value_type
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 
252 private:
256  VectorType diagonal;
257 };
258 
259 /* ---------------------------------- Inline functions ------------------- */
260 
261 #ifndef DOXYGEN
262 
263 template <typename VectorType>
264 DiagonalMatrix<VectorType>::DiagonalMatrix(const VectorType &vec)
265  : diagonal(vec)
266 {}
267 
268 
269 
270 template <typename VectorType>
271 void
273 {
274  diagonal.reinit(0);
275 }
276 
277 
278 
279 template <typename VectorType>
280 std::size_t
282 {
283  return diagonal.memory_consumption();
284 }
285 
286 
287 
288 template <typename VectorType>
289 void
290 DiagonalMatrix<VectorType>::reinit(const VectorType &vec)
291 {
292  diagonal = vec;
293 }
294 
295 
296 
297 template <typename VectorType>
298 void
300 {
301  dst.reinit(diagonal);
302 }
303 
304 
305 
306 template <typename VectorType>
307 void
309 {
310  diagonal.compress(operation);
311 }
312 
313 
314 
315 template <typename VectorType>
316 VectorType &
318 {
319  return diagonal;
320 }
321 
322 
323 
324 template <typename VectorType>
325 const VectorType &
327 {
328  return diagonal;
329 }
330 
331 
332 
333 template <typename VectorType>
334 typename VectorType::size_type
336 {
337  return diagonal.size();
338 }
339 
340 
341 
342 template <typename VectorType>
343 typename VectorType::size_type
345 {
346  return diagonal.size();
347 }
348 
349 
350 
351 template <typename VectorType>
352 typename 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 
363 template <typename VectorType>
364 typename VectorType::value_type &
366 {
367  Assert(i == j, ExcIndexRange(j, i, i + 1));
368  (void)j;
369  return diagonal(i);
370 }
371 
372 
373 
374 template <typename VectorType>
375 template <typename number2>
376 void
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 
391 template <typename VectorType>
392 void
394  const size_type j,
395  const value_type value)
396 {
397  if (i == j)
398  diagonal(i) += value;
399 }
400 
401 
402 
403 namespace 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  auto *const dst_ptr = dst.begin();
426  const auto *const src_ptr = src.begin();
427  const auto *const 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 
438 template <typename VectorType>
439 void
440 DiagonalMatrix<VectorType>::vmult(VectorType &dst, const VectorType &src) const
441 {
442  internal::DiagonalMatrix::assign_and_scale(dst, src, diagonal);
443 }
444 
445 
446 
447 template <typename VectorType>
448 void
449 DiagonalMatrix<VectorType>::Tvmult(VectorType &dst, const VectorType &src) const
450 {
451  vmult(dst, src);
452 }
453 
454 
455 
456 template <typename VectorType>
457 void
459  const VectorType &src) const
460 {
461  VectorType tmp(src);
462  tmp.scale(diagonal);
463  dst += tmp;
464 }
465 
466 
467 
468 template <typename VectorType>
469 void
471  const VectorType &src) const
472 {
473  vmult_add(dst, src);
474 }
475 
476 
477 
478 template <typename VectorType>
479 typename VectorType::value_type
480 DiagonalMatrix<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 
489 template <typename VectorType>
490 void
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
VectorType & get_vector()
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
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
value_type & operator()(const size_type i, const size_type j)
typename VectorType::size_type size_type
void compress(VectorOperation::values operation)
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
const VectorType & get_vector() const
size_type locally_owned_size() const
Definition: vector.h:110
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:144
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
#define Assert(cond, exc)
Definition: exceptions.h:1631
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
@ diagonal
Matrix is diagonal.
types::global_dof_index size_type
Definition: cuda_kernels.h:45