Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
Ifpack2_Details_InverseDiagonalKernel_def.hpp
1/*
2//@HEADER
3// ***********************************************************************
4//
5// Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
6// Copyright (2009) Sandia Corporation
7//
8// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9// license for use of this work by or on behalf of the U.S. Government.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// ***********************************************************************
39//@HEADER
40*/
41
42#ifndef IFPACK2_DETAILS_INVERSEDIAGONALKERNEL_DEF_HPP
43#define IFPACK2_DETAILS_INVERSEDIAGONALKERNEL_DEF_HPP
44
45#include "Tpetra_CrsMatrix.hpp"
46#include "Tpetra_MultiVector.hpp"
47#include "Tpetra_Operator.hpp"
48#include "Tpetra_Vector.hpp"
49#include "Tpetra_Export_decl.hpp"
50#include "Tpetra_Import_decl.hpp"
51#include "Kokkos_ArithTraits.hpp"
52#include "Teuchos_Assert.hpp"
53#include <type_traits>
54#include "KokkosSparse_spmv_impl.hpp"
55
56namespace Ifpack2 {
57namespace Details {
58namespace Impl {
59
63template<class DVector,
64 class AMatrix,
65 class DiagOffsetType,
66 bool do_L1,
67 bool fix_tiny>
69
70 using execution_space = typename AMatrix::execution_space;
71 using LO = typename AMatrix::non_const_ordinal_type;
72 using value_type = typename AMatrix::non_const_value_type;
73 using team_policy = typename Kokkos::TeamPolicy<execution_space>;
74 using team_member = typename team_policy::member_type;
75 using ATV = Kokkos::ArithTraits<value_type>;
76 // using IST = typename vector_type::impl_scalar_type;
77 using magnitude_type = typename ATV::mag_type;
78 using MATV = Kokkos::ArithTraits<magnitude_type>;
79
80 DVector m_d;
81 AMatrix m_A;
82 DiagOffsetType m_offsets;
83 magnitude_type m_L1Eta;
84 magnitude_type m_MinDiagonalValue;
85
86 InverseDiagonalWithExtraction (const DVector& m_d_,
87 const AMatrix& m_A_,
88 const DiagOffsetType& m_offsets_,
89 const magnitude_type m_L1Eta_,
90 const magnitude_type m_MinDiagonalValue_) :
91 m_d (m_d_),
92 m_A (m_A_),
93 m_offsets (m_offsets_),
94 m_L1Eta (m_L1Eta_),
95 m_MinDiagonalValue (m_MinDiagonalValue_)
96 {
97 const size_t numRows = m_A.numRows ();
98
99 TEUCHOS_ASSERT( numRows == size_t (m_d.extent (0)) );
100 TEUCHOS_ASSERT( numRows == size_t (m_offsets.extent (0)) );
101 }
102
103 KOKKOS_INLINE_FUNCTION
104 void operator() (const LO lclRow) const
105 {
106 const size_t INV = Tpetra::Details::OrdinalTraits<size_t>::invalid ();
107 const value_type one = ATV::one();
108
109 // In case the row has no entries
110 m_d(lclRow,0) = ATV::zero();
111
112 if (m_offsets(lclRow) != INV) {
113 auto curRow = m_A.rowConst (lclRow);
114 value_type d = curRow.value(m_offsets(lclRow));
115
116 if (do_L1) {
117 // Setup for L1 Methods.
118 // Here we add half the value of the off-processor entries in the row,
119 // but only if diagonal isn't sufficiently large.
120 //
121 // This follows from Equation (6.5) in: Baker, Falgout, Kolev and
122 // Yang. "Multigrid Smoothers for Ultraparallel Computing." SIAM
123 // J. Sci. Comput., Vol. 33, No. 5. (2011), pp. 2864-2887.
124
125 const magnitude_type half = MATV::one () / (MATV::one () + MATV::one ());
126 const LO numRows = static_cast<LO> (m_A.numRows ());
127 const LO row_length = static_cast<LO> (curRow.length);
128 magnitude_type diagonal_boost = MATV::zero();
129 for (LO iEntry = 0; iEntry < row_length; iEntry++) {
130 if (curRow.colidx(iEntry) >= numRows)
131 diagonal_boost += ATV::magnitude(curRow.value(iEntry));
132 }
133 diagonal_boost *= half;
134 if (ATV::magnitude(d) < m_L1Eta * diagonal_boost)
135 d += diagonal_boost;
136 }
137
138 if (fix_tiny) {
139 // Replace diagonal entries that are too small.
140
141 if (ATV::magnitude(d) <= m_MinDiagonalValue)
142 d = m_MinDiagonalValue;
143 }
144
145 // invert diagonal entries
146 m_d(lclRow,0) = one / d;
147 }
148 }
149
150};
151
152} // namespace Impl
153
154
155template<class TpetraOperatorType>
157InverseDiagonalKernel (const Teuchos::RCP<const operator_type>& A)
158{
159 setMatrix (A);
160}
161
162template<class TpetraOperatorType>
163void
164InverseDiagonalKernel<TpetraOperatorType>::
165setMatrix (const Teuchos::RCP<const operator_type>& A)
166{
167 if (A_op_.get () != A.get ()) {
168 A_op_ = A;
169
170 using Teuchos::rcp_dynamic_cast;
171 A_crs_ = rcp_dynamic_cast<const crs_matrix_type> (A);
172
173 TEUCHOS_TEST_FOR_EXCEPTION
174 (A_crs_.is_null(), std::logic_error,
175 "Ifpack2::Details::InverseDiagonalKernel: operator A must be a Tpetra::CrsMatrix.");
176
177 const size_t lclNumRows = A_crs_->getRowMap ()->getLocalNumElements ();
178
179 if (offsets_.extent (0) < lclNumRows) {
180 using Kokkos::view_alloc;
181 using Kokkos::WithoutInitializing;
182 using offsets_view_type = Kokkos::View<size_t*, device_type>;
183
184 offsets_ = offsets_view_type (); // clear 1st to save mem
185 auto howAlloc = view_alloc ("offsets", WithoutInitializing);
186 offsets_ = offsets_view_type (howAlloc, lclNumRows);
187 }
188
189 A_crs_->getCrsGraph ()->getLocalDiagOffsets (offsets_);
190 }
191}
192
193template<class TpetraOperatorType>
194void
195InverseDiagonalKernel<TpetraOperatorType>::
196compute (vector_type& D_inv,
197 bool do_l1, magnitude_type L1Eta,
198 bool fixTinyDiagEntries, magnitude_type MinDiagonalValue)
199{
200
201
202 // Canonicalize template arguments to avoid redundant instantiations.
203 using d_type = typename vector_type::dual_view_type::t_dev;
204 // using h_matrix_type = typename crs_matrix_type::local_matrix_host_type;
205 using d_matrix_type = typename crs_matrix_type::local_matrix_device_type;
206
207 const char kernel_label[] = "inverse_diagonal_kernel";
208 using execution_space = typename NT::execution_space;
209 using range_type = Kokkos::RangePolicy<execution_space, LO>;
210 const size_t lclNumRows = A_crs_->getRowMap ()->getLocalNumElements ();
211 auto policy = range_type(0, lclNumRows);
212
213 d_type d = D_inv.getLocalViewDevice(Tpetra::Access::OverwriteAll);
214 d_matrix_type a = A_crs_->getLocalMatrixDevice();
215
216 if (do_l1) {
217 constexpr bool do_l1_template = true;
218 if (fixTinyDiagEntries) {
219 constexpr bool fix_tiny_template = true;
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 } else {
229 constexpr bool fix_tiny_template = false;
230 using functor_type =
231 Impl::InverseDiagonalWithExtraction<d_type,
232 d_matrix_type,
233 offset_type,
234 do_l1_template,
235 fix_tiny_template>;
236 functor_type func (d, a, offsets_, L1Eta, MinDiagonalValue);
237 Kokkos::parallel_for (kernel_label, policy, func);
238 }
239 } else {
240 constexpr bool do_l1_template = false;
241 if (fixTinyDiagEntries) {
242 constexpr bool fix_tiny_template = true;
243 using functor_type =
244 Impl::InverseDiagonalWithExtraction<d_type,
245 d_matrix_type,
246 offset_type,
247 do_l1_template,
248 fix_tiny_template>;
249 functor_type func (d, a, offsets_, L1Eta, MinDiagonalValue);
250 Kokkos::parallel_for (kernel_label, policy, func);
251 } else {
252 constexpr bool fix_tiny_template = false;
253 using functor_type =
254 Impl::InverseDiagonalWithExtraction<d_type,
255 d_matrix_type,
256 offset_type,
257 do_l1_template,
258 fix_tiny_template>;
259 functor_type func (d, a, offsets_, L1Eta, MinDiagonalValue);
260 Kokkos::parallel_for (kernel_label, policy, func);
261 }
262 }
263}
264
265} // namespace Details
266} // namespace Ifpack2
267
268#define IFPACK2_DETAILS_INVERSEDIAGONALKERNEL_INSTANT(SC,LO,GO,NT) \
269 template class Ifpack2::Details::InverseDiagonalKernel<Tpetra::Operator<SC, LO, GO, NT> >;
270
271#endif // IFPACK2_DETAILS_INVERSEDIAGONALKERNEL_DEF_HPP
Compute scaled damped residual for Chebyshev.
Definition Ifpack2_Details_InverseDiagonalKernel_decl.hpp:77
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition Ifpack2_AdditiveSchwarz_decl.hpp:74
Functor for extracting the inverse diagonal of a matrix.
Definition Ifpack2_Details_InverseDiagonalKernel_def.hpp:68