Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Stokhos_Tpetra_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 STOKHOS_TPETRA_UQ_PCE_HPP
43 #define STOKHOS_TPETRA_UQ_PCE_HPP
44 
45 // This header file should be included whenever compiling any Tpetra
46 // code with Stokhos scalar types
47 
48 // MP includes and specializations
50 
51 // Kokkos includes
52 #include "Tpetra_ConfigDefs.hpp"
53 #include "Kokkos_Core.hpp"
54 #include "Kokkos_BufferMacros.hpp"
55 #include "KokkosCompat_ClassicNodeAPI_Wrapper.hpp"
56 #include "KokkosCompat_View.hpp"
57 #include "KokkosCompat_View_def.hpp"
58 
59 // Hack for mean-based prec-setup where we get the PCE size from the first
60 // entry in the Teuchos::ArrayView. This is clearly quite dangerous and is
61 // likely to bite us in the ass at some point!
62 namespace Kokkos {
63  namespace Compat {
64  template <typename D, typename S>
65  Kokkos::View<Sacado::UQ::PCE<S>*,D>
66  getKokkosViewDeepCopy(const Teuchos::ArrayView< Sacado::UQ::PCE<S> >& a) {
67  typedef Sacado::UQ::PCE<S> T;
68  typedef typename Kokkos::Impl::if_c<
69  ::Kokkos::Impl::VerifyExecutionCanAccessMemorySpace< D, Kokkos::HostSpace>::value,
70  typename D::execution_space, Kokkos::HostSpace>::type
71  HostDevice;
72  typedef Kokkos::View<T*,D> view_type;
73  typedef Kokkos::View<T*,typename view_type::array_layout,HostDevice,Kokkos::MemoryUnmanaged> unmanaged_host_view_type;
74  if (a.size() == 0)
75  return view_type();
76  view_type v("", a.size(), a[0].size());
77  unmanaged_host_view_type hv(a.getRawPtr(), a.size(), a[0].size());
78  Kokkos::deep_copy(v,hv);
79  return v;
80  }
81 
82  template <typename D, typename S>
83  Kokkos::View<const Sacado::UQ::PCE<S>*,D>
84  getKokkosViewDeepCopy(const Teuchos::ArrayView<const Sacado::UQ::PCE<S> >& a) {
85  typedef Sacado::UQ::PCE<S> T;
86  typedef typename Kokkos::Impl::if_c<
87  ::Kokkos::Impl::VerifyExecutionCanAccessMemorySpace< D, Kokkos::HostSpace>::value,
88  typename D::execution_space, Kokkos::HostSpace>::type
89  HostDevice;
90  typedef Kokkos::View<T*,D> view_type;
91  typedef Kokkos::View<const T*,typename view_type::array_layout,HostDevice,Kokkos::MemoryUnmanaged> unmanaged_host_view_type;
92  if (a.size() == 0)
93  return view_type();
94  view_type v("", a.size(), a[0].size());
95  unmanaged_host_view_type hv(a.getRawPtr(), a.size(), a[0].size());
96  Kokkos::deep_copy(v,hv);
97  return v;
98  }
99  }
100 }
101 
102 // Kokkos-Linalg
105 #include "Kokkos_MV_UQ_PCE.hpp"
111 #include "Kokkos_Random_UQ_PCE.hpp"
112 
113 namespace Stokhos {
114 
115 // Traits for determining device type from node type
116 template <typename Node>
118  // Prefer Serial execution space as the default, but if that's not
119  // available, use the Host memory space's default execution space.
120 #if defined(KOKKOS_HAVE_SERIAL)
121  typedef Kokkos::Serial type;
122 #else
124 #endif // defined(KOKKOS_HAVE_SERIAL)
125 };
126 
127 template <typename Device>
128 struct DeviceForNode2< Kokkos::Compat::KokkosDeviceWrapperNode<Device> > {
129  typedef Device type;
130 };
131 
132 }
133 
134 #include "Tpetra_Details_PackTraits.hpp"
135 
136 namespace Tpetra {
137 namespace Details {
138 
143 template<typename S, typename D>
144 struct PackTraits< Sacado::UQ::PCE<S>, D > {
147  typedef D device_type;
148  typedef typename execution_space::size_type size_type;
149 
152  static const bool compileTimeSize = false;
153 
154  typedef Kokkos::View<const char*, device_type, Kokkos::MemoryUnmanaged> input_buffer_type;
155  typedef Kokkos::View<char*, device_type, Kokkos::MemoryUnmanaged> output_buffer_type;
156  typedef Kokkos::View<const value_type*, device_type, Kokkos::MemoryUnmanaged> input_array_type;
157  typedef Kokkos::View<value_type*, device_type, Kokkos::MemoryUnmanaged> output_array_type;
158 
160  typedef PackTraits< scalar_value_type, device_type > SPT;
161  typedef typename SPT::input_array_type scalar_input_array_type;
162  typedef typename SPT::output_array_type scalar_output_array_type;
163 
164  static size_t numValuesPerScalar (const value_type& x) {
165  return x.size ();
166  }
167 
168  static Kokkos::View<value_type*, device_type>
169  allocateArray (const value_type& x, const size_t numEnt, const std::string& label = "")
170  {
171  typedef Kokkos::View<value_type*, device_type> view_type;
172 
173  const size_type numVals = numValuesPerScalar (x);
174  return view_type (label, static_cast<size_type> (numEnt), numVals);
175  }
176 
177  static size_t
179  const input_array_type& inBuf,
180  const size_t numEnt)
181  {
182 #ifdef HAVE_TPETRA_DEBUG
183  TEUCHOS_TEST_FOR_EXCEPTION(
184  static_cast<size_t> (inBuf.dimension_0 ()) < numEnt,
185  std::invalid_argument, "PackTraits::packArray: inBuf.dimension_0() = "
186  << inBuf.dimension_0 () << " < numEnt = " << numEnt << ".");
187 #endif // HAVE_TPETRA_DEBUG
188 
189  if (numEnt == 0) {
190  return 0;
191  }
192  else {
193  // Check whether input array is contiguously allocated based on the size
194  // of the first entry. We can only pack contiguously allocated data
195  // since that is the only way we can guarrantee all of the PCE arrays
196  // are the same size and the buffer will allocated correctly.
197  const size_t scalar_size = numValuesPerScalar(inBuf(0));
198  const size_t in_dim = inBuf.dimension_0();
199  const scalar_value_type* last_coeff = inBuf(in_dim-1).coeff();
200  const scalar_value_type* last_coeff_expected =
201  inBuf(0).coeff() + (in_dim-1)*scalar_size;
202  const bool is_contiguous = (last_coeff == last_coeff_expected);
203  TEUCHOS_TEST_FOR_EXCEPTION(
204  !is_contiguous, std::logic_error,
205  "Cannot pack non-contiguous PCE array since buffer size calculation" <<
206  " is likely wrong.");
207 
208  // Check we are packing length-1 PCE arrays (mean-based preconditioner).
209  // We can technically pack length > 1, but the unpack assumes the
210  // output array is sized appropriately. Currently this is not the case
211  // in Tpetra::CrsMatrix::transferAndFillComplete() which allocates a
212  // local Teuchos::Array for the CSR values, which will only be length-1
213  // by default.
214  TEUCHOS_TEST_FOR_EXCEPTION(
215  scalar_size != 1, std::logic_error,
216  "Cannot pack PCE array with pce_size > 1 since unpack array" <<
217  " may not be allocated correctly.");
218 
219  const size_t flat_numEnt = numEnt * scalar_size;
220  scalar_input_array_type flat_inBuf(inBuf(0).coeff(), flat_numEnt);
221  return SPT::packArray(outBuf, flat_inBuf, flat_numEnt);
222  }
223  }
224 
225  static size_t
227  const input_buffer_type& inBuf,
228  const size_t numEnt)
229  {
230 #ifdef HAVE_TPETRA_DEBUG
231  TEUCHOS_TEST_FOR_EXCEPTION(
232  static_cast<size_t> (outBuf.dimension_0 ()) < numEnt,
233  std::invalid_argument,
234  "PackTraits::unpackArray: outBuf.dimension_0 () = " <<
235  outBuf.dimension_0 () << " < numEnt = " << numEnt << ".");
236 #endif // HAVE_TPETRA_DEBUG
237 
238  if (numEnt == 0) {
239  return static_cast<size_t> (0);
240  }
241  else {
242  // Check whether output array is contiguously allocated based on the size
243  // of the first entry. We have a simpler method to unpack in this case
244  const size_type scalar_size = numValuesPerScalar(outBuf(0));
245  const size_type out_dim = outBuf.dimension_0();
246  const scalar_value_type* last_coeff = outBuf(out_dim-1).coeff();
247  const scalar_value_type* last_coeff_expected =
248  outBuf(0).coeff() + (out_dim-1)*scalar_size;
249  const bool is_contiguous = (last_coeff == last_coeff_expected);
250 
251  if (is_contiguous) {
252  // Unpack all of the PCE coefficients for the whole array
253  const size_t flat_numEnt = numEnt * scalar_size;
254  scalar_output_array_type flat_outBuf(outBuf(0).coeff(), flat_numEnt);
255  return SPT::unpackArray(flat_outBuf, inBuf, flat_numEnt);
256  }
257  else {
258  // Unpack one entry at a time. This assumes each entry of outBuf
259  // is the correct size based on the packing. This is is only
260  // guarranteed to be true for pce_size == 1, hence the check in
261  // packArray().
262  size_t numBytesTotal = 0;
263  const size_type in_dim = inBuf.dimension_0();
264  for (size_t i=0; i<numEnt; ++i) {
265  input_buffer_type val_inBuf(inBuf.ptr_on_device()+numBytesTotal,
266  in_dim-numBytesTotal);
267  const size_t numBytes = unpackValue(outBuf(i), val_inBuf);
268  numBytesTotal += numBytes;
269  }
270  return numBytesTotal;
271  }
272  }
273  }
274 
275  static size_t
276  packValueCount (const value_type& inVal)
277  {
278  return inVal.size () * SPT::packValueCount (inVal.val ());
279  }
280 
281  static size_t
283  const value_type& inVal)
284  {
285  const size_t numBytes = packValueCount (inVal);
286  memcpy (outBuf.ptr_on_device (), inVal.coeff (), numBytes);
287  return numBytes;
288  }
289 
290  static size_t
291  unpackValue (value_type& outVal, const input_buffer_type& inBuf)
292  {
293  const size_t numBytes = packValueCount (outVal);
294  memcpy (outVal.coeff (), inBuf.ptr_on_device (), numBytes);
295  return numBytes;
296  }
297 }; // struct PackTraits
298 
299 } // namespace Details
300 } // namespace Tpetra
301 
302 namespace Tpetra {
303  template <class S, class L, class G, class N, bool> class MultiVector;
304  template <class S, class L, class G, class N, bool> class Vector;
305 }
306 
307 namespace Kokkos {
308  template <class S, class L, class G, class N, bool c>
311  typedef typename MV::dual_view_type dual_view_type;
312  typedef typename dual_view_type::t_dev device_type;
313  typedef typename dual_view_type::t_host host_type;
314  dual_view_type dual_view = mv.getDualView();
315  if (dual_view.modified_host() > dual_view.modified_device())
316  return dimension_scalar(dual_view.template view<device_type>());
317  return dimension_scalar(dual_view.template view<host_type>());
318  }
319 
320  template <class S, class L, class G, class N, bool c>
322  typedef Tpetra::Vector<S,L,G,N,c> V;
323  typedef typename V::dual_view_type dual_view_type;
324  typedef typename dual_view_type::t_dev device_type;
325  typedef typename dual_view_type::t_host host_type;
326  dual_view_type dual_view = v.getDualView();
327  if (dual_view.modified_host() > dual_view.modified_device())
328  return dimension_scalar(dual_view.template view<device_type>());
329  return dimension_scalar(dual_view.template view<host_type>());
330  }
331 }
332 
333 #endif // STOKHOS_TPETRA_UQ_PCE_HPP
Kokkos::View< const char *, device_type, Kokkos::MemoryUnmanaged > input_buffer_type
PackTraits< scalar_value_type, device_type > SPT
static size_t unpackArray(const output_array_type &outBuf, const input_buffer_type &inBuf, const size_t numEnt)
Kokkos::DefaultExecutionSpace execution_space
static size_t unpackValue(value_type &outVal, const input_buffer_type &inBuf)
Kokkos::View< Sacado::UQ::PCE< S > *, D > getKokkosViewDeepCopy(const Teuchos::ArrayView< Sacado::UQ::PCE< S > > &a)
Kokkos::HostSpace::execution_space type
static size_t packValue(const output_buffer_type &outBuf, const value_type &inVal)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
Kokkos::View< const value_type *, device_type, Kokkos::MemoryUnmanaged > input_array_type
Top-level namespace for Stokhos classes and functions.
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< View< T, P... > >::value, unsigned >::type dimension_scalar(const View< T, P... > &view)
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
static size_t packArray(const output_buffer_type &outBuf, const input_array_type &inBuf, const size_t numEnt)
static Kokkos::View< value_type *, device_type > allocateArray(const value_type &x, const size_t numEnt, const std::string &label="")
Kokkos::View< char *, device_type, Kokkos::MemoryUnmanaged > output_buffer_type
Kokkos::View< value_type *, device_type, Kokkos::MemoryUnmanaged > output_array_type