Sacado Package Browser (Single Doxygen Collection)  Version of the Day
NestedFadUnitTests.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Sacado Package
5 // Copyright (2006) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
23 // USA
24 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
25 // (etphipp@sandia.gov).
26 //
27 // ***********************************************************************
28 // @HEADER
29 
30 #ifndef NESTED_FADUNITTESTS_HPP
31 #define NESTED_FADUNITTESTS_HPP
32 
33 // Sacado includes
34 #include "Sacado_No_Kokkos.hpp"
35 #include "Sacado_Random.hpp"
36 
37 // Fad includes
38 #include "Fad/fad.h"
39 
40 // Cppunit includes
41 #include <cppunit/extensions/HelperMacros.h>
42 
43 #define COMPARE_VALUES(a, b) \
44  CPPUNIT_ASSERT( std::abs(a-b) < this->tol_a + this->tol_r*std::abs(a) );
45 
46 #define COMPARE_FADS(a, b) \
47 CPPUNIT_ASSERT(a.size() == b.size()); \
48 CPPUNIT_ASSERT(a.hasFastAccess() == b.hasFastAccess()); \
49 COMPARE_VALUES(a.val(), b.val()); \
50 for (int zz=0; zz<a.size(); zz++) { \
51  COMPARE_VALUES(a.dx(zz), b.dx(zz)); \
52  COMPARE_VALUES(a.fastAccessDx(zz), b.fastAccessDx(zz)); \
53  } \
54  ;
55 
56 #define COMPARE_NESTED_FADS(a, b) \
57 CPPUNIT_ASSERT(a.size() == b.size()); \
58 CPPUNIT_ASSERT(a.hasFastAccess() == b.hasFastAccess()); \
59 COMPARE_FADS(a.val(), b.val()); \
60 for (int z=0; z<a.size(); z++) { \
61  COMPARE_FADS(a.dx(z), b.dx(z)); \
62  COMPARE_FADS(a.fastAccessDx(z), b.fastAccessDx(z)); \
63  } \
64  ;
65 
66 #define BINARY_OP_TEST(TESTNAME,OP) \
67  void TESTNAME () { \
68  c_dfad = a_dfad OP b_dfad; \
69  c_fad = a_fad OP b_fad; \
70  COMPARE_NESTED_FADS(c_dfad, c_fad); \
71  \
72  double val = urand.number(); \
73  c_dfad = a_dfad OP val; \
74  c_fad = a_fad OP FadType(val); \
75  COMPARE_NESTED_FADS(c_dfad, c_fad); \
76  \
77  c_dfad = val OP b_dfad; \
78  c_fad = FadType(val) OP b_fad; \
79  COMPARE_NESTED_FADS(c_dfad, c_fad); \
80  }
81 
82 #define RELOP_TEST(TESTNAME,OP) \
83  void TESTNAME () { \
84  bool r1 = a_dfad OP b_dfad; \
85  bool r2 = a_fad OP b_fad; \
86  CPPUNIT_ASSERT(r1 == r2); \
87  \
88  double val = urand.number(); \
89  r1 = a_dfad OP val; \
90  r2 = a_fad OP FadType(val); \
91  CPPUNIT_ASSERT(r1 == r2); \
92  \
93  r1 = val OP b_dfad; \
94  r2 = FadType(val) OP b_fad; \
95  CPPUNIT_ASSERT(r1 == r2); \
96  }
97 
98 #define BINARY_FUNC_TEST(TESTNAME,FUNC) \
99  void TESTNAME () { \
100  c_dfad = FUNC (a_dfad,b_dfad); \
101  c_fad = FUNC (a_fad,b_fad); \
102  COMPARE_NESTED_FADS(c_dfad, c_fad); \
103  \
104  double val = urand.number(); \
105  c_dfad = FUNC (a_dfad,val); \
106  c_fad = FUNC (a_fad,FadType(val)); \
107  COMPARE_NESTED_FADS(c_dfad, c_fad); \
108  \
109  c_dfad = FUNC (val,b_dfad); \
110  c_fad = FUNC (FadType(val),b_fad); \
111  COMPARE_NESTED_FADS(c_dfad, c_fad); \
112  }
113 
114 #define UNARY_OP_TEST(TESTNAME,OP) \
115  void TESTNAME () { \
116  c_dfad = OP a_dfad; \
117  c_fad = OP a_fad; \
118  COMPARE_NESTED_FADS(c_dfad, c_fad); \
119  }
120 
121 #define UNARY_FUNC_TEST(TESTNAME,FUNC) \
122  void TESTNAME () { \
123  c_dfad = FUNC (a_dfad); \
124  c_fad = FUNC (a_fad); \
125  COMPARE_NESTED_FADS(c_dfad, c_fad); \
126  }
127 
128 #define UNARY_ASSIGNOP_TEST(TESTNAME,OP) \
129  void TESTNAME () { \
130  c_dfad OP a_dfad; \
131  c_fad OP a_fad; \
132  COMPARE_NESTED_FADS(c_dfad, c_fad); \
133  \
134  double val = urand.number(); \
135  c_dfad OP val; \
136  c_fad OP FadType(val); \
137  COMPARE_NESTED_FADS(c_dfad, c_fad); \
138  }
139 
140 // A class for testing each Fad operation
141 template <class FadFadType, class ScalarType>
142 class FadFadOpsUnitTest : public CppUnit::TestFixture {
143 
145 
146  CPPUNIT_TEST(testAddition);
147  CPPUNIT_TEST(testSubtraction);
148  CPPUNIT_TEST(testMultiplication);
149  CPPUNIT_TEST(testDivision);
150 
151  CPPUNIT_TEST(testEquals);
152  CPPUNIT_TEST(testNotEquals);
153  CPPUNIT_TEST(testLessThanOrEquals);
154  CPPUNIT_TEST(testGreaterThanOrEquals);
155  CPPUNIT_TEST(testLessThan);
156  CPPUNIT_TEST(testGreaterThan);
157 
158  CPPUNIT_TEST(testPow);
161 
162  CPPUNIT_TEST(testUnaryPlus);
163  CPPUNIT_TEST(testUnaryMinus);
164 
165  CPPUNIT_TEST(testExp);
166  CPPUNIT_TEST(testLog);
167  CPPUNIT_TEST(testLog10);
168  CPPUNIT_TEST(testSqrt);
169  CPPUNIT_TEST(testCos);
170  CPPUNIT_TEST(testSin);
171  CPPUNIT_TEST(testTan);
172  CPPUNIT_TEST(testACos);
173  CPPUNIT_TEST(testASin);
174  CPPUNIT_TEST(testATan);
175  CPPUNIT_TEST(testCosh);
176  CPPUNIT_TEST(testSinh);
177  CPPUNIT_TEST(testTanh);
178  CPPUNIT_TEST(testAbs);
179  CPPUNIT_TEST(testFAbs);
180 
181  CPPUNIT_TEST(testPlusEquals);
182  CPPUNIT_TEST(testMinusEquals);
183  CPPUNIT_TEST(testTimesEquals);
184  CPPUNIT_TEST(testDivideEquals);
185 
187 
192 
194 
195 public:
196 
197  typedef typename FadFadType::value_type FadType;
198 
200 
201  FadFadOpsUnitTest(int numComponents1, int numComponents2,
202  ScalarType absolute_tolerance,
203  ScalarType relative_tolerance);
204 
205  void setUp();
206 
207  void tearDown();
208 
209  BINARY_OP_TEST(testAddition, +);
210  BINARY_OP_TEST(testSubtraction, -);
211  BINARY_OP_TEST(testMultiplication, *);
212  BINARY_OP_TEST(testDivision, /);
213 
214  RELOP_TEST(testEquals, ==);
215  RELOP_TEST(testNotEquals, !=);
216  RELOP_TEST(testLessThanOrEquals, <=);
217  RELOP_TEST(testGreaterThanOrEquals, >=);
218  RELOP_TEST(testLessThan, <);
219  RELOP_TEST(testGreaterThan, >);
220 
221  BINARY_FUNC_TEST(testPow, pow);
222 
223  UNARY_OP_TEST(testUnaryPlus, +);
224  UNARY_OP_TEST(testUnaryMinus, -);
225 
226  UNARY_FUNC_TEST(testExp, exp);
227  UNARY_FUNC_TEST(testLog, log);
228  UNARY_FUNC_TEST(testLog10, log10);
229  UNARY_FUNC_TEST(testSqrt, sqrt);
230  UNARY_FUNC_TEST(testCos, cos);
231  UNARY_FUNC_TEST(testSin, sin);
232  UNARY_FUNC_TEST(testTan, tan);
233  UNARY_FUNC_TEST(testACos, acos);
234  UNARY_FUNC_TEST(testASin, asin);
235  UNARY_FUNC_TEST(testATan, atan);
236  UNARY_FUNC_TEST(testCosh, cosh);
237  UNARY_FUNC_TEST(testSinh, sinh);
238  UNARY_FUNC_TEST(testTanh, tanh);
239  UNARY_FUNC_TEST(testAbs, abs);
240  UNARY_FUNC_TEST(testFAbs, fabs);
241 
242  UNARY_ASSIGNOP_TEST(testPlusEquals, +=);
243  UNARY_ASSIGNOP_TEST(testMinusEquals, -=);
244  UNARY_ASSIGNOP_TEST(testTimesEquals, *=);
245  UNARY_ASSIGNOP_TEST(testDivideEquals, /=);
246 
247  void testMax();
248  void testMin();
249 
250  template <typename ScalarT>
251  ScalarT composite1(const ScalarT& a, const ScalarT& b) {
252  ScalarT t1 = 3. * a + sin(b) / log(fabs(a - b * 7.));
253  ScalarT t2 = 1.0e3;
254  ScalarT t3 = 5.7e4;
255  ScalarT t4 = 3.2e5;
256  t1 *= cos(a + exp(t1)) / 6. - tan(t1*sqrt(abs(a * log10(abs(b)))));
257  t1 -= acos((6.+asin(pow(fabs(a),b)/t2))/t3) * asin(pow(fabs(b),2.)*1.0/t4) * atan((b*pow(2.,log(abs(a))))/(t3*t4));
258  t1 /= cosh(b - 0.7) + 7.*sinh(t1 + 0.8)*tanh(9./a) - 9.;
259  t1 += pow(abs(a*4.),b-8.)/cos(a*b*a);
260 
261  return t1;
262  }
263 
264  template <typename ScalarT>
265  ScalarT composite1_fad(const ScalarT& a, const ScalarT& b) {
266  ScalarT t1 = FadType(3.) * a + sin(b) / log(fabs(a - b * FadType(7.)));
267  ScalarT t2 = FadType(1.0e3);
268  ScalarT t3 = FadType(7e4);
269  ScalarT t4 = FadType(3.2e5);
270  t1 *= cos(a + exp(t1)) / FadType(6.) - tan(t1*sqrt(abs(a * log10(abs(b)))));
271  t1 -= acos((FadType(6.)+asin(pow(fabs(a),b)/t2))/t3) * asin(pow(fabs(b),FadType(2.))*FadType(1.0)/t4) * atan((b*pow(FadType(2.),log(abs(a))))/(t3*t4));
272  t1 /= cosh(b - FadType(0.7)) + FadType(7.)*sinh(t1 + FadType(0.8))*tanh(FadType(9.)/a) - FadType(9.);
273  t1 += pow(abs(a*FadType(4.)),b-FadType(8.))/cos(a*b*a);
274 
275  return t1;
276  }
277 
278  void testComposite1() {
282  }
283 
284  void testPlusLR() {
285  FadFadType aa_dfad = a_dfad;
286  FAD::Fad< FadType > aa_fad = a_fad;
287  aa_dfad = 1.0;
288  aa_fad = FadType(1.0);
289  aa_dfad = aa_dfad + b_dfad;
290  aa_fad = aa_fad + b_fad;
291  COMPARE_NESTED_FADS(aa_dfad, aa_fad);
292  }
293 
294  void testMinusLR() {
295  FadFadType aa_dfad = a_dfad;
296  FAD::Fad< FadType > aa_fad = a_fad;
297  aa_dfad = 1.0;
298  aa_fad = FadType(1.0);
299  aa_dfad = aa_dfad - b_dfad;
300  aa_fad = aa_fad - b_fad;
301  COMPARE_NESTED_FADS(aa_dfad, aa_fad);
302  }
303 
304  void testTimesLR() {
305  FadFadType aa_dfad = a_dfad;
306  FAD::Fad< FadType > aa_fad = a_fad;
307  aa_dfad = 2.0;
308  aa_fad = FadType(2.0);
309  aa_dfad = aa_dfad * b_dfad;
310  aa_fad = aa_fad * b_fad;
311  COMPARE_NESTED_FADS(aa_dfad, aa_fad);
312  }
313 
314  void testDivideLR() {
315  FadFadType aa_dfad = a_dfad;
316  FAD::Fad< FadType > aa_fad = a_fad;
317  aa_dfad = 2.0;
318  aa_fad = FadType(2.0);
319  aa_dfad = aa_dfad / b_dfad;
320  aa_fad = aa_fad / b_fad;
321  COMPARE_NESTED_FADS(aa_dfad, aa_fad);
322  }
323 
324 protected:
325 
326  // DFad variables
327  FadFadType a_dfad, b_dfad, c_dfad;
328 
329  // Fad variables
330  FAD::Fad<FadType> a_fad, b_fad, c_fad;
331 
332  // Random number generator
334 
335  // Number of derivative components
336  int n1, n2;
337 
338  // Tolerances to which fad objects should be the same
339  ScalarType tol_a, tol_r;
340 
341 }; // class FadFadOpsUnitTest
342 
343 // Set the random number generator with a specific seed to prevent NaN's
344 // in the second derivatives, likely due to overflow
345 
346 template <class FadFadType, class ScalarType>
349  urand(0.0, 1.0, 123456), n1(5), n2(3), tol_a(1.0e-15), tol_r(1.0e-14) {}
350 
351 template <class FadFadType, class ScalarType>
353 FadFadOpsUnitTest(int numComponents1, int numComponents2,
354  ScalarType absolute_tolerance,
355  ScalarType relative_tolerance) :
356  urand(0.0, 1.0, 123456),
357  n1(numComponents1),
358  n2(numComponents2),
359  tol_a(absolute_tolerance),
360  tol_r(relative_tolerance) {}
361 
362 template <class FadFadType, class ScalarType>
364  ScalarType val;
365 
366  val = urand.number();
367  a_dfad = FadFadType(n1,FadType(n2,val));
368  a_fad = FAD::Fad<FadType>(n1,FadType(n2,val));
369 
370  val = urand.number();
371  b_dfad = FadFadType(n1,FadType(n2,val));
372  b_fad = FAD::Fad<FadType>(n1,FadType(n2,val));
373 
374  for (int j=0; j<n2; j++) {
375  ScalarType val2;
376  val2 = urand.number();
377  a_dfad.val().fastAccessDx(j) = val2;
378  a_fad.val().fastAccessDx(j) = val2;
379 
380  val2 = urand.number();
381  b_dfad.val().fastAccessDx(j) = val2;
382  b_fad.val().fastAccessDx(j) = val2;
383  }
384 
385  for (int i=0; i<n1; i++) {
386  val = urand.number();
387  a_dfad.fastAccessDx(i) = FadType(n2,val);
388  a_fad.fastAccessDx(i) = FadType(n2,val);
389 
390  val = urand.number();
391  b_dfad.fastAccessDx(i) = FadType(n2,val);
392  b_fad.fastAccessDx(i) = FadType(n2,val);
393 
394  for (int j=0; j<n2; j++) {
395  ScalarType val2;
396  val2 = urand.number();
397  a_dfad.fastAccessDx(i).fastAccessDx(j) = val2;
398  a_fad.fastAccessDx(i).fastAccessDx(j) = val2;
399 
400  val2 = urand.number();
401  b_dfad.fastAccessDx(i).fastAccessDx(j) = val2;
402  b_fad.fastAccessDx(i).fastAccessDx(j) = val2;
403  }
404  }
405 }
406 
407 template <class FadFadType, class ScalarType>
410 
411 template <class FadFadType, class ScalarType>
414  FadType val;
415 
416  FadFadType aa_dfad = a_dfad + 1.0;
417  c_dfad = max(aa_dfad, a_dfad);
418  COMPARE_FADS(c_dfad.val(), aa_dfad.val());
419  for (int i=0; i<n1; i++) {
420  COMPARE_FADS(c_dfad.dx(i), aa_dfad.dx(i));
421  COMPARE_FADS(c_dfad.fastAccessDx(i), aa_dfad.fastAccessDx(i));
422  }
423 
424  c_dfad = max(a_dfad, aa_dfad);
425  COMPARE_FADS(c_dfad.val(), aa_dfad.val());
426  for (int i=0; i<n1; i++) {
427  COMPARE_FADS(c_dfad.dx(i), aa_dfad.dx(i));
428  COMPARE_FADS(c_dfad.fastAccessDx(i), aa_dfad.fastAccessDx(i));
429  }
430 
431  c_dfad = max(a_dfad+1.0, a_dfad);
432  COMPARE_FADS(c_dfad.val(), aa_dfad.val());
433  for (int i=0; i<n1; i++) {
434  COMPARE_FADS(c_dfad.dx(i), aa_dfad.dx(i));
435  COMPARE_FADS(c_dfad.fastAccessDx(i), aa_dfad.fastAccessDx(i));
436  }
437 
438  c_dfad = max(a_dfad, a_dfad+1.0);
439  COMPARE_FADS(c_dfad.val(), aa_dfad.val());
440  for (int i=0; i<n1; i++) {
441  COMPARE_FADS(c_dfad.dx(i), aa_dfad.dx(i));
442  COMPARE_FADS(c_dfad.fastAccessDx(i), aa_dfad.fastAccessDx(i));
443  }
444 
445  val = a_dfad.val() + 1;
446  c_dfad = max(a_dfad, val);
447  COMPARE_FADS(c_dfad.val(), val);
448  for (int i=0; i<n1; i++) {
449  COMPARE_FADS(c_dfad.dx(i), FadType(0.0));
450  }
451 
452  val = a_dfad.val() - 1;
453  c_dfad = max(a_dfad, val);
454  COMPARE_FADS(c_dfad.val(), a_dfad.val());
455  for (int i=0; i<n1; i++) {
456  COMPARE_FADS(c_dfad.dx(i), a_dfad.dx(i));
457  COMPARE_FADS(c_dfad.fastAccessDx(i), a_dfad.fastAccessDx(i));
458  }
459 
460  val = b_dfad.val() + 1;
461  c_dfad = max(val, b_dfad);
462  COMPARE_FADS(c_dfad.val(), val);
463  for (int i=0; i<n1; i++) {
464  COMPARE_FADS(c_dfad.dx(i), FadType(0.0));
465  }
466 
467  val = b_dfad.val() - 1;
468  c_dfad = max(val, b_dfad);
469  COMPARE_FADS(c_dfad.val(), b_dfad.val());
470  for (int i=0; i<n1; i++) {
471  COMPARE_FADS(c_dfad.dx(i), b_dfad.dx(i));
472  COMPARE_FADS(c_dfad.fastAccessDx(i), b_dfad.fastAccessDx(i));
473  }
474 }
475 
476 template <class FadFadType, class ScalarType>
479  FadType val;
480 
481  FadFadType aa_dfad = a_dfad - 1.0;
482  c_dfad = min(aa_dfad, a_dfad);
483  COMPARE_FADS(c_dfad.val(), aa_dfad.val());
484  for (int i=0; i<n1; i++) {
485  COMPARE_FADS(c_dfad.dx(i), aa_dfad.dx(i));
486  COMPARE_FADS(c_dfad.fastAccessDx(i), aa_dfad.fastAccessDx(i));
487  }
488 
489  c_dfad = min(a_dfad, aa_dfad);
490  COMPARE_FADS(c_dfad.val(), aa_dfad.val());
491  for (int i=0; i<n1; i++) {
492  COMPARE_FADS(c_dfad.dx(i), aa_dfad.dx(i));
493  COMPARE_FADS(c_dfad.fastAccessDx(i), aa_dfad.fastAccessDx(i));
494  }
495 
496  val = a_dfad.val() - 1;
497  c_dfad = min(a_dfad, val);
498  COMPARE_FADS(c_dfad.val(), val);
499  for (int i=0; i<n1; i++) {
500  COMPARE_FADS(c_dfad.dx(i), FadType(0.0));
501  }
502 
503  val = a_dfad.val() + 1;
504  c_dfad = min(a_dfad, val);
505  COMPARE_FADS(c_dfad.val(), a_dfad.val());
506  for (int i=0; i<n1; i++) {
507  COMPARE_FADS(c_dfad.dx(i), a_dfad.dx(i));
508  COMPARE_FADS(c_dfad.fastAccessDx(i), a_dfad.fastAccessDx(i));
509  }
510 
511  val = b_dfad.val() - 1;
512  c_dfad = min(val, b_dfad);
513  COMPARE_FADS(c_dfad.val(), val);
514  for (int i=0; i<n1; i++) {
515  COMPARE_FADS(c_dfad.dx(i), FadType(0.0));
516  }
517 
518  val = b_dfad.val() + 1;
519  c_dfad = min(val, b_dfad);
520  COMPARE_FADS(c_dfad.val(), b_dfad.val());
521  for (int i=0; i<n1; i++) {
522  COMPARE_FADS(c_dfad.dx(i), b_dfad.dx(i));
523  COMPARE_FADS(c_dfad.fastAccessDx(i), b_dfad.fastAccessDx(i));
524  }
525 }
526 
527 #undef COMPARE_VALUES
528 #undef COMPARE_FADS
529 #undef COMPARE_NESTED_FADS
530 
531 #endif // NESETD_FADUNITTESTS_HPP
ScalarT composite1(const ScalarT &a, const ScalarT &b)
asin(expr.val())
cosh(expr.val())
abs(expr.val())
#define COMPARE_FADS(a, b)
UNARY_OP_TEST(testUnaryPlus,+)
FAD::Fad< FadType > a_fad
Sacado::Fad::DFad< double > FadType
pow(expr1.val(), expr2.val())
BINARY_FUNC_TEST(testPow, pow)
atan(expr.val())
FAD::Fad< FadType > b_fad
BINARY_OP_TEST(testAddition,+)
FAD::Fad< FadType > c_fad
expr val()
tanh(expr.val())
#define COMPARE_NESTED_FADS(a, b)
SimpleFad< ValueT > min(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
CPPUNIT_TEST(testAddition)
UNARY_FUNC_TEST(testExp, exp)
sqrt(expr.val())
sinh(expr.val())
tan(expr.val())
Sacado::Random< ScalarType > urand
UNARY_ASSIGNOP_TEST(testPlusEquals,+=)
FadFadType::value_type FadType
sin(expr.val())
log(expr.val())
acos(expr.val())
SimpleFad< ValueT > max(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
ScalarT composite1_fad(const ScalarT &a, const ScalarT &b)
exp(expr.val())
RELOP_TEST(testEquals,==)
fabs(expr.val())
CPPUNIT_TEST_SUITE(FadFadOpsUnitTest)
log10(expr.val())
cos(expr.val())