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