Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Kokkos_MV_UQ_PCE.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef KOKKOS_MV_UQ_PCE_HPP
43#define KOKKOS_MV_UQ_PCE_HPP
44
45#include "Sacado_UQ_PCE.hpp"
49/*
50//----------------------------------------------------------------------------
51// Specializations of Kokkos Vector/MultiVector math functions
52//----------------------------------------------------------------------------
53
54namespace Kokkos {
55
56// Rank-1 vector add with Sacado::UQ::PCE scalar type, constant a, b
57template <typename RS, typename RL, typename RD, typename RM,
58 typename XS, typename XL, typename XD, typename XM,
59 typename YS, typename YL, typename YD, typename YM>
60Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM>
61V_Add( const Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM >& r,
62 const typename Sacado::UQ::PCE<XS>::value_type& av,
63 const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x,
64 const typename Sacado::UQ::PCE<XS>::value_type& bv,
65 const Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM >& y,
66 int n = -1)
67{
68 typedef Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM > RVector;
69 typedef Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM > XVector;
70 typedef Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM > YVector;
71
72 typename RVector::flat_array_type r_flat = r;
73 typename XVector::flat_array_type x_flat = x;
74 typename YVector::flat_array_type y_flat = y;
75 if (n != -1) n = n * r.sacado_size();
76
77 V_Add( r_flat, av, x_flat, bv, y_flat, n );
78
79 return r;
80}
81
82// Rank-1 vector add with Sacado::UQ::PCE scalar type, non-constant a, b
83template <typename RS, typename RL, typename RD, typename RM,
84 typename XS, typename XL, typename XD, typename XM,
85 typename YS, typename YL, typename YD, typename YM>
86Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM>
87V_Add( const Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM >& r,
88 const Sacado::UQ::PCE<XS>& av,
89 const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x,
90 const Sacado::UQ::PCE<XS>& bv,
91 const Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM >& y,
92 int n = -1)
93{
94 if (Sacado::is_constant(av) && Sacado::is_constant(bv)) {
95 return V_Add( r, av.fastAccessCoeff(0), x, bv.fastAccessCoeff(0), y, n );
96 }
97 else {
98 Impl::raise_error("V_Add not implemented for non-constant a or b");
99 }
100 return r;
101}
102
103// Rank-2 vector add with Sacado::UQ::PCE scalar type, constant a, b
104template <typename RS, typename RL, typename RD, typename RM,
105 typename XS, typename XL, typename XD, typename XM,
106 typename YS, typename YL, typename YD, typename YM>
107Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM>
108MV_Add( const Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM >& r,
109 const typename Sacado::UQ::PCE<XS>::value_type& av,
110 const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x,
111 const typename Sacado::UQ::PCE<XS>::value_type& bv,
112 const Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM >& y,
113 int n = -1)
114{
115 typedef Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM > RVector;
116 typedef Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM > XVector;
117 typedef Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM > YVector;
118
119 typename RVector::flat_array_type r_flat = r;
120 typename XVector::flat_array_type x_flat = x;
121 typename YVector::flat_array_type y_flat = y;
122 if (n != -1) n = n * r.sacado_size();
123
124 MV_Add( r_flat, av, x_flat, bv, y_flat, n );
125
126 return r;
127}
128
129// Rank-2 vector add with Sacado::UQ::PCE scalar type, non-constant a, b
130template <typename RS, typename RL, typename RD, typename RM,
131 typename XS, typename XL, typename XD, typename XM,
132 typename YS, typename YL, typename YD, typename YM>
133Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM>
134MV_Add( const Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM >& r,
135 const Sacado::UQ::PCE<XS>& av,
136 const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x,
137 const Sacado::UQ::PCE<XS>& bv,
138 const Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM >& y,
139 int n = -1)
140{
141 if (Sacado::is_constant(av) && Sacado::is_constant(bv)) {
142 return MV_Add( r, av.fastAccessCoeff(0), x, bv.fastAccessCoeff(0), y, n );
143 }
144 else {
145 Impl::raise_error("MV_Add not implemented for non-constant a or b");
146 }
147 return r;
148}
149
150// Rank-1 dot product
151template <typename XS, typename XL, typename XD, typename XM,
152 typename YS, typename YL, typename YD, typename YM>
153typename Details::InnerProductSpaceTraits< Sacado::UQ::PCE<XS> >::dot_type
154V_Dot( const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x,
155 const Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM >& y,
156 int n = -1 )
157{
158 typedef Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM > XVector;
159 typedef Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM > YVector;
160
161 typename XVector::flat_array_type x_flat = x;
162 typename YVector::flat_array_type y_flat = y;
163 if (n != -1) n = n * x.sacado_size();
164
165 return V_Dot( x_flat, y_flat, n );
166}
167
168// Rank-2 dot product
169template <typename rVector,
170 typename XS, typename XL, typename XD, typename XM,
171 typename YS, typename YL, typename YD, typename YM>
172void
173MV_Dot( const rVector& r,
174 const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x,
175 const Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM >& y,
176 int n = -1 )
177{
178 typedef Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM > XVector;
179 typedef Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM > YVector;
180
181 typename XVector::flat_array_type x_flat = x;
182 typename YVector::flat_array_type y_flat = y;
183 if (n != -1) n = n * x.sacado_size();
184
185 MV_Dot( r, x_flat, y_flat, n );
186}
187
188template<class VT1, class VT2, class VT3>
189struct MV_ElementWiseMultiplyFunctor;
190
191template <typename CT, typename CD, typename CM,
192 typename AT, typename AD, typename AM,
193 typename BT, typename BD, typename BM>
194struct MV_ElementWiseMultiplyFunctor<
195 View< CT,LayoutLeft,CD,CM,Impl::ViewPCEContiguous >,
196 View< AT,LayoutLeft,AD,AM,Impl::ViewPCEContiguous >,
197 View< BT,LayoutLeft,BD,BM,Impl::ViewPCEContiguous > >
198{
199 typedef View< CT,LayoutLeft,CD,CM,Impl::ViewPCEContiguous > CVector;
200 typedef View< AT,LayoutLeft,AD,AM,Impl::ViewPCEContiguous > AVector;
201 typedef View< BT,LayoutLeft,BD,BM,Impl::ViewPCEContiguous > BVector;
202
203 typedef typename CVector::array_type CArray;
204 typedef typename AVector::array_type AArray;
205 typedef typename BVector::array_type BArray;
206
207 typedef typename CArray::execution_space execution_space;
208 typedef typename CArray::size_type size_type;
209
210 typename CArray::const_value_type m_c;
211 CArray m_C;
212 typename AArray::const_value_type m_ab;
213 typename AArray::const_type m_A ;
214 typename BArray::const_type m_B ;
215 const size_type m_n;
216 const size_type m_n_pce;
217
218 MV_ElementWiseMultiplyFunctor(typename CVector::const_value_type c,
219 CVector C,
220 typename AVector::const_value_type ab,
221 typename AVector::const_type A,
222 typename BVector::const_type B,
223 const size_type n) :
224 m_c(c.fastAccessCoeff(0)),
225 m_C(C),
226 m_ab(ab.fastAccessCoeff(0)),
227 m_A(A),
228 m_B(B),
229 m_n(n),
230 m_n_pce(C.sacado_size()) {}
231
232 KOKKOS_INLINE_FUNCTION
233 void operator()( const size_type i) const
234 {
235 // Currently specialized for use case where A is degree-0
236 typename AArray::const_value_type Ai = m_A(0,i);
237 for (size_type k=0; k<m_n; ++k) {
238#ifdef KOKKOS_ENABLE_PRAGMA_IVDEP
239#pragma ivdep
240#endif
241 for (size_type l=0; l<m_n_pce; ++l)
242 m_C(l,i,k) = m_c*m_C(l,i,k) + m_ab*Ai*m_B(l,i,k);
243 }
244 }
245};
246
247template<class VT1, class VT2, class VT3>
248struct V_ElementWiseMultiplyFunctor;
249
250template <typename CT, typename CD, typename CM,
251 typename AT, typename AD, typename AM,
252 typename BT, typename BD, typename BM>
253struct V_ElementWiseMultiplyFunctor<
254 View< CT,LayoutLeft,CD,CM,Impl::ViewPCEContiguous >,
255 View< AT,LayoutLeft,AD,AM,Impl::ViewPCEContiguous >,
256 View< BT,LayoutLeft,BD,BM,Impl::ViewPCEContiguous > >
257{
258 typedef View< CT,LayoutLeft,CD,CM,Impl::ViewPCEContiguous > CVector;
259 typedef View< AT,LayoutLeft,AD,AM,Impl::ViewPCEContiguous > AVector;
260 typedef View< BT,LayoutLeft,BD,BM,Impl::ViewPCEContiguous > BVector;
261
262 typedef typename CVector::array_type CArray;
263 typedef typename AVector::array_type AArray;
264 typedef typename BVector::array_type BArray;
265
266 typedef typename CArray::execution_space execution_space;
267 typedef typename CArray::size_type size_type;
268
269 typename CArray::const_value_type m_c;
270 CArray m_C;
271 typename AArray::const_value_type m_ab;
272 typename AArray::const_type m_A ;
273 typename BArray::const_type m_B ;
274 const size_type m_n_pce;
275
276 V_ElementWiseMultiplyFunctor(typename CVector::const_value_type c,
277 CVector C,
278 typename AVector::const_value_type ab,
279 typename AVector::const_type A,
280 typename BVector::const_type B) :
281 m_c(c.fastAccessCoeff(0)),
282 m_C(C),
283 m_ab(ab.fastAccessCoeff(0)),
284 m_A(A),
285 m_B(B),
286 m_n_pce(C.sacado_size()) {}
287
288 KOKKOS_INLINE_FUNCTION
289 void operator()( const size_type i) const
290 {
291 // Currently specialized for use case where A is degree-0
292 typename AArray::const_value_type Ai = m_A(0,i);
293#ifdef KOKKOS_ENABLE_PRAGMA_IVDEP
294#pragma ivdep
295#endif
296 for (size_type l=0; l<m_n_pce; ++l)
297 m_C(l,i) = m_c*m_C(l,i) + m_ab*Ai*m_B(l,i);
298 }
299};
300
301// Rank-1 vector multiply with Sacado::UQ::PCE scalar type, constant c, ab
302template <typename CS, typename CL, typename CD, typename CM,
303 typename AS, typename AL, typename AD, typename AM,
304 typename BS, typename BL, typename BD, typename BM>
305Kokkos::View< Sacado::UQ::PCE<CS>*, CL, CD, CM>
306V_ElementWiseMultiply(
307 const typename Sacado::UQ::PCE<CS>::value_type& c,
308 const Kokkos::View< Sacado::UQ::PCE<CS>*, CL, CD, CM >& C,
309 const typename Sacado::UQ::PCE<AS>::value_type& ab,
310 const Kokkos::View< Sacado::UQ::PCE<AS>*, AL, AD, AM >& A,
311 const Kokkos::View< Sacado::UQ::PCE<BS>*, BL, BD, BM >& B )
312{
313 typedef View< Sacado::UQ::PCE<CS>*, CL, CD, CM > CVector;
314 typedef View< Sacado::UQ::PCE<AS>*, AL, AD, AM > AVector;
315 typedef View< Sacado::UQ::PCE<BS>*, BL, BD, BM > BVector;
316
317 typedef View< typename CVector::data_type, typename CVector::array_layout, typename CVector::execution_space, typename CVector::memory_traits > CView;
318 typedef View< typename AVector::data_type, typename AVector::array_layout, typename AVector::execution_space, typename AVector::memory_traits > AView;
319 typedef View< typename BVector::data_type, typename BVector::array_layout, typename BVector::execution_space, typename BVector::memory_traits > BView;
320
321 V_ElementWiseMultiplyFunctor<CView,AView,BView> op(c,C,ab,A,B) ;
322 Kokkos::parallel_for( C.extent(0) , op );
323 return C;
324}
325
326// Rank-1 vector multiply with Sacado::UQ::PCE scalar type, non-constant c, ab
327template <typename CS, typename CL, typename CD, typename CM,
328 typename AS, typename AL, typename AD, typename AM,
329 typename BS, typename BL, typename BD, typename BM>
330Kokkos::View< Sacado::UQ::PCE<CS>*, CL, CD, CM>
331V_ElementWiseMultiply(
332 const Sacado::UQ::PCE<CS>& c,
333 const Kokkos::View< Sacado::UQ::PCE<CS>*, CL, CD, CM >& C,
334 const Sacado::UQ::PCE<AS>& ab,
335 const Kokkos::View< Sacado::UQ::PCE<AS>*, AL, AD, AM >& A,
336 const Kokkos::View< Sacado::UQ::PCE<BS>*, BL, BD, BM >& B )
337{
338 if (Sacado::is_constant(c) && Sacado::is_constant(ab)) {
339 return V_ElementWiseMultiply( c.fastAccessCoeff(0), C,
340 ab.fastAccessCoeff(0), A, B );
341 }
342 else {
343 Impl::raise_error(
344 "V_ElementWiseMultiply not implemented for non-constant c or ab");
345 }
346 return C;
347}
348
349// Rank-2 vector multiply with Sacado::UQ::PCE scalar type, constant c, ab
350template <typename CS, typename CL, typename CD, typename CM,
351 typename AS, typename AL, typename AD, typename AM,
352 typename BS, typename BL, typename BD, typename BM>
353Kokkos::View< Sacado::UQ::PCE<CS>**, CL, CD, CM>
354MV_ElementWiseMultiply(
355 const typename Sacado::UQ::PCE<CS>::value_type& c,
356 const Kokkos::View< Sacado::UQ::PCE<CS>**, CL, CD, CM >& C,
357 const typename Sacado::UQ::PCE<AS>::value_type& ab,
358 const Kokkos::View< Sacado::UQ::PCE<AS>*, AL, AD, AM >& A,
359 const Kokkos::View< Sacado::UQ::PCE<BS>**, BL, BD, BM >& B )
360{
361 typedef View< Sacado::UQ::PCE<CS>**, CL, CD, CM > CVector;
362 typedef View< Sacado::UQ::PCE<AS>*, AL, AD, AM > AVector;
363 typedef View< Sacado::UQ::PCE<BS>**, BL, BD, BM > BVector;
364
365 typedef View< typename CVector::data_type, typename CVector::array_layout, typename CVector::execution_space, typename CVector::memory_traits > CView;
366 typedef View< typename AVector::data_type, typename AVector::array_layout, typename AVector::execution_space, typename AVector::memory_traits > AView;
367 typedef View< typename BVector::data_type, typename BVector::array_layout, typename BVector::execution_space, typename BVector::memory_traits > BView;
368
369 MV_ElementWiseMultiplyFunctor<CView,AView,BView> op(c,C,ab,A,B,C.extent(1)) ;
370 Kokkos::parallel_for( C.extent(0) , op );
371 return C;
372}
373
374// Rank-2 vector multiply with Sacado::UQ::PCE scalar type, non-constant c, ab
375template <typename CS, typename CL, typename CD, typename CM,
376 typename AS, typename AL, typename AD, typename AM,
377 typename BS, typename BL, typename BD, typename BM>
378Kokkos::View< Sacado::UQ::PCE<CS>**, CL, CD, CM>
379MV_ElementWiseMultiply(
380 const Sacado::UQ::PCE<CS>& c,
381 const Kokkos::View< Sacado::UQ::PCE<CS>**, CL, CD, CM >& C,
382 const Sacado::UQ::PCE<AS>& ab,
383 const Kokkos::View< Sacado::UQ::PCE<AS>*, AL, AD, AM >& A,
384 const Kokkos::View< Sacado::UQ::PCE<BS>**, BL, BD, BM >& B )
385{
386 if (Sacado::is_constant(c) && Sacado::is_constant(ab)) {
387 return MV_ElementWiseMultiply( c.fastAccessCoeff(0), C,
388 ab.fastAccessCoeff(0), A, B );
389 }
390 else {
391 Impl::raise_error(
392 "MV_ElementWiseMultiply not implemented for non-constant c or ab");
393 }
394 return C;
395}
396
397// Rank-1 vector scale with Sacado::UQ::PCE scalar type, constant a, b
398template <typename RS, typename RL, typename RD, typename RM,
399 typename XS, typename XL, typename XD, typename XM>
400Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM>
401V_MulScalar( const Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM >& r,
402 const typename Sacado::UQ::PCE<XS>::value_type& a,
403 const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x )
404{
405 typedef Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM > RVector;
406 typedef Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM > XVector;
407
408 typename RVector::flat_array_type r_flat = r;
409 typename XVector::flat_array_type x_flat = x;
410
411 V_MulScalar( r_flat, a, x_flat );
412
413 return r;
414}
415
416// Rank-1 vector scale with Sacado::UQ::PCE scalar type, non-constant a, b
417template <typename RS, typename RL, typename RD, typename RM,
418 typename XS, typename XL, typename XD, typename XM>
419Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM>
420V_MulScalar( const Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM >& r,
421 const Sacado::UQ::PCE<XS>& a,
422 const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x )
423{
424 if (Sacado::is_constant(a)) {
425 return V_MulScalar( r, a.fastAccessCoeff(0), x );
426 }
427 else {
428 Impl::raise_error("V_MulScalar not implemented for non-constant a");
429 }
430 return r;
431}
432
433// Rank-2 vector scale with Sacado::UQ::PCE scalar type, constant a, b
434template <typename RS, typename RL, typename RD, typename RM,
435 typename XS, typename XL, typename XD, typename XM>
436Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM>
437MV_MulScalar( const Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM >& r,
438 const typename Sacado::UQ::PCE<XS>::value_type& a,
439 const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x )
440{
441 typedef Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM > RVector;
442 typedef Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM > XVector;
443
444 typename RVector::flat_array_type r_flat = r;
445 typename XVector::flat_array_type x_flat = x;
446
447 MV_MulScalar( r_flat, a, x_flat );
448
449 return r;
450}
451
452// Rank-2 vector scale with Sacado::UQ::PCE scalar type, non-constant a, b
453template <typename RS, typename RL, typename RD, typename RM,
454 typename XS, typename XL, typename XD, typename XM>
455Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM>
456MV_MulScalar( const Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM >& r,
457 const Sacado::UQ::PCE<XS>& a,
458 const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x )
459{
460 if (Sacado::is_constant(a)) {
461 return MV_MulScalar( r, a.fastAccessCoeff(0), x );
462 }
463 else {
464 Impl::raise_error("MV_MulScalar not implemented for non-constant a or b");
465 }
466 return r;
467}
468
469template <typename T>
470struct V_ReciprocalThresholdSelfFunctor;
471
472template <typename T, typename D, typename M>
473struct V_ReciprocalThresholdSelfFunctor<
474 View< T,LayoutLeft,D,M,Impl::ViewPCEContiguous > >
475{
476 typedef View< T,LayoutLeft,D,M,Impl::ViewPCEContiguous > XVector;
477 typedef typename XVector::array_type array_type;
478
479 typedef typename array_type::execution_space execution_space;
480 typedef typename array_type::size_type size_type;
481 typedef typename array_type::non_const_value_type value_type;
482 typedef Kokkos::Details::ArithTraits<value_type> KAT;
483 typedef typename KAT::mag_type mag_type;
484
485 const array_type m_x;
486 const value_type m_min_val;
487 const value_type m_min_val_mag;
488 const size_type m_n_pce;
489
490 V_ReciprocalThresholdSelfFunctor(
491 const XVector& x,
492 const typename XVector::non_const_value_type& min_val) :
493 m_x(x),
494 m_min_val(min_val.fastAccessCoeff(0)),
495 m_min_val_mag(KAT::abs(m_min_val)),
496 m_n_pce(x.sacado_size()) {}
497 //--------------------------------------------------------------------------
498
499 KOKKOS_INLINE_FUNCTION
500 void operator()( const size_type i) const
501 {
502 if (KAT::abs(m_x(0,i)) < m_min_val_mag)
503 m_x(0,i) = m_min_val;
504 else
505 m_x(0,i) = KAT::one() / m_x(0,i);
506 for (size_type l=1; l<m_n_pce; ++l)
507 m_x(l,i) = KAT::zero();
508 }
509};
510
511template <typename T>
512struct MV_ReciprocalThresholdSelfFunctor;
513
514template <typename T, typename D, typename M>
515struct MV_ReciprocalThresholdSelfFunctor<
516 View< T,LayoutLeft,D,M,Impl::ViewPCEContiguous > >
517{
518 typedef View< T,LayoutLeft,D,M,Impl::ViewPCEContiguous > XVector;
519 typedef typename XVector::array_type array_type;
520
521 typedef typename array_type::execution_space execution_space;
522 typedef typename array_type::size_type size_type;
523 typedef typename array_type::non_const_value_type value_type;
524 typedef Kokkos::Details::ArithTraits<value_type> KAT;
525 typedef typename KAT::mag_type mag_type;
526
527 const array_type m_x;
528 const value_type m_min_val;
529 const value_type m_min_val_mag;
530 const size_type m_n;
531 const size_type m_n_pce;
532
533 MV_ReciprocalThresholdSelfFunctor(
534 const XVector& x,
535 const typename XVector::non_const_value_type& min_val,
536 const size_type& n) :
537 m_x(x),
538 m_min_val(min_val.fastAccessCoeff(0)),
539 m_min_val_mag(KAT::abs(m_min_val)),
540 m_n(n),
541 m_n_pce(x.sacado_size()) {}
542 //--------------------------------------------------------------------------
543
544 KOKKOS_INLINE_FUNCTION
545 void operator()( const size_type i) const
546 {
547 for (size_type k=0; k<m_n; ++k) {
548 if (KAT::abs(m_x(0,i,k)) < m_min_val_mag)
549 m_x(0,i,k) = m_min_val;
550 else
551 m_x(0,i,k) = KAT::one() / m_x(0,i,k);
552 for (size_type l=1; l<m_n_pce; ++l)
553 m_x(l,i,k) = KAT::zero();
554 }
555 }
556};
557
558} // namespace Kokkos
559*/
560#endif /* #ifndef KOKKOS_MV_UQ_PCE_HPP */