Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Kokkos_MV_UQ_PCE.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef KOKKOS_MV_UQ_PCE_HPP
11 #define KOKKOS_MV_UQ_PCE_HPP
12 
13 #include "Sacado_UQ_PCE.hpp"
14 #include "Kokkos_View_UQ_PCE.hpp"
16 #include "Kokkos_Blas1_UQ_PCE.hpp"
17 /*
18 //----------------------------------------------------------------------------
19 // Specializations of Kokkos Vector/MultiVector math functions
20 //----------------------------------------------------------------------------
21 
22 namespace Kokkos {
23 
24 // Rank-1 vector add with Sacado::UQ::PCE scalar type, constant a, b
25 template <typename RS, typename RL, typename RD, typename RM,
26  typename XS, typename XL, typename XD, typename XM,
27  typename YS, typename YL, typename YD, typename YM>
28 Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM>
29 V_Add( const Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM >& r,
30  const typename Sacado::UQ::PCE<XS>::value_type& av,
31  const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x,
32  const typename Sacado::UQ::PCE<XS>::value_type& bv,
33  const Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM >& y,
34  int n = -1)
35 {
36  typedef Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM > RVector;
37  typedef Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM > XVector;
38  typedef Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM > YVector;
39 
40  typename RVector::flat_array_type r_flat = r;
41  typename XVector::flat_array_type x_flat = x;
42  typename YVector::flat_array_type y_flat = y;
43  if (n != -1) n = n * r.sacado_size();
44 
45  V_Add( r_flat, av, x_flat, bv, y_flat, n );
46 
47  return r;
48 }
49 
50 // Rank-1 vector add with Sacado::UQ::PCE scalar type, non-constant a, b
51 template <typename RS, typename RL, typename RD, typename RM,
52  typename XS, typename XL, typename XD, typename XM,
53  typename YS, typename YL, typename YD, typename YM>
54 Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM>
55 V_Add( const Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM >& r,
56  const Sacado::UQ::PCE<XS>& av,
57  const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x,
58  const Sacado::UQ::PCE<XS>& bv,
59  const Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM >& y,
60  int n = -1)
61 {
62  if (Sacado::is_constant(av) && Sacado::is_constant(bv)) {
63  return V_Add( r, av.fastAccessCoeff(0), x, bv.fastAccessCoeff(0), y, n );
64  }
65  else {
66  Impl::raise_error("V_Add not implemented for non-constant a or b");
67  }
68  return r;
69 }
70 
71 // Rank-2 vector add with Sacado::UQ::PCE scalar type, constant a, b
72 template <typename RS, typename RL, typename RD, typename RM,
73  typename XS, typename XL, typename XD, typename XM,
74  typename YS, typename YL, typename YD, typename YM>
75 Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM>
76 MV_Add( const Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM >& r,
77  const typename Sacado::UQ::PCE<XS>::value_type& av,
78  const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x,
79  const typename Sacado::UQ::PCE<XS>::value_type& bv,
80  const Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM >& y,
81  int n = -1)
82 {
83  typedef Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM > RVector;
84  typedef Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM > XVector;
85  typedef Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM > YVector;
86 
87  typename RVector::flat_array_type r_flat = r;
88  typename XVector::flat_array_type x_flat = x;
89  typename YVector::flat_array_type y_flat = y;
90  if (n != -1) n = n * r.sacado_size();
91 
92  MV_Add( r_flat, av, x_flat, bv, y_flat, n );
93 
94  return r;
95 }
96 
97 // Rank-2 vector add with Sacado::UQ::PCE scalar type, non-constant a, b
98 template <typename RS, typename RL, typename RD, typename RM,
99  typename XS, typename XL, typename XD, typename XM,
100  typename YS, typename YL, typename YD, typename YM>
101 Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM>
102 MV_Add( const Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM >& r,
103  const Sacado::UQ::PCE<XS>& av,
104  const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x,
105  const Sacado::UQ::PCE<XS>& bv,
106  const Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM >& y,
107  int n = -1)
108 {
109  if (Sacado::is_constant(av) && Sacado::is_constant(bv)) {
110  return MV_Add( r, av.fastAccessCoeff(0), x, bv.fastAccessCoeff(0), y, n );
111  }
112  else {
113  Impl::raise_error("MV_Add not implemented for non-constant a or b");
114  }
115  return r;
116 }
117 
118 // Rank-1 dot product
119 template <typename XS, typename XL, typename XD, typename XM,
120  typename YS, typename YL, typename YD, typename YM>
121 typename Details::InnerProductSpaceTraits< Sacado::UQ::PCE<XS> >::dot_type
122 V_Dot( const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x,
123  const Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM >& y,
124  int n = -1 )
125 {
126  typedef Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM > XVector;
127  typedef Kokkos::View< Sacado::UQ::PCE<YS>*, YL, YD, YM > YVector;
128 
129  typename XVector::flat_array_type x_flat = x;
130  typename YVector::flat_array_type y_flat = y;
131  if (n != -1) n = n * x.sacado_size();
132 
133  return V_Dot( x_flat, y_flat, n );
134 }
135 
136 // Rank-2 dot product
137 template <typename rVector,
138  typename XS, typename XL, typename XD, typename XM,
139  typename YS, typename YL, typename YD, typename YM>
140 void
141 MV_Dot( const rVector& r,
142  const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x,
143  const Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM >& y,
144  int n = -1 )
145 {
146  typedef Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM > XVector;
147  typedef Kokkos::View< Sacado::UQ::PCE<YS>**, YL, YD, YM > YVector;
148 
149  typename XVector::flat_array_type x_flat = x;
150  typename YVector::flat_array_type y_flat = y;
151  if (n != -1) n = n * x.sacado_size();
152 
153  MV_Dot( r, x_flat, y_flat, n );
154 }
155 
156 template<class VT1, class VT2, class VT3>
157 struct MV_ElementWiseMultiplyFunctor;
158 
159 template <typename CT, typename CD, typename CM,
160  typename AT, typename AD, typename AM,
161  typename BT, typename BD, typename BM>
162 struct MV_ElementWiseMultiplyFunctor<
163  View< CT,LayoutLeft,CD,CM,Impl::ViewPCEContiguous >,
164  View< AT,LayoutLeft,AD,AM,Impl::ViewPCEContiguous >,
165  View< BT,LayoutLeft,BD,BM,Impl::ViewPCEContiguous > >
166 {
167  typedef View< CT,LayoutLeft,CD,CM,Impl::ViewPCEContiguous > CVector;
168  typedef View< AT,LayoutLeft,AD,AM,Impl::ViewPCEContiguous > AVector;
169  typedef View< BT,LayoutLeft,BD,BM,Impl::ViewPCEContiguous > BVector;
170 
171  typedef typename CVector::array_type CArray;
172  typedef typename AVector::array_type AArray;
173  typedef typename BVector::array_type BArray;
174 
175  typedef typename CArray::execution_space execution_space;
176  typedef typename CArray::size_type size_type;
177 
178  typename CArray::const_value_type m_c;
179  CArray m_C;
180  typename AArray::const_value_type m_ab;
181  typename AArray::const_type m_A ;
182  typename BArray::const_type m_B ;
183  const size_type m_n;
184  const size_type m_n_pce;
185 
186  MV_ElementWiseMultiplyFunctor(typename CVector::const_value_type c,
187  CVector C,
188  typename AVector::const_value_type ab,
189  typename AVector::const_type A,
190  typename BVector::const_type B,
191  const size_type n) :
192  m_c(c.fastAccessCoeff(0)),
193  m_C(C),
194  m_ab(ab.fastAccessCoeff(0)),
195  m_A(A),
196  m_B(B),
197  m_n(n),
198  m_n_pce(C.sacado_size()) {}
199 
200  KOKKOS_INLINE_FUNCTION
201  void operator()( const size_type i) const
202  {
203  // Currently specialized for use case where A is degree-0
204  typename AArray::const_value_type Ai = m_A(0,i);
205  for (size_type k=0; k<m_n; ++k) {
206 #ifdef KOKKOS_ENABLE_PRAGMA_IVDEP
207 #pragma ivdep
208 #endif
209  for (size_type l=0; l<m_n_pce; ++l)
210  m_C(l,i,k) = m_c*m_C(l,i,k) + m_ab*Ai*m_B(l,i,k);
211  }
212  }
213 };
214 
215 template<class VT1, class VT2, class VT3>
216 struct V_ElementWiseMultiplyFunctor;
217 
218 template <typename CT, typename CD, typename CM,
219  typename AT, typename AD, typename AM,
220  typename BT, typename BD, typename BM>
221 struct V_ElementWiseMultiplyFunctor<
222  View< CT,LayoutLeft,CD,CM,Impl::ViewPCEContiguous >,
223  View< AT,LayoutLeft,AD,AM,Impl::ViewPCEContiguous >,
224  View< BT,LayoutLeft,BD,BM,Impl::ViewPCEContiguous > >
225 {
226  typedef View< CT,LayoutLeft,CD,CM,Impl::ViewPCEContiguous > CVector;
227  typedef View< AT,LayoutLeft,AD,AM,Impl::ViewPCEContiguous > AVector;
228  typedef View< BT,LayoutLeft,BD,BM,Impl::ViewPCEContiguous > BVector;
229 
230  typedef typename CVector::array_type CArray;
231  typedef typename AVector::array_type AArray;
232  typedef typename BVector::array_type BArray;
233 
234  typedef typename CArray::execution_space execution_space;
235  typedef typename CArray::size_type size_type;
236 
237  typename CArray::const_value_type m_c;
238  CArray m_C;
239  typename AArray::const_value_type m_ab;
240  typename AArray::const_type m_A ;
241  typename BArray::const_type m_B ;
242  const size_type m_n_pce;
243 
244  V_ElementWiseMultiplyFunctor(typename CVector::const_value_type c,
245  CVector C,
246  typename AVector::const_value_type ab,
247  typename AVector::const_type A,
248  typename BVector::const_type B) :
249  m_c(c.fastAccessCoeff(0)),
250  m_C(C),
251  m_ab(ab.fastAccessCoeff(0)),
252  m_A(A),
253  m_B(B),
254  m_n_pce(C.sacado_size()) {}
255 
256  KOKKOS_INLINE_FUNCTION
257  void operator()( const size_type i) const
258  {
259  // Currently specialized for use case where A is degree-0
260  typename AArray::const_value_type Ai = m_A(0,i);
261 #ifdef KOKKOS_ENABLE_PRAGMA_IVDEP
262 #pragma ivdep
263 #endif
264  for (size_type l=0; l<m_n_pce; ++l)
265  m_C(l,i) = m_c*m_C(l,i) + m_ab*Ai*m_B(l,i);
266  }
267 };
268 
269 // Rank-1 vector multiply with Sacado::UQ::PCE scalar type, constant c, ab
270 template <typename CS, typename CL, typename CD, typename CM,
271  typename AS, typename AL, typename AD, typename AM,
272  typename BS, typename BL, typename BD, typename BM>
273 Kokkos::View< Sacado::UQ::PCE<CS>*, CL, CD, CM>
274 V_ElementWiseMultiply(
275  const typename Sacado::UQ::PCE<CS>::value_type& c,
276  const Kokkos::View< Sacado::UQ::PCE<CS>*, CL, CD, CM >& C,
277  const typename Sacado::UQ::PCE<AS>::value_type& ab,
278  const Kokkos::View< Sacado::UQ::PCE<AS>*, AL, AD, AM >& A,
279  const Kokkos::View< Sacado::UQ::PCE<BS>*, BL, BD, BM >& B )
280 {
281  typedef View< Sacado::UQ::PCE<CS>*, CL, CD, CM > CVector;
282  typedef View< Sacado::UQ::PCE<AS>*, AL, AD, AM > AVector;
283  typedef View< Sacado::UQ::PCE<BS>*, BL, BD, BM > BVector;
284 
285  typedef View< typename CVector::data_type, typename CVector::array_layout, typename CVector::execution_space, typename CVector::memory_traits > CView;
286  typedef View< typename AVector::data_type, typename AVector::array_layout, typename AVector::execution_space, typename AVector::memory_traits > AView;
287  typedef View< typename BVector::data_type, typename BVector::array_layout, typename BVector::execution_space, typename BVector::memory_traits > BView;
288 
289  V_ElementWiseMultiplyFunctor<CView,AView,BView> op(c,C,ab,A,B) ;
290  Kokkos::parallel_for( C.extent(0) , op );
291  return C;
292 }
293 
294 // Rank-1 vector multiply with Sacado::UQ::PCE scalar type, non-constant c, ab
295 template <typename CS, typename CL, typename CD, typename CM,
296  typename AS, typename AL, typename AD, typename AM,
297  typename BS, typename BL, typename BD, typename BM>
298 Kokkos::View< Sacado::UQ::PCE<CS>*, CL, CD, CM>
299 V_ElementWiseMultiply(
300  const Sacado::UQ::PCE<CS>& c,
301  const Kokkos::View< Sacado::UQ::PCE<CS>*, CL, CD, CM >& C,
302  const Sacado::UQ::PCE<AS>& ab,
303  const Kokkos::View< Sacado::UQ::PCE<AS>*, AL, AD, AM >& A,
304  const Kokkos::View< Sacado::UQ::PCE<BS>*, BL, BD, BM >& B )
305 {
306  if (Sacado::is_constant(c) && Sacado::is_constant(ab)) {
307  return V_ElementWiseMultiply( c.fastAccessCoeff(0), C,
308  ab.fastAccessCoeff(0), A, B );
309  }
310  else {
311  Impl::raise_error(
312  "V_ElementWiseMultiply not implemented for non-constant c or ab");
313  }
314  return C;
315 }
316 
317 // Rank-2 vector multiply with Sacado::UQ::PCE scalar type, constant c, ab
318 template <typename CS, typename CL, typename CD, typename CM,
319  typename AS, typename AL, typename AD, typename AM,
320  typename BS, typename BL, typename BD, typename BM>
321 Kokkos::View< Sacado::UQ::PCE<CS>**, CL, CD, CM>
322 MV_ElementWiseMultiply(
323  const typename Sacado::UQ::PCE<CS>::value_type& c,
324  const Kokkos::View< Sacado::UQ::PCE<CS>**, CL, CD, CM >& C,
325  const typename Sacado::UQ::PCE<AS>::value_type& ab,
326  const Kokkos::View< Sacado::UQ::PCE<AS>*, AL, AD, AM >& A,
327  const Kokkos::View< Sacado::UQ::PCE<BS>**, BL, BD, BM >& B )
328 {
329  typedef View< Sacado::UQ::PCE<CS>**, CL, CD, CM > CVector;
330  typedef View< Sacado::UQ::PCE<AS>*, AL, AD, AM > AVector;
331  typedef View< Sacado::UQ::PCE<BS>**, BL, BD, BM > BVector;
332 
333  typedef View< typename CVector::data_type, typename CVector::array_layout, typename CVector::execution_space, typename CVector::memory_traits > CView;
334  typedef View< typename AVector::data_type, typename AVector::array_layout, typename AVector::execution_space, typename AVector::memory_traits > AView;
335  typedef View< typename BVector::data_type, typename BVector::array_layout, typename BVector::execution_space, typename BVector::memory_traits > BView;
336 
337  MV_ElementWiseMultiplyFunctor<CView,AView,BView> op(c,C,ab,A,B,C.extent(1)) ;
338  Kokkos::parallel_for( C.extent(0) , op );
339  return C;
340 }
341 
342 // Rank-2 vector multiply with Sacado::UQ::PCE scalar type, non-constant c, ab
343 template <typename CS, typename CL, typename CD, typename CM,
344  typename AS, typename AL, typename AD, typename AM,
345  typename BS, typename BL, typename BD, typename BM>
346 Kokkos::View< Sacado::UQ::PCE<CS>**, CL, CD, CM>
347 MV_ElementWiseMultiply(
348  const Sacado::UQ::PCE<CS>& c,
349  const Kokkos::View< Sacado::UQ::PCE<CS>**, CL, CD, CM >& C,
350  const Sacado::UQ::PCE<AS>& ab,
351  const Kokkos::View< Sacado::UQ::PCE<AS>*, AL, AD, AM >& A,
352  const Kokkos::View< Sacado::UQ::PCE<BS>**, BL, BD, BM >& B )
353 {
354  if (Sacado::is_constant(c) && Sacado::is_constant(ab)) {
355  return MV_ElementWiseMultiply( c.fastAccessCoeff(0), C,
356  ab.fastAccessCoeff(0), A, B );
357  }
358  else {
359  Impl::raise_error(
360  "MV_ElementWiseMultiply not implemented for non-constant c or ab");
361  }
362  return C;
363 }
364 
365 // Rank-1 vector scale with Sacado::UQ::PCE scalar type, constant a, b
366 template <typename RS, typename RL, typename RD, typename RM,
367  typename XS, typename XL, typename XD, typename XM>
368 Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM>
369 V_MulScalar( const Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM >& r,
370  const typename Sacado::UQ::PCE<XS>::value_type& a,
371  const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x )
372 {
373  typedef Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM > RVector;
374  typedef Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM > XVector;
375 
376  typename RVector::flat_array_type r_flat = r;
377  typename XVector::flat_array_type x_flat = x;
378 
379  V_MulScalar( r_flat, a, x_flat );
380 
381  return r;
382 }
383 
384 // Rank-1 vector scale with Sacado::UQ::PCE scalar type, non-constant a, b
385 template <typename RS, typename RL, typename RD, typename RM,
386  typename XS, typename XL, typename XD, typename XM>
387 Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM>
388 V_MulScalar( const Kokkos::View< Sacado::UQ::PCE<RS>*, RL, RD, RM >& r,
389  const Sacado::UQ::PCE<XS>& a,
390  const Kokkos::View< Sacado::UQ::PCE<XS>*, XL, XD, XM >& x )
391 {
392  if (Sacado::is_constant(a)) {
393  return V_MulScalar( r, a.fastAccessCoeff(0), x );
394  }
395  else {
396  Impl::raise_error("V_MulScalar not implemented for non-constant a");
397  }
398  return r;
399 }
400 
401 // Rank-2 vector scale with Sacado::UQ::PCE scalar type, constant a, b
402 template <typename RS, typename RL, typename RD, typename RM,
403  typename XS, typename XL, typename XD, typename XM>
404 Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM>
405 MV_MulScalar( const Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM >& r,
406  const typename Sacado::UQ::PCE<XS>::value_type& a,
407  const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x )
408 {
409  typedef Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM > RVector;
410  typedef Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM > XVector;
411 
412  typename RVector::flat_array_type r_flat = r;
413  typename XVector::flat_array_type x_flat = x;
414 
415  MV_MulScalar( r_flat, a, x_flat );
416 
417  return r;
418 }
419 
420 // Rank-2 vector scale with Sacado::UQ::PCE scalar type, non-constant a, b
421 template <typename RS, typename RL, typename RD, typename RM,
422  typename XS, typename XL, typename XD, typename XM>
423 Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM>
424 MV_MulScalar( const Kokkos::View< Sacado::UQ::PCE<RS>**, RL, RD, RM >& r,
425  const Sacado::UQ::PCE<XS>& a,
426  const Kokkos::View< Sacado::UQ::PCE<XS>**, XL, XD, XM >& x )
427 {
428  if (Sacado::is_constant(a)) {
429  return MV_MulScalar( r, a.fastAccessCoeff(0), x );
430  }
431  else {
432  Impl::raise_error("MV_MulScalar not implemented for non-constant a or b");
433  }
434  return r;
435 }
436 
437 template <typename T>
438 struct V_ReciprocalThresholdSelfFunctor;
439 
440 template <typename T, typename D, typename M>
441 struct V_ReciprocalThresholdSelfFunctor<
442  View< T,LayoutLeft,D,M,Impl::ViewPCEContiguous > >
443 {
444  typedef View< T,LayoutLeft,D,M,Impl::ViewPCEContiguous > XVector;
445  typedef typename XVector::array_type array_type;
446 
447  typedef typename array_type::execution_space execution_space;
448  typedef typename array_type::size_type size_type;
449  typedef typename array_type::non_const_value_type value_type;
450  typedef Kokkos::ArithTraits<value_type> KAT;
451  typedef typename KAT::mag_type mag_type;
452 
453  const array_type m_x;
454  const value_type m_min_val;
455  const value_type m_min_val_mag;
456  const size_type m_n_pce;
457 
458  V_ReciprocalThresholdSelfFunctor(
459  const XVector& x,
460  const typename XVector::non_const_value_type& min_val) :
461  m_x(x),
462  m_min_val(min_val.fastAccessCoeff(0)),
463  m_min_val_mag(KAT::abs(m_min_val)),
464  m_n_pce(x.sacado_size()) {}
465  //--------------------------------------------------------------------------
466 
467  KOKKOS_INLINE_FUNCTION
468  void operator()( const size_type i) const
469  {
470  if (KAT::abs(m_x(0,i)) < m_min_val_mag)
471  m_x(0,i) = m_min_val;
472  else
473  m_x(0,i) = KAT::one() / m_x(0,i);
474  for (size_type l=1; l<m_n_pce; ++l)
475  m_x(l,i) = KAT::zero();
476  }
477 };
478 
479 template <typename T>
480 struct MV_ReciprocalThresholdSelfFunctor;
481 
482 template <typename T, typename D, typename M>
483 struct MV_ReciprocalThresholdSelfFunctor<
484  View< T,LayoutLeft,D,M,Impl::ViewPCEContiguous > >
485 {
486  typedef View< T,LayoutLeft,D,M,Impl::ViewPCEContiguous > XVector;
487  typedef typename XVector::array_type array_type;
488 
489  typedef typename array_type::execution_space execution_space;
490  typedef typename array_type::size_type size_type;
491  typedef typename array_type::non_const_value_type value_type;
492  typedef Kokkos::ArithTraits<value_type> KAT;
493  typedef typename KAT::mag_type mag_type;
494 
495  const array_type m_x;
496  const value_type m_min_val;
497  const value_type m_min_val_mag;
498  const size_type m_n;
499  const size_type m_n_pce;
500 
501  MV_ReciprocalThresholdSelfFunctor(
502  const XVector& x,
503  const typename XVector::non_const_value_type& min_val,
504  const size_type& n) :
505  m_x(x),
506  m_min_val(min_val.fastAccessCoeff(0)),
507  m_min_val_mag(KAT::abs(m_min_val)),
508  m_n(n),
509  m_n_pce(x.sacado_size()) {}
510  //--------------------------------------------------------------------------
511 
512  KOKKOS_INLINE_FUNCTION
513  void operator()( const size_type i) const
514  {
515  for (size_type k=0; k<m_n; ++k) {
516  if (KAT::abs(m_x(0,i,k)) < m_min_val_mag)
517  m_x(0,i,k) = m_min_val;
518  else
519  m_x(0,i,k) = KAT::one() / m_x(0,i,k);
520  for (size_type l=1; l<m_n_pce; ++l)
521  m_x(l,i,k) = KAT::zero();
522  }
523  }
524 };
525 
526 } // namespace Kokkos
527 */
528 #endif /* #ifndef KOKKOS_MV_UQ_PCE_HPP */