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 
67  typedef typename Teuchos::ScalarTraits<value_type>::magnitudeType value_mag_type;
68  typedef typename Teuchos::ScalarTraits<value_type>::halfPrecision value_half_type;
69  typedef typename Teuchos::ScalarTraits<value_type>::doublePrecision value_double_type;
70 
71  typedef typename Sacado::mpl::apply<storage_type,ordinal_type,value_mag_type>::type storage_mag_type;
72  typedef typename Sacado::mpl::apply<storage_type,ordinal_type,value_half_type>::type storage_half_type;
73  typedef typename Sacado::mpl::apply<storage_type,ordinal_type,value_double_type>::type storage_double_type;
74 
75  typedef typename Sacado::mpl::apply<PCEType, storage_mag_type>::type magnitudeType;
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 = Teuchos::ScalarTraits<value_type>::isComplex;
83  static const bool isOrdinal = Teuchos::ScalarTraits<value_type>::isOrdinal;
84  static const bool isComparable =
85  Teuchos::ScalarTraits<value_type>::isComparable;
86  static const bool hasMachineParameters =
87  Teuchos::ScalarTraits<value_type>::hasMachineParameters;
88 
89  static typename Teuchos::ScalarTraits<value_type>::magnitudeType eps() {
90  return Teuchos::ScalarTraits<value_type>::eps();
91  }
92 
93  static typename Teuchos::ScalarTraits<value_type>::magnitudeType sfmin() {
94  return Teuchos::ScalarTraits<value_type>::sfmin();
95  }
96 
97  static typename Teuchos::ScalarTraits<value_type>::magnitudeType base() {
98  return Teuchos::ScalarTraits<value_type>::base();
99  }
100 
101  static typename Teuchos::ScalarTraits<value_type>::magnitudeType prec() {
102  return Teuchos::ScalarTraits<value_type>::prec();
103  }
104 
105  static typename Teuchos::ScalarTraits<value_type>::magnitudeType t() {
106  return Teuchos::ScalarTraits<value_type>::t();
107  }
108 
109  static typename Teuchos::ScalarTraits<value_type>::magnitudeType rnd() {
111  }
112 
113  static typename Teuchos::ScalarTraits<value_type>::magnitudeType emin() {
114  return Teuchos::ScalarTraits<value_type>::emin();
115  }
116 
117  static typename Teuchos::ScalarTraits<value_type>::magnitudeType rmin() {
118  return Teuchos::ScalarTraits<value_type>::rmin();
119  }
120 
121  static typename Teuchos::ScalarTraits<value_type>::magnitudeType emax() {
122  return Teuchos::ScalarTraits<value_type>::emax();
123  }
124 
125  static typename Teuchos::ScalarTraits<value_type>::magnitudeType rmax() {
126  return Teuchos::ScalarTraits<value_type>::rmax();
127  }
128 
129  static magnitudeType magnitude(const PCEType& a) {
130  return a.two_norm();
131  }
132 
133  static innerProductType innerProduct(const PCEType& a, const PCEType& b) {
134  return a.inner_product(b);
135  }
136 
137  static PCEType zero() {
138  return PCEType(0.0);
139  }
140 
141  static PCEType one() {
142  return PCEType(1.0);
143  }
144 
145  // Conjugate is only defined for real derivative components
146 
147  static PCEType conjugate(const PCEType& x) {
148  PCEType y = x;
149  y.copyForWrite();
150  y.val() = Teuchos::ScalarTraits<value_type>::conjugate(x.val());
151  return y;
152  }
153 
154  // Real part is only defined for real derivative components
155 
156  static PCEType real(const PCEType& x) {
157  PCEType y = x;
158  y.copyForWrite();
159  y.val() = Teuchos::ScalarTraits<value_type>::real(x.val());
160  return y;
161  }
162 
163  // Imaginary part is only defined for real derivative components
164 
165  static PCEType imag(const PCEType& x) {
166  return PCEType(Teuchos::ScalarTraits<value_type>::imag(x.val()));
167  }
168 
169 
170  static value_type nan() {
171  return Teuchos::ScalarTraits<value_type>::nan();
172  }
173 
174  static bool isnaninf(const PCEType& x) {
175  for (int i=0; i<x.size(); i++)
176  if (Teuchos::ScalarTraits<value_type>::isnaninf(x.fastAccessCoeff(i)))
177  return true;
178  return false;
179  }
180 
181  static void seedrandom(unsigned int s) {
182  Teuchos::ScalarTraits<value_type>::seedrandom(s);
183  }
184 
185  static value_type random() {
186  return Teuchos::ScalarTraits<value_type>::random();
187  }
188 
189  static const char * name() {
190  return "Sacado::UQ::PCE<>";
191  }
192 
193  static PCEType squareroot(const PCEType& x) {
194  return std::sqrt(x);
195  }
196 
197  static PCEType pow(const PCEType& x, const PCEType& y) {
198  return std::pow(x,y);
199  }
200 
201  static PCEType log(const PCEType& x) {
202  return std::log(x);
203  }
204 
205  static PCEType log10(const PCEType& x) {
206  return std::log10(x);
207  }
208 
209  // Helper function to determine whether a complex value is real
210 
211  static bool is_complex_real(const value_type& x) {
212  return
213  Teuchos::ScalarTraits<value_type>::magnitude(x-Teuchos::ScalarTraits<value_type>::real(x)) == 0;
214  }
215 
216  // Helper function to determine whether a Fad type is real
217 
218  static bool is_pce_real(const PCEType& x) {
219  if (x.size() == 0)
220  return true;
221  if (Teuchos::ScalarTraits<value_type>::isComplex) {
222  for (int i=0; i<x.size(); i++)
223  if (!is_complex_real(x.fastAccessCoeff(i)))
224  return false;
225  }
226  return true;
227  }
228 
229  }; // class PCEScalarTraitsImp
230 
232  template <typename TypeTo, typename PCEType>
233  struct PCEValueTypeConversionTraitsImp {
234  typedef typename Sacado::ValueType<PCEType>::type ValueT;
235  typedef Teuchos::ValueTypeConversionTraits<TypeTo,ValueT> VTCT;
236  static TypeTo convert( const PCEType t ) {
237  return VTCT::convert(t.val());
238  }
239  static TypeTo safeConvert( const PCEType t ) {
240  return VTCT::safeConvert(t.val());
241  }
242  };
243 
245  template <typename Ordinal, typename PCEType>
246  class PCESerializationTraitsImp {
247  typedef typename Sacado::ValueType<PCEType>::type ValueT;
248  typedef Teuchos::SerializationTraits<Ordinal,ValueT> vSerT;
249  typedef Teuchos::SerializationTraits<Ordinal,int> iSerT;
250  typedef Teuchos::SerializationTraits<Ordinal,Ordinal> oSerT;
251 
252  public:
253 
255  static const bool supportsDirectSerialization = false;
256 
258 
259 
261  static Ordinal fromCountToIndirectBytes(const Ordinal count,
262  const PCEType buffer[]) {
263  Ordinal bytes = 0;
264  for (Ordinal i=0; i<count; i++) {
265  int sz = buffer[i].size();
266  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &sz);
267  Ordinal b2 = vSerT::fromCountToIndirectBytes(sz, buffer[i].coeff());
268  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
269  bytes += b1+b2+b3;
270  }
271  return bytes;
272  }
273 
275  static void serialize (const Ordinal count,
276  const PCEType buffer[],
277  const Ordinal bytes,
278  char charBuffer[]) {
279  for (Ordinal i=0; i<count; i++) {
280  // First serialize size
281  int sz = buffer[i].size();
282  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &sz);
283  iSerT::serialize(1, &sz, b1, charBuffer);
284  charBuffer += b1;
285 
286  // Next serialize PCE coefficients
287  Ordinal b2 = vSerT::fromCountToIndirectBytes(sz, buffer[i].coeff());
288  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
289  oSerT::serialize(1, &b2, b3, charBuffer);
290  charBuffer += b3;
291  vSerT::serialize(sz, buffer[i].coeff(), b2, charBuffer);
292  charBuffer += b2;
293  }
294  }
295 
297  static Ordinal fromIndirectBytesToCount(const Ordinal bytes,
298  const char charBuffer[]) {
299  Ordinal count = 0;
300  Ordinal bytes_used = 0;
301  while (bytes_used < bytes) {
302 
303  // Bytes for size
304  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
305  bytes_used += b1;
306  charBuffer += b1;
307 
308  // Bytes for PCE coefficients
309  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
310  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
311  bytes_used += b3;
312  charBuffer += b3;
313  bytes_used += *b2;
314  charBuffer += *b2;
315 
316  ++count;
317  }
318  return count;
319  }
320 
322  static void deserialize (const Ordinal bytes,
323  const char charBuffer[],
324  const Ordinal count,
325  PCEType buffer[]) {
326  for (Ordinal i=0; i<count; i++) {
327 
328  // Deserialize size
329  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
330  const int *sz = iSerT::convertFromCharPtr(charBuffer);
331  charBuffer += b1;
332 
333  // Make sure PCE object is ready to receive values
334  // We assume it has already been initialized with the proper
335  // cijk object
336  if (buffer[i].size() != *sz)
337  buffer[i].reset(buffer[i].cijk(), *sz);
338  buffer[i].copyForWrite();
339 
340  // Deserialize PCE coefficients
341  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
342  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
343  charBuffer += b3;
344  vSerT::deserialize(*b2, charBuffer, *sz, buffer[i].coeff());
345  charBuffer += *b2;
346  }
347 
348  }
349 
351 
352  };
353 
354 
356  template <typename Ordinal, typename PCEType, typename ValueSerializer>
357  class PCESerializerImp {
358 
359  public:
360 
362  typedef ValueSerializer value_serializer_type;
363 
365  typedef typename PCEType::cijk_type cijk_type;
366 
367 
368  protected:
369  typedef typename Sacado::ValueType<PCEType>::type ValueT;
370  typedef Teuchos::SerializationTraits<Ordinal,int> iSerT;
371  typedef Teuchos::SerializationTraits<Ordinal,Ordinal> oSerT;
372 
373  cijk_type cijk;
374  Teuchos::RCP<const ValueSerializer> vs;
375  int sz;
376 
377  public:
378 
380  static const bool supportsDirectSerialization = false;
381 
382  PCESerializerImp(const cijk_type& cijk_,
383  const Teuchos::RCP<const ValueSerializer>& vs_) :
384  cijk(cijk_), vs(vs_), sz(cijk.dimension()) {}
385 
387  cijk_type getSerializerCijk() const { return cijk; }
388 
390  Teuchos::RCP<const value_serializer_type> getValueSerializer() const {
391  return vs; }
392 
394 
395 
397  Ordinal fromCountToIndirectBytes(const Ordinal count,
398  const PCEType buffer[]) const {
399  Ordinal bytes = 0;
400  PCEType *x = NULL;
401  const PCEType *cx;
402  for (Ordinal i=0; i<count; i++) {
403  int my_sz = buffer[i].size();
404  if (sz != my_sz) {
405  if (x == NULL)
406  x = new PCEType;
407  *x = buffer[i];
408  x->reset(cijk);
409  cx = x;
410  }
411  else
412  cx = &(buffer[i]);
413  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &sz);
414  Ordinal b2 = vs->fromCountToIndirectBytes(sz, cx->coeff());
415  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
416  bytes += b1+b2+b3;
417  }
418  if (x != NULL)
419  delete x;
420  return bytes;
421  }
422 
424  void serialize (const Ordinal count,
425  const PCEType buffer[],
426  const Ordinal bytes,
427  char charBuffer[]) const {
428  PCEType *x = NULL;
429  const PCEType *cx;
430  for (Ordinal i=0; i<count; i++) {
431  // First serialize size
432  int my_sz = buffer[i].size();
433  if (sz != my_sz) {
434  if (x == NULL)
435  x = new PCEType(cijk);
436  *x = buffer[i];
437  x->reset(cijk);
438  cx = x;
439  }
440  else
441  cx = &(buffer[i]);
442  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &sz);
443  iSerT::serialize(1, &sz, b1, charBuffer);
444  charBuffer += b1;
445 
446  // Next serialize PCE coefficients
447  Ordinal b2 = vs->fromCountToIndirectBytes(sz, cx->coeff());
448  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
449  oSerT::serialize(1, &b2, b3, charBuffer);
450  charBuffer += b3;
451  vs->serialize(sz, cx->coeff(), b2, charBuffer);
452  charBuffer += b2;
453  }
454  if (x != NULL)
455  delete x;
456  }
457 
459  Ordinal fromIndirectBytesToCount(const Ordinal bytes,
460  const char charBuffer[]) const {
461  Ordinal count = 0;
462  Ordinal bytes_used = 0;
463  while (bytes_used < bytes) {
464 
465  // Bytes for size
466  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
467  bytes_used += b1;
468  charBuffer += b1;
469 
470  // Bytes for PCE coefficients
471  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
472  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
473  bytes_used += b3;
474  charBuffer += b3;
475  bytes_used += *b2;
476  charBuffer += *b2;
477 
478  ++count;
479  }
480  return count;
481  }
482 
484  void deserialize (const Ordinal bytes,
485  const char charBuffer[],
486  const Ordinal count,
487  PCEType buffer[]) const {
488  for (Ordinal i=0; i<count; i++) {
489 
490  // Deserialize size
491  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
492  const int *my_sz = iSerT::convertFromCharPtr(charBuffer);
493  charBuffer += b1;
494 
495  // Create empty PCE object of given size
496  buffer[i].reset(cijk);
497 
498  // Deserialize PCE coefficients
499  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
500  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
501  charBuffer += b3;
502  vs->deserialize(*b2, charBuffer, *my_sz, buffer[i].coeff());
503  charBuffer += *b2;
504  }
505 
506  }
507 
509 
510  };
511 
512  } // namespace UQ
513 
514 } // namespace Sacado
515 
516 #endif // HAVE_SACADO_TEUCHOS
517 
518 #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