Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Sacado_UQ_PCE_ScalarTraitsImp.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 SACADO_UQ_PCE_SCALARTRAITSIMP_HPP
43 #define SACADO_UQ_PCE_SCALARTRAITSIMP_HPP
44 
45 #ifdef HAVE_SACADO_TEUCHOS
46 
47 #include "Teuchos_ScalarTraits.hpp"
48 #include "Teuchos_SerializationTraits.hpp"
49 #include "Teuchos_RCP.hpp"
50 #include "Teuchos_Assert.hpp"
51 #include "Sacado_mpl_apply.hpp"
52 #include "Teuchos_as.hpp"
53 
54 #include <iterator>
55 
56 namespace Sacado {
57 
58  namespace UQ {
59 
61  template <typename PCEType>
62  struct PCEScalarTraitsImp {
63  typedef typename PCEType::storage_type storage_type;
64  typedef typename storage_type::value_type value_type;
66  typedef Teuchos::ScalarTraits<value_type> TVT;
67 
68  typedef typename TVT::magnitudeType value_mag_type;
69  typedef typename TVT::halfPrecision value_half_type;
70  typedef typename TVT::doublePrecision value_double_type;
71 
72  typedef typename Sacado::mpl::apply<storage_type,ordinal_type,value_mag_type>::type storage_mag_type;
73  typedef typename Sacado::mpl::apply<storage_type,ordinal_type,value_half_type>::type storage_half_type;
74  typedef typename Sacado::mpl::apply<storage_type,ordinal_type,value_double_type>::type storage_double_type;
75 
76  typedef value_mag_type magnitudeType;
77  typedef typename Sacado::mpl::apply<PCEType, storage_half_type>::type halfPrecision;
78  typedef typename Sacado::mpl::apply<PCEType, storage_double_type>::type doublePrecision;
79 
80  typedef value_type innerProductType;
81 
82  static const bool isComplex = TVT::isComplex;
83  static const bool isOrdinal = TVT::isOrdinal;
84  static const bool isComparable = TVT::isComparable;
85  static const bool hasMachineParameters = TVT::hasMachineParameters;
86 
87  static value_mag_type eps() { return TVT::eps(); }
88 
89  static value_mag_type sfmin() { return TVT::sfmin(); }
90 
91  static value_mag_type base() { return TVT::base(); }
92 
93  static value_mag_type prec() { return TVT::prec(); }
94 
95  static value_mag_type t() { return TVT::t(); }
96 
97  static value_mag_type rnd() { return TVT::rnd(); }
98 
99  static value_mag_type emin() { return TVT::emin(); }
100 
101  static value_mag_type rmin() { return TVT::rmin(); }
102 
103  static value_mag_type emax() { return TVT::emax(); }
104 
105  static value_mag_type rmax() { return TVT::rmax(); }
106 
107  static magnitudeType magnitude(const PCEType& a) {
108  return a.two_norm();
109  }
110 
111  static innerProductType innerProduct(const PCEType& a, const PCEType& b) {
112  return a.inner_product(b);
113  }
114 
115  static PCEType zero() { return PCEType(0.0); }
116 
117  static PCEType one() { return PCEType(1.0); }
118 
119  static PCEType conjugate(const PCEType& x) {
120  PCEType y = x;
121  y.copyForWrite();
122  y.val() = TVT::conjugate(x.val());
123  return y;
124  }
125 
126  static magnitudeType real(const PCEType& x) {
127  magnitudeType m = magnitudeType(0.0);
128  const ordinal_type sz = x.size();
129  for (ordinal_type i=0; i<sz; ++i) {
130  value_mag_type t = TVT::real(x.fastAccessCoeff(i));
131  m +=t*t;
132  }
133  return std::sqrt(m);
134  }
135 
136  static magnitudeType imag(const PCEType& x) {
137  magnitudeType m = magnitudeType(0.0);
138  const ordinal_type sz = x.size();
139  for (ordinal_type i=0; i<sz; ++i) {
140  value_mag_type t = TVT::imag(x.fastAccessCoeff(i));
141  m +=t*t;
142  }
143  return std::sqrt(m);
144  }
145 
146 
147  static value_type nan() { return TVT::nan(); }
148 
149  static bool isnaninf(const PCEType& x) {
150  for (int i=0; i<x.size(); i++)
151  if (TVT::isnaninf(x.fastAccessCoeff(i)))
152  return true;
153  return false;
154  }
155 
156  static void seedrandom(unsigned int s) { TVT::seedrandom(s); }
157 
158  static PCEType random() { return PCEType(TVT::random()); }
159 
160  static const char * name() { return "Sacado::UQ::PCE<>"; }
161 
162  static PCEType squareroot(const PCEType& x) { return std::sqrt(x); }
163 
164  static PCEType pow(const PCEType& x, const PCEType& y) {
165  return std::pow(x,y);
166  }
167 
168  static PCEType log(const PCEType& x) { return std::log(x); }
169 
170  static PCEType log10(const PCEType& x) { return std::log10(x); }
171 
172  }; // class PCEScalarTraitsImp
173 
175  template <typename TypeTo, typename PCEType>
176  struct PCEValueTypeConversionTraitsImp {
177  typedef typename Sacado::ValueType<PCEType>::type ValueT;
178  typedef Teuchos::ValueTypeConversionTraits<TypeTo,ValueT> VTCT;
179  static TypeTo convert( const PCEType t ) {
180  return VTCT::convert(t.val());
181  }
182  static TypeTo safeConvert( const PCEType t ) {
183  return VTCT::safeConvert(t.val());
184  }
185  };
186 
188  template <typename Ordinal, typename PCEType>
189  class PCESerializationTraitsImp {
190  typedef typename Sacado::ValueType<PCEType>::type ValueT;
191  typedef Teuchos::SerializationTraits<Ordinal,ValueT> vSerT;
192  typedef Teuchos::SerializationTraits<Ordinal,int> iSerT;
193  typedef Teuchos::SerializationTraits<Ordinal,Ordinal> oSerT;
194 
195  public:
196 
198  static const bool supportsDirectSerialization = false;
199 
201 
202 
204  static Ordinal fromCountToIndirectBytes(const Ordinal count,
205  const PCEType buffer[]) {
206  Ordinal bytes = 0;
207  for (Ordinal i=0; i<count; i++) {
208  int sz = buffer[i].size();
209  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &sz);
210  Ordinal b2 = vSerT::fromCountToIndirectBytes(sz, buffer[i].coeff());
211  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
212  bytes += b1+b2+b3;
213  }
214  return bytes;
215  }
216 
218  static void serialize (const Ordinal count,
219  const PCEType buffer[],
220  const Ordinal bytes,
221  char charBuffer[]) {
222  for (Ordinal i=0; i<count; i++) {
223  // First serialize size
224  int sz = buffer[i].size();
225  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &sz);
226  iSerT::serialize(1, &sz, b1, charBuffer);
227  charBuffer += b1;
228 
229  // Next serialize PCE coefficients
230  Ordinal b2 = vSerT::fromCountToIndirectBytes(sz, buffer[i].coeff());
231  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
232  oSerT::serialize(1, &b2, b3, charBuffer);
233  charBuffer += b3;
234  vSerT::serialize(sz, buffer[i].coeff(), b2, charBuffer);
235  charBuffer += b2;
236  }
237  }
238 
240  static Ordinal fromIndirectBytesToCount(const Ordinal bytes,
241  const char charBuffer[]) {
242  Ordinal count = 0;
243  Ordinal bytes_used = 0;
244  while (bytes_used < bytes) {
245 
246  // Bytes for size
247  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
248  bytes_used += b1;
249  charBuffer += b1;
250 
251  // Bytes for PCE coefficients
252  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
253  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
254  bytes_used += b3;
255  charBuffer += b3;
256  bytes_used += *b2;
257  charBuffer += *b2;
258 
259  ++count;
260  }
261  return count;
262  }
263 
265  static void deserialize (const Ordinal bytes,
266  const char charBuffer[],
267  const Ordinal count,
268  PCEType buffer[]) {
269  for (Ordinal i=0; i<count; i++) {
270 
271  // Deserialize size
272  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
273  const int *sz = iSerT::convertFromCharPtr(charBuffer);
274  charBuffer += b1;
275 
276  // Make sure PCE object is ready to receive values
277  // We assume it has already been initialized with the proper
278  // cijk object
279  if (buffer[i].size() != *sz)
280  buffer[i].reset(buffer[i].cijk(), *sz);
281  buffer[i].copyForWrite();
282 
283  // Deserialize PCE coefficients
284  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
285  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
286  charBuffer += b3;
287  vSerT::deserialize(*b2, charBuffer, *sz, buffer[i].coeff());
288  charBuffer += *b2;
289  }
290 
291  }
292 
294 
295  };
296 
297 
299  template <typename Ordinal, typename PCEType, typename ValueSerializer>
300  class PCESerializerImp {
301 
302  public:
303 
305  typedef ValueSerializer value_serializer_type;
306 
308  typedef typename PCEType::cijk_type cijk_type;
309 
310 
311  protected:
312  typedef typename Sacado::ValueType<PCEType>::type ValueT;
313  typedef Teuchos::SerializationTraits<Ordinal,int> iSerT;
314  typedef Teuchos::SerializationTraits<Ordinal,Ordinal> oSerT;
315 
316  cijk_type cijk;
317  Teuchos::RCP<const ValueSerializer> vs;
318  int sz;
319 
320  public:
321 
323  static const bool supportsDirectSerialization = false;
324 
325  PCESerializerImp(const cijk_type& cijk_,
326  const Teuchos::RCP<const ValueSerializer>& vs_) :
327  cijk(cijk_), vs(vs_), sz(cijk.dimension()) {}
328 
330  cijk_type getSerializerCijk() const { return cijk; }
331 
333  Teuchos::RCP<const value_serializer_type> getValueSerializer() const {
334  return vs; }
335 
337 
338 
340  Ordinal fromCountToIndirectBytes(const Ordinal count,
341  const PCEType buffer[]) const {
342  Ordinal bytes = 0;
343  PCEType *x = NULL;
344  const PCEType *cx;
345  for (Ordinal i=0; i<count; i++) {
346  int my_sz = buffer[i].size();
347  if (sz != my_sz) {
348  if (x == NULL)
349  x = new PCEType;
350  *x = buffer[i];
351  x->reset(cijk);
352  cx = x;
353  }
354  else
355  cx = &(buffer[i]);
356  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &sz);
357  Ordinal b2 = vs->fromCountToIndirectBytes(sz, cx->coeff());
358  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
359  bytes += b1+b2+b3;
360  }
361  if (x != NULL)
362  delete x;
363  return bytes;
364  }
365 
367  void serialize (const Ordinal count,
368  const PCEType buffer[],
369  const Ordinal bytes,
370  char charBuffer[]) const {
371  PCEType *x = NULL;
372  const PCEType *cx;
373  for (Ordinal i=0; i<count; i++) {
374  // First serialize size
375  int my_sz = buffer[i].size();
376  if (sz != my_sz) {
377  if (x == NULL)
378  x = new PCEType(cijk);
379  *x = buffer[i];
380  x->reset(cijk);
381  cx = x;
382  }
383  else
384  cx = &(buffer[i]);
385  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &sz);
386  iSerT::serialize(1, &sz, b1, charBuffer);
387  charBuffer += b1;
388 
389  // Next serialize PCE coefficients
390  Ordinal b2 = vs->fromCountToIndirectBytes(sz, cx->coeff());
391  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
392  oSerT::serialize(1, &b2, b3, charBuffer);
393  charBuffer += b3;
394  vs->serialize(sz, cx->coeff(), b2, charBuffer);
395  charBuffer += b2;
396  }
397  if (x != NULL)
398  delete x;
399  }
400 
402  Ordinal fromIndirectBytesToCount(const Ordinal bytes,
403  const char charBuffer[]) const {
404  Ordinal count = 0;
405  Ordinal bytes_used = 0;
406  while (bytes_used < bytes) {
407 
408  // Bytes for size
409  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
410  bytes_used += b1;
411  charBuffer += b1;
412 
413  // Bytes for PCE coefficients
414  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
415  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
416  bytes_used += b3;
417  charBuffer += b3;
418  bytes_used += *b2;
419  charBuffer += *b2;
420 
421  ++count;
422  }
423  return count;
424  }
425 
427  void deserialize (const Ordinal bytes,
428  const char charBuffer[],
429  const Ordinal count,
430  PCEType buffer[]) const {
431  for (Ordinal i=0; i<count; i++) {
432 
433  // Deserialize size
434  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
435  const int *my_sz = iSerT::convertFromCharPtr(charBuffer);
436  charBuffer += b1;
437 
438  // Create empty PCE object of given size
439  buffer[i].reset(cijk);
440 
441  // Deserialize PCE coefficients
442  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
443  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
444  charBuffer += b3;
445  vs->deserialize(*b2, charBuffer, *my_sz, buffer[i].coeff());
446  charBuffer += *b2;
447  }
448 
449  }
450 
452 
453  };
454 
455  } // namespace UQ
456 
457 } // namespace Sacado
458 
459 #endif // HAVE_SACADO_TEUCHOS
460 
461 #endif // SACADO_FAD_SCALARTRAITSIMP_HPP
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Stokhos::StandardStorage< int, double > storage_type
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
Sacado::Random< double > rnd
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Sacado::UQ::PCE< storage_type > PCEType
KOKKOS_INLINE_FUNCTION PCE< Storage > log(const PCE< Storage > &a)
KOKKOS_INLINE_FUNCTION PCE< Storage > log10(const PCE< Storage > &a)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267