Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Details_InverseDiagonalKernel_def.hpp
1// @HEADER
2// *****************************************************************************
3// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
4//
5// Copyright 2009 NTESS and the Ifpack2 contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef IFPACK2_DETAILS_INVERSEDIAGONALKERNEL_DEF_HPP
11#define IFPACK2_DETAILS_INVERSEDIAGONALKERNEL_DEF_HPP
12
13#include "Tpetra_CrsMatrix.hpp"
14#include "Tpetra_MultiVector.hpp"
15#include "Tpetra_Operator.hpp"
16#include "Tpetra_Vector.hpp"
17#include "Tpetra_Export_decl.hpp"
18#include "Tpetra_Import_decl.hpp"
19#if KOKKOS_VERSION >= 40799
20#include "KokkosKernels_ArithTraits.hpp"
21#else
22#include "Kokkos_ArithTraits.hpp"
23#endif
24#include "Teuchos_Assert.hpp"
25#include <type_traits>
26#include "KokkosSparse_spmv_impl.hpp"
27
28namespace Ifpack2 {
29namespace Details {
30namespace Impl {
31
35template <class DVector,
36 class AMatrix,
37 class DiagOffsetType,
38 bool do_L1,
39 bool fix_tiny>
40struct InverseDiagonalWithExtraction {
41 using execution_space = typename AMatrix::execution_space;
42 using LO = typename AMatrix::non_const_ordinal_type;
43 using value_type = typename AMatrix::non_const_value_type;
44 using team_policy = typename Kokkos::TeamPolicy<execution_space>;
45 using team_member = typename team_policy::member_type;
46#if KOKKOS_VERSION >= 40799
47 using ATV = KokkosKernels::ArithTraits<value_type>;
48#else
49 using ATV = Kokkos::ArithTraits<value_type>;
50#endif
51 // using IST = typename vector_type::impl_scalar_type;
52 using magnitude_type = typename ATV::mag_type;
53#if KOKKOS_VERSION >= 40799
54 using MATV = KokkosKernels::ArithTraits<magnitude_type>;
55#else
56 using MATV = Kokkos::ArithTraits<magnitude_type>;
57#endif
58
59 DVector m_d;
60 AMatrix m_A;
61 DiagOffsetType m_offsets;
62 magnitude_type m_L1Eta;
63 magnitude_type m_MinDiagonalValue;
64
65 InverseDiagonalWithExtraction(const DVector& m_d_,
66 const AMatrix& m_A_,
67 const DiagOffsetType& m_offsets_,
68 const magnitude_type m_L1Eta_,
69 const magnitude_type m_MinDiagonalValue_)
70 : m_d(m_d_)
71 , m_A(m_A_)
72 , m_offsets(m_offsets_)
73 , m_L1Eta(m_L1Eta_)
74 , m_MinDiagonalValue(m_MinDiagonalValue_) {
75 const size_t numRows = m_A.numRows();
76
77 TEUCHOS_ASSERT(numRows == size_t(m_d.extent(0)));
78 TEUCHOS_ASSERT(numRows == size_t(m_offsets.extent(0)));
79 }
80
81 KOKKOS_INLINE_FUNCTION
82 void operator()(const LO lclRow) const {
83 const size_t INV = Tpetra::Details::OrdinalTraits<size_t>::invalid();
84 const value_type one = ATV::one();
85
86 // In case the row has no entries
87 m_d(lclRow, 0) = ATV::zero();
88
89 if (m_offsets(lclRow) != INV) {
90 auto curRow = m_A.rowConst(lclRow);
91 value_type d = curRow.value(m_offsets(lclRow));
92
93 if (do_L1) {
94 // Setup for L1 Methods.
95 // Here we add half the value of the off-processor entries in the row,
96 // but only if diagonal isn't sufficiently large.
97 //
98 // This follows from Equation (6.5) in: Baker, Falgout, Kolev and
99 // Yang. "Multigrid Smoothers for Ultraparallel Computing." SIAM
100 // J. Sci. Comput., Vol. 33, No. 5. (2011), pp. 2864-2887.
101
102 const magnitude_type half = MATV::one() / (MATV::one() + MATV::one());
103 const LO numRows = static_cast<LO>(m_A.numRows());
104 const LO row_length = static_cast<LO>(curRow.length);
105 magnitude_type diagonal_boost = MATV::zero();
106 for (LO iEntry = 0; iEntry < row_length; iEntry++) {
107 if (curRow.colidx(iEntry) >= numRows)
108 diagonal_boost += ATV::magnitude(curRow.value(iEntry));
109 }
110 diagonal_boost *= half;
111 if (ATV::magnitude(d) < m_L1Eta * diagonal_boost)
112 d += diagonal_boost;
113 }
114
115 if (fix_tiny) {
116 // Replace diagonal entries that are too small.
117
118 if (ATV::magnitude(d) <= m_MinDiagonalValue)
119 d = m_MinDiagonalValue;
120 }
121
122 // invert diagonal entries
123 m_d(lclRow, 0) = one / d;
124 }
125 }
126};
127
128} // namespace Impl
129
130template <class TpetraOperatorType>
131InverseDiagonalKernel<TpetraOperatorType>::
132 InverseDiagonalKernel(const Teuchos::RCP<const operator_type>& A) {
133 setMatrix(A);
134}
135
136template <class TpetraOperatorType>
137void InverseDiagonalKernel<TpetraOperatorType>::
138 setMatrix(const Teuchos::RCP<const operator_type>& A) {
139 if (A_op_.get() != A.get()) {
140 A_op_ = A;
141
142 using Teuchos::rcp_dynamic_cast;
143 A_crs_ = rcp_dynamic_cast<const crs_matrix_type>(A);
144
145 TEUCHOS_TEST_FOR_EXCEPTION(A_crs_.is_null(), std::logic_error,
146 "Ifpack2::Details::InverseDiagonalKernel: operator A must be a Tpetra::CrsMatrix.");
147
148 const size_t lclNumRows = A_crs_->getRowMap()->getLocalNumElements();
149
150 if (offsets_.extent(0) < lclNumRows) {
151 using Kokkos::view_alloc;
152 using Kokkos::WithoutInitializing;
153 using offsets_view_type = Kokkos::View<size_t*, device_type>;
154
155 offsets_ = offsets_view_type(); // clear 1st to save mem
156 auto howAlloc = view_alloc("offsets", WithoutInitializing);
157 offsets_ = offsets_view_type(howAlloc, lclNumRows);
158 }
159
160 A_crs_->getCrsGraph()->getLocalDiagOffsets(offsets_);
161 }
162}
163
164template <class TpetraOperatorType>
165void InverseDiagonalKernel<TpetraOperatorType>::
166 compute(vector_type& D_inv,
167 bool do_l1, magnitude_type L1Eta,
168 bool fixTinyDiagEntries, magnitude_type MinDiagonalValue) {
169 // Canonicalize template arguments to avoid redundant instantiations.
170 using d_type = typename vector_type::dual_view_type::t_dev;
171 // using h_matrix_type = typename crs_matrix_type::local_matrix_host_type;
172 using d_matrix_type = typename crs_matrix_type::local_matrix_device_type;
173
174 const char kernel_label[] = "inverse_diagonal_kernel";
175 using execution_space = typename NT::execution_space;
176 using range_type = Kokkos::RangePolicy<execution_space, LO>;
177 const size_t lclNumRows = A_crs_->getRowMap()->getLocalNumElements();
178 auto policy = range_type(0, lclNumRows);
179
180 d_type d = D_inv.getLocalViewDevice(Tpetra::Access::OverwriteAll);
181 d_matrix_type a = A_crs_->getLocalMatrixDevice();
182
183 if (do_l1) {
184 constexpr bool do_l1_template = true;
185 if (fixTinyDiagEntries) {
186 constexpr bool fix_tiny_template = true;
187 using functor_type =
188 Impl::InverseDiagonalWithExtraction<d_type,
189 d_matrix_type,
190 offset_type,
191 do_l1_template,
192 fix_tiny_template>;
193 functor_type func(d, a, offsets_, L1Eta, MinDiagonalValue);
194 Kokkos::parallel_for(kernel_label, policy, func);
195 } else {
196 constexpr bool fix_tiny_template = false;
197 using functor_type =
198 Impl::InverseDiagonalWithExtraction<d_type,
199 d_matrix_type,
200 offset_type,
201 do_l1_template,
202 fix_tiny_template>;
203 functor_type func(d, a, offsets_, L1Eta, MinDiagonalValue);
204 Kokkos::parallel_for(kernel_label, policy, func);
205 }
206 } else {
207 constexpr bool do_l1_template = false;
208 if (fixTinyDiagEntries) {
209 constexpr bool fix_tiny_template = true;
210 using functor_type =
211 Impl::InverseDiagonalWithExtraction<d_type,
212 d_matrix_type,
213 offset_type,
214 do_l1_template,
215 fix_tiny_template>;
216 functor_type func(d, a, offsets_, L1Eta, MinDiagonalValue);
217 Kokkos::parallel_for(kernel_label, policy, func);
218 } else {
219 constexpr bool fix_tiny_template = false;
220 using functor_type =
221 Impl::InverseDiagonalWithExtraction<d_type,
222 d_matrix_type,
223 offset_type,
224 do_l1_template,
225 fix_tiny_template>;
226 functor_type func(d, a, offsets_, L1Eta, MinDiagonalValue);
227 Kokkos::parallel_for(kernel_label, policy, func);
228 }
229 }
230}
231
232} // namespace Details
233} // namespace Ifpack2
234
235#define IFPACK2_DETAILS_INVERSEDIAGONALKERNEL_INSTANT(SC, LO, GO, NT) \
236 template class Ifpack2::Details::InverseDiagonalKernel<Tpetra::Operator<SC, LO, GO, NT> >;
237
238#endif // IFPACK2_DETAILS_INVERSEDIAGONALKERNEL_DEF_HPP
Ifpack2::Details::InverseDiagonalKernel
Compute scaled damped residual for Chebyshev.
Definition Ifpack2_Details_InverseDiagonalKernel_decl.hpp:45
Ifpack2::Details::LinearSolver
Ifpack2's implementation of Trilinos::Details::LinearSolver interface.
Definition Ifpack2_Details_LinearSolver_decl.hpp:75
Ifpack2::Details::LinearSolver::setMatrix
void setMatrix(const Teuchos::RCP< const OP > &A)
Set the solver's matrix.
Definition Ifpack2_Details_LinearSolver_def.hpp:53
Details
Ifpack2 implementation details.
Ifpack2
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:40
Ifpack2::Details::Impl::InverseDiagonalWithExtraction
Functor for extracting the inverse diagonal of a matrix.
Definition Ifpack2_Details_InverseDiagonalKernel_def.hpp:40
Generated by doxygen 1.9.8