Sacado Package Browser (Single Doxygen Collection)  Version of the Day
Sacado_Tay_ScalarTraitsImp.hpp
Go to the documentation of this file.
1 // $Id$
2 // $Source$
3 // @HEADER
4 // ***********************************************************************
5 //
6 // Sacado Package
7 // Copyright (2006) Sandia Corporation
8 //
9 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
10 // the U.S. Government retains certain rights in this software.
11 //
12 // This library is free software; you can redistribute it and/or modify
13 // it under the terms of the GNU Lesser General Public License as
14 // published by the Free Software Foundation; either version 2.1 of the
15 // License, or (at your option) any later version.
16 //
17 // This library is distributed in the hope that it will be useful, but
18 // WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 // Lesser General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public
23 // License along with this library; if not, write to the Free Software
24 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
25 // USA
26 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
27 // (etphipp@sandia.gov).
28 //
29 // ***********************************************************************
30 // @HEADER
31 
32 #ifndef SACADO_TAY_SCALARTRAITSIMP_HPP
33 #define SACADO_TAY_SCALARTRAITSIMP_HPP
34 
35 #ifdef HAVE_SACADO_TEUCHOS
36 
37 #include "Teuchos_ScalarTraits.hpp"
38 #include "Teuchos_SerializationTraits.hpp"
39 #include "Teuchos_SerializationTraitsHelpers.hpp"
40 #include "Teuchos_Assert.hpp"
41 #include "Sacado_mpl_apply.hpp"
42 
43 #include <iterator>
44 
45 namespace Sacado {
46 
47  namespace Tay {
48 
50  template <typename TayType>
51  struct ScalarTraitsImp {
52  typedef typename Sacado::ValueType<TayType>::type ValueT;
53 
54  typedef typename mpl::apply<TayType,typename Teuchos::ScalarTraits<ValueT>::magnitudeType>::type magnitudeType;
55  typedef typename mpl::apply<TayType,typename Teuchos::ScalarTraits<ValueT>::halfPrecision>::type halfPrecision;
56  typedef typename mpl::apply<TayType,typename Teuchos::ScalarTraits<ValueT>::doublePrecision>::type doublePrecision;
57 
58  static const bool isComplex = Teuchos::ScalarTraits<ValueT>::isComplex;
59  static const bool isOrdinal = Teuchos::ScalarTraits<ValueT>::isOrdinal;
60  static const bool isComparable =
61  Teuchos::ScalarTraits<ValueT>::isComparable;
62  static const bool hasMachineParameters =
63  Teuchos::ScalarTraits<ValueT>::hasMachineParameters;
64  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType eps() {
65  return Teuchos::ScalarTraits<ValueT>::eps();
66  }
67  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType sfmin() {
68  return Teuchos::ScalarTraits<ValueT>::sfmin();
69  }
70  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType base() {
71  return Teuchos::ScalarTraits<ValueT>::base();
72  }
73  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType prec() {
74  return Teuchos::ScalarTraits<ValueT>::prec();
75  }
76  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType t() {
77  return Teuchos::ScalarTraits<ValueT>::t();
78  }
79  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType rnd() {
81  }
82  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType emin() {
83  return Teuchos::ScalarTraits<ValueT>::emin();
84  }
85  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType rmin() {
86  return Teuchos::ScalarTraits<ValueT>::rmin();
87  }
88  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType emax() {
89  return Teuchos::ScalarTraits<ValueT>::emax();
90  }
91  static typename Teuchos::ScalarTraits<ValueT>::magnitudeType rmax() {
92  return Teuchos::ScalarTraits<ValueT>::rmax();
93  }
94  static magnitudeType magnitude(const TayType& a) {
95 #ifdef TEUCHOS_DEBUG
96  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
97  a, "Error, the input value to magnitude(...) a = " << a <<
98  " can not be NaN!" );
99  TEUCHOS_TEST_FOR_EXCEPTION(is_tay_real(a) == false, std::runtime_error,
100  "Complex magnitude is not a differentiable "
101  "function of complex inputs.");
102 #endif
103  //return std::fabs(a);
104  magnitudeType b(a.degree(),
105  Teuchos::ScalarTraits<ValueT>::magnitude(a.val()));
106  if (Teuchos::ScalarTraits<ValueT>::real(a.val()) >= 0)
107  for (int i=1; i<=a.degree(); i++)
108  b.fastAccessCoeff(i) =
109  Teuchos::ScalarTraits<ValueT>::magnitude(a.fastAccessCoeff(i));
110  else
111  for (int i=1; i<=a.degree(); i++)
112  b.fastAccessCoeff(i) =
113  -Teuchos::ScalarTraits<ValueT>::magnitude(a.fastAccessCoeff(i));
114  return b;
115  }
116  static ValueT zero() {
117  return ValueT(0.0);
118  }
119  static ValueT one() {
120  return ValueT(1.0);
121  }
122 
123  // Conjugate is only defined for real derivative components
124  static TayType conjugate(const TayType& x) {
125 #ifdef TEUCHOS_DEBUG
126  TEUCHOS_TEST_FOR_EXCEPTION(is_tay_real(x) == false, std::runtime_error,
127  "Complex conjugate is not a differentiable "
128  "function of complex inputs.");
129 #endif
130  TayType y = x;
131  y.copyForWrite();
132  y.val() = Teuchos::ScalarTraits<ValueT>::conjugate(x.val());
133  return y;
134  }
135 
136  // Real part is only defined for real derivative components
137  static TayType real(const TayType& x) {
138 #ifdef TEUCHOS_DEBUG
139  TEUCHOS_TEST_FOR_EXCEPTION(is_tay_real(x) == false, std::runtime_error,
140  "Real component is not a differentiable "
141  "function of complex inputs.");
142 #endif
143  TayType y = x;
144  y.copyForWrite();
145  y.val() = Teuchos::ScalarTraits<ValueT>::real(x.val());
146  return y;
147  }
148 
149  // Imaginary part is only defined for real derivative components
150  static TayType imag(const TayType& x) {
151 #ifdef TEUCHOS_DEBUG
152  TEUCHOS_TEST_FOR_EXCEPTION(is_tay_real(x) == false, std::runtime_error,
153  "Imaginary component is not a differentiable "
154  "function of complex inputs.");
155 #endif
156  return TayType(Teuchos::ScalarTraits<ValueT>::imag(x.val()));
157  }
158 
159  static ValueT nan() {
160  return Teuchos::ScalarTraits<ValueT>::nan();
161  }
162  static bool isnaninf(const TayType& x) {
163  for (int i=0; i<=x.degree(); i++)
164  if (Teuchos::ScalarTraits<ValueT>::isnaninf(x.fastAccessCoeff(i)))
165  return true;
166  return false;
167  }
168  static void seedrandom(unsigned int s) {
169  Teuchos::ScalarTraits<ValueT>::seedrandom(s);
170  }
171  static ValueT random() {
172  return Teuchos::ScalarTraits<ValueT>::random();
173  }
174  static std::string name() {
176  }
177  static TayType squareroot(const TayType& x) {
178 #ifdef TEUCHOS_DEBUG
179  TEUCHOS_SCALAR_TRAITS_NAN_INF_ERR(
180  x, "Error, the input value to squareroot(...) a = " << x <<
181  " can not be NaN!" );
182 #endif
183  return std::sqrt(x);
184  }
185  static TayType pow(const TayType& x, const TayType& y) {
186  return std::pow(x,y);
187  }
188 
189  // Helper function to determine whether a complex value is real
190  static bool is_complex_real(const ValueT& x) {
191  return
192  Teuchos::ScalarTraits<ValueT>::magnitude(x-Teuchos::ScalarTraits<ValueT>::real(x)) == 0;
193  }
194 
195  // Helper function to determine whether a Fad type is real
196  static bool is_tay_real(const TayType& x) {
197  if (x.size() == 0)
198  return true;
199  if (Teuchos::ScalarTraits<ValueT>::isComplex) {
200  for (int i=0; i<=x.degree(); i++)
201  if (!is_complex_real(x.fastAccessCoeff(i)))
202  return false;
203  }
204  return true;
205  }
206 
207  }; // class ScalarTraitsImp
208 
210  template <typename Ordinal, typename TayType, typename Serializer>
211  struct SerializationImp {
212 
213  private:
214 
216  typedef Teuchos::SerializationTraits<Ordinal,unsigned int> iSerT;
217 
219  typedef Teuchos::SerializationTraits<Ordinal,Ordinal> oSerT;
220 
221  public:
222 
224  static const bool supportsDirectSerialization = false;
225 
227 
228 
230  static Ordinal fromCountToIndirectBytes(const Serializer& vs,
231  const Ordinal count,
232  const TayType buffer[],
233  const Ordinal sz = 0) {
234  Ordinal bytes = 0;
235  TayType *x = NULL;
236  const TayType *cx;
237  for (Ordinal i=0; i<count; i++) {
238  unsigned int my_sz = buffer[i].degree()+1;
239  unsigned int tot_sz = sz;
240  if (sz == 0) tot_sz = my_sz;
241  if (tot_sz != my_sz) {
242  x = new TayType(buffer[i]);
243  x->resize(tot_sz-1, true);
244  cx = x;
245  }
246  else
247  cx = &(buffer[i]);
248  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &tot_sz);
249  Ordinal b2 = vs.fromCountToIndirectBytes(
250  tot_sz, &(cx->fastAccessCoeff(0)));
251  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
252  bytes += b1+b2+b3;
253  if (x != NULL) {
254  delete x;
255  x = NULL;
256  }
257  }
258  return bytes;
259  }
260 
262  static void serialize (const Serializer& vs,
263  const Ordinal count,
264  const TayType buffer[],
265  const Ordinal bytes,
266  char charBuffer[],
267  const Ordinal sz = 0) {
268  TayType *x = NULL;
269  const TayType *cx;
270  for (Ordinal i=0; i<count; i++) {
271  // First serialize degree
272  unsigned int my_sz = buffer[i].degree()+1;
273  unsigned int tot_sz = sz;
274  if (sz == 0) tot_sz = my_sz;
275  Ordinal b1 = iSerT::fromCountToIndirectBytes(1, &tot_sz);
276  iSerT::serialize(1, &tot_sz, b1, charBuffer);
277  charBuffer += b1;
278 
279  // Next serialize taylor coefficients
280  if (tot_sz != my_sz) {
281  x = new TayType(buffer[i]);
282  x->resize(tot_sz-1, true);
283  cx = x;
284  }
285  else
286  cx = &(buffer[i]);
287  Ordinal b2 = vs.fromCountToIndirectBytes(
288  tot_sz, &(cx->fastAccessCoeff(0)));
289  Ordinal b3 = oSerT::fromCountToIndirectBytes(1, &b2);
290  oSerT::serialize(1, &b2, b3, charBuffer);
291  charBuffer += b3;
292  vs.serialize(tot_sz, &(cx->fastAccessCoeff(0)), b2, charBuffer);
293  charBuffer += b2;
294  if (x != NULL) {
295  delete x;
296  x = NULL;
297  }
298  }
299  }
300 
302  static Ordinal fromIndirectBytesToCount(const Serializer& vs,
303  const Ordinal bytes,
304  const char charBuffer[],
305  const Ordinal sz = 0) {
306  Ordinal count = 0;
307  Ordinal bytes_used = 0;
308  while (bytes_used < bytes) {
309 
310  // Bytes for degree
311  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
312  bytes_used += b1;
313  charBuffer += b1;
314 
315  // Bytes for taylor coefficients
316  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
317  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
318  bytes_used += b3;
319  charBuffer += b3;
320  bytes_used += *b2;
321  charBuffer += *b2;
322 
323  ++count;
324  }
325  return count;
326  }
327 
329  static void deserialize (const Serializer& vs,
330  const Ordinal bytes,
331  const char charBuffer[],
332  const Ordinal count,
333  TayType buffer[],
334  const Ordinal sz = 0) {
335  for (Ordinal i=0; i<count; i++) {
336 
337  // Deserialize degree
338  Ordinal b1 = iSerT::fromCountToDirectBytes(1);
339  const unsigned int *my_sz = iSerT::convertFromCharPtr(charBuffer);
340  charBuffer += b1;
341 
342  // Create empty Taylor object of given size
343  unsigned int tot_sz = sz;
344  if (sz == 0) tot_sz = *my_sz;
345  buffer[i] = TayType(tot_sz-1, 0.0);
346 
347  // Deserialize taylor coefficients
348  Ordinal b3 = oSerT::fromCountToDirectBytes(1);
349  const Ordinal *b2 = oSerT::convertFromCharPtr(charBuffer);
350  charBuffer += b3;
351  vs.deserialize(*b2, charBuffer, *my_sz,
352  &(buffer[i].fastAccessCoeff(0)));
353  charBuffer += *b2;
354  }
355 
356  }
357 
359 
360  };
361 
363  template <typename Ordinal, typename TayType>
364  struct SerializationTraitsImp {
365 
366  private:
367 
369  typedef typename Sacado::ValueType<TayType>::type ValueT;
370 
372  typedef Teuchos::DefaultSerializer<Ordinal,ValueT> DS;
373 
375  typedef typename DS::DefaultSerializerType ValueSerializer;
376 
378  typedef SerializationImp<Ordinal,TayType,ValueSerializer> Imp;
379 
380  public:
381 
383  static const bool supportsDirectSerialization =
384  Imp::supportsDirectSerialization;
385 
387 
388 
390  static Ordinal fromCountToIndirectBytes(const Ordinal count,
391  const TayType buffer[]) {
392  return Imp::fromCountToIndirectBytes(
393  DS::getDefaultSerializer(), count, buffer);
394  }
395 
397  static void serialize (const Ordinal count,
398  const TayType buffer[],
399  const Ordinal bytes,
400  char charBuffer[]) {
401  Imp::serialize(
402  DS::getDefaultSerializer(), count, buffer, bytes, charBuffer);
403  }
404 
406  static Ordinal fromIndirectBytesToCount(const Ordinal bytes,
407  const char charBuffer[]) {
408  return Imp::fromIndirectBytesToCount(
409  DS::getDefaultSerializer(), bytes, charBuffer);
410  }
411 
413  static void deserialize (const Ordinal bytes,
414  const char charBuffer[],
415  const Ordinal count,
416  TayType buffer[]) {
417  Imp::deserialize(
418  DS::getDefaultSerializer(), bytes, charBuffer, count, buffer);
419  }
420 
422 
423  };
424 
426  template <typename Ordinal, typename TayType, typename ValueSerializer>
427  class SerializerImp {
428 
429  private:
430 
432  typedef SerializationImp<Ordinal,TayType,ValueSerializer> Imp;
433 
435  Teuchos::RCP<const ValueSerializer> vs;
436 
438  Ordinal sz;
439 
440  public:
441 
443  typedef ValueSerializer value_serializer_type;
444 
446  static const bool supportsDirectSerialization =
447  Imp::supportsDirectSerialization;
448 
450  SerializerImp(const Teuchos::RCP<const ValueSerializer>& vs_,
451  Ordinal sz_ = 0) :
452  vs(vs_), sz(sz_) {}
453 
455  Ordinal getSerializerSize() const { return sz; }
456 
458  Teuchos::RCP<const value_serializer_type> getValueSerializer() const {
459  return vs; }
460 
462 
463 
465  Ordinal fromCountToIndirectBytes(const Ordinal count,
466  const TayType buffer[]) const {
467  return Imp::fromCountToIndirectBytes(*vs, count, buffer, sz);
468  }
469 
471  void serialize (const Ordinal count,
472  const TayType buffer[],
473  const Ordinal bytes,
474  char charBuffer[]) const {
475  Imp::serialize(*vs, count, buffer, bytes, charBuffer, sz);
476  }
477 
479  Ordinal fromIndirectBytesToCount(const Ordinal bytes,
480  const char charBuffer[]) const {
481  return Imp::fromIndirectBytesToCount(*vs, bytes, charBuffer, sz);
482  }
483 
485  void deserialize (const Ordinal bytes,
486  const char charBuffer[],
487  const Ordinal count,
488  TayType buffer[]) const {
489  return Imp::deserialize(*vs, bytes, charBuffer, count, buffer, sz);
490  }
491 
493 
494  };
495 
496  } // namespace Tay
497 
498 } // namespace Sacado
499 
500 #endif // HAVE_SACADO_TEUCHOS
501 
502 #endif // SACADO_FAD_SCALARTRAITSIMP_HPP
static std::string eval()
PowExprType< Expr< T1 >, Expr< T2 > >::expr_type pow(const Expr< T1 > &expr1, const Expr< T2 > &expr2)
pow(expr1.val(), expr2.val())
sqrt(expr.val())
int Ordinal
Sacado::Random< double > rnd