Sacado Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hesopcheck.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Sacado Package
5 // Copyright (2007) 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 // Individually check concistency among Hv products via rad2.h, trad2.hpp,
31 // Fad<Rad> and Rad<Fad> for all unary and binary ops.
32 
33 #include <cstdio>
34 #include <float.h> // for DBL_MAX
35 #include "Sacado_DisableKokkosCuda.hpp" // Disable Cuda stuff that fails
36 #include "Sacado_MathFunctions.hpp"
37 #include "Sacado_trad.hpp"
38 #include "Sacado_trad2.hpp"
39 #include "Sacado_rad2.hpp"
40 #include "Sacado_Fad_SFad.hpp"
41 
48 
49 #define UF(T,f,fn) T fn(T &x) { return f(x); }
50 
51 #define UF4(f) UF(ADVar,f,f ## _ADV) UF(DADVar,f,f ## _DADV) UF(FDReal,f,f ## _FDR) UF(AFReal,f,f ## _AFR)
52 
53 UF4(abs)
54 UF4(acos)
55 UF4(acosh)
56 UF4(asin)
57 UF4(asinh)
58 UF4(atan)
59 UF4(atanh)
60 UF4(cos)
61 UF4(cosh)
62 UF4(exp)
63 UF4(fabs)
64 UF4(log)
65 UF4(log10)
66 UF4(sin)
67 UF4(sinh)
68 UF4(sqrt)
69 UF4(tan)
70 UF4(tanh)
71 
72 #undef UF
73 #define UF(T,f,fn) T fn(T &x, T &y) { return f(x,y); }
74 
75 UF4(atan2)
76 UF4(pow)
77 
78 typedef ADVar (*ADVar_uf )(ADVar &);
79 typedef DADVar (*DADVar_uf)(DADVar &);
80 typedef FDReal (*FDReal_uf)(FDReal &);
81 typedef AFReal (*AFReal_uf)(AFReal &);
82 
83 typedef ADVar (*ADVar_uf2 )(ADVar &, ADVar &);
84 typedef DADVar (*DADVar_uf2)(DADVar &, DADVar &);
85 typedef FDReal (*FDReal_uf2)(FDReal &, FDReal &);
86 typedef AFReal (*AFReal_uf2)(AFReal &, AFReal &);
87 
88 static double
89  dom_acosh[] = { 1., DBL_MAX },
90  dom_all[] = { -DBL_MAX, DBL_MAX },
91  dom_invtrig[] = { -1., 1. },
92  dom_log[] = { DBL_MIN, DBL_MAX };
93 
94  struct
95 Func4 {
96  const char *name;
97  const double *dom;
102  };
103 
104 #define U4f(f,d) { #f, d, f ## _ADV, f ## _DADV, f ## _FDR, f ## _AFR }
105 
106 static Func4 UT[] = {
107  U4f(abs, dom_all),
108  U4f(acos, dom_invtrig),
109  U4f(acosh, dom_acosh),
110  U4f(asin, dom_invtrig),
111  U4f(asinh, dom_all),
112  U4f(atan, dom_all),
114  U4f(cos, dom_all),
115  U4f(cosh, dom_all),
116  U4f(exp, dom_all),
117  U4f(fabs, dom_all),
118  U4f(log, dom_log),
119  U4f(log10, dom_log),
120  U4f(sin, dom_all),
121  U4f(sinh, dom_all),
122  U4f(sqrt, dom_log),
123  U4f(tan, dom_all),
124  U4f(tanh, dom_all)
125  };
126 
127  struct
128 Func42 {
129  const char *name;
130  const double *xdom;
135  };
136 
137 static Func42 UT2[] = {
138  U4f(atan2, dom_all),
139  U4f(pow, dom_log)
140  };
141 
142 #undef U4f
143 #undef UF4
144 #undef UF
145 
146  static int
147 differ(double a, double b)
148 {
149  double d = a - b;;
150  if (d == 0.)
151  return 0;
152  if (d < 0.)
153  d = -d;
154  if (a < 0.)
155  a = -a;
156  if (b < 0.)
157  b = - b;
158  return d > 5e-15 * (a + b);
159  }
160 
161  static char *progname;
162 
163  static int
164 usage(int rc)
165 {
166  fprintf(rc ? stderr : stdout, "Usage: %s [-v]\n\
167  to compare consistency among several schemes for Hessian-vector computations.\n\
168  -v ==> verbose (show results)\n\
169  No -v ==> just print OK if all goes well; otherwise complain.\n", progname);
170  return rc;
171  }
172 
173  int
174 main(int argc, char **argv)
175 {
176  Func4 *f, *fe;
177  Func42 *f2, *f2e;
178  char *s;
179  double a, *ap, c, h[4], h1[4], v[3];
180  int b0, b1, b2, k, verbose;
181  static double unopargs[] = { .7, -.9, 2.5, -4.2, .1, .3, 1.2, -2.4, 0. };
182  static double binargs[] = { .1,.2, -.3,.4, -.5,-.6, .7,-.8,
183  1.1,2.2, -2.3,3.4, -2.5,-3.6, 2.7,-3.8, 0.};
184 
185 #define Dcl(T,T1) T x ## T1, y ## T1, z ## T1
186  Dcl(ADVar, ADV);
187  Dcl(DADVar, DADV);
188  Dcl(FDReal, FDR);
189  Dcl(AFReal, AFR);
190 #undef Dcl
191  ADVar *pADV[2];
192  DADVar *pDADV[2];
193  DReal tDR;
194  FDReal tfFDR;
195  FReal tfFR;
196 
197  progname = argv[0];
198  verbose = 0;
199  if (argc > 1) {
200  if (argc > 2 || *(s = argv[1]) != '-')
201  return usage(1);
202  switch(s[1]) {
203  case 'h':
204  case '?':
205  return usage(s[2] != 0);
206  case '-':
207  if (!s[2])
208  break;
209  return usage(strcmp(s,"--help") ? 1 : 0);
210  case 'v':
211  if (!s[2]) {
212  verbose = 1;
213  break;
214  }
215  default:
216  return usage(1);
217  }
218  }
219 
220  v[0] = v[2] = 1.;
221  v[1] = 0.;
222 #define In(T) p ## T[0] = &x ## T; p ## T[1] = &z ## T
223  In(ADV);
224  In(DADV);
225 #undef In
226  b0 = b1 = b2 = 0;
227  fe = UT + sizeof(UT)/sizeof(UT[0]);
228  for(ap = unopargs; (a = *ap++) != 0.; )
229  for(f = UT; f < fe; f++) {
230  if (a < f->dom[0] || a > f->dom[1])
231  continue;
232  xADV = a;
233  yADV = f->f1(xADV);
234  ADVar::Gradcomp();
235  ADVar::Hvprod(1,pADV,v,h);
236  if (verbose) {
237  printf("%s(%g) = %g\n", f->name, xADV.val(), yADV.val());
238  printf("%s' = %g\n", f->name, xADV.adj());
239  printf("%s\" = %g\n", f->name, h[0]);
240  }
241  xDADV = a;
242  yDADV = f->f2(xDADV);
244  if (differ(yDADV.val(), yADV.val())) {
245  ++b0;
246  printf("%s(%g): yADV = %g, yDADV = %g, diff = %.2g\n",
247  f->name, a, yADV.val(), yDADV.val(),
248  yADV.val() - yDADV.val());
249  }
250  if (differ(xADV.adj(), xDADV.adj())) {
251  ++b1;
252  printf("%s'(%g): ADV ==> %g, DADV ==> %g, diff = %.2g\n",
253  f->name, a, xADV.adj(), xDADV.adj(),
254  xADV.adj() - xDADV.adj());
255  }
256  DADVar::Hvprod(1,pDADV,v,h1);
257  if (differ(h[0], h1[0])) {
258  ++b2;
259  printf("%s\"(%g): ADV ==> %g, DADV ==> %g, diff = %.2g\n",
260  f->name, a, h[0], h1[0], h[0] - h1[0]);
261  }
262 
263  tfFDR = 0.;
264  tfFDR.diff(0,1);
265  xFDR = a + tfFDR*v[0];
266  yFDR = f->f3(xFDR);
267  tDR = yFDR.dx(0);
268  DReal::Gradcomp();
269 
270  if (differ(yFDR.val().val(), yADV.val())) {
271  ++b0;
272  printf("%s(%g): yFDR = %g, yADV = %g, diff = %.2g\n",
273  f->name, a, yFDR.val().val(), yADV.val(),
274  yFDR.val().val() - yADV.val());
275  }
276  if (differ(xFDR.val().adj(), h[0])) {
277  ++b2;
278  printf("%s\"(%g): FDR ==> %g, ADV ==> %g, diff = %.2g\n",
279  f->name, a, xFDR.val().adj(), h[0],
280  xFDR.val().adj() - h[0]);
281  }
282 
283  tfFR = 0.;
284  tfFR.diff(0,1);
285  xAFR = a + tfFR*v[0];
286  yAFR = f->f4(xAFR);
288  if (differ(yAFR.val().val(), yADV.val())) {
289  ++b0;
290  printf("%s(%g): yAFR = %g, yADV = %g, diff = %.2g\n",
291  f->name, a, yAFR.val().val(), yADV.val(),
292  yAFR.val().val() - yADV.val());
293  }
294  if (differ(xAFR.adj().val(), xADV.adj())) {
295  ++b1;
296  printf("%s'(%g): AFR ==> %g, ADV ==> %g, diff = %.2g\n",
297  f->name, a, xAFR.adj().val(), xADV.adj(),
298  xAFR.adj().val() - xADV.adj());
299  }
300  if (differ(xAFR.adj().dx(0), h[0])) {
301  ++b2;
302  printf("%s\"(%g): AFR ==> %g, ADV ==> %g, diff = %.2g\n",
303  f->name, a, xAFR.adj().dx(0), h[0],
304  xAFR.adj().dx(0) - h[0]);
305  }
306  }
307  f2e = UT2 + sizeof(UT2)/sizeof(UT2[0]);
308  for(ap = binargs; (a = *ap) != 0.; ap += 2) {
309  c = ap[1];
310  for(f2 = UT2; f2 < f2e; f2++) {
311  if (a < f2->xdom[0] || a > f2->xdom[1])
312  continue;
313  xADV = a;
314  zADV = c;
315  yADV = f2->f1(xADV, zADV);
316  ADVar::Gradcomp();
317  ADVar::Hvprod(2,pADV,v,h);
318  ADVar::Hvprod(2,pADV,v+1,h+2);
319  if (verbose) {
320  printf("%s(%g,%g) = %g\n", f2->name, xADV.val(), zADV.val(), yADV.val());
321  printf("%s' = (%g, %g)\n", f2->name, xADV.adj(), zADV.adj());
322  printf("%s\" = (%8g\t%g)\n(%8g\t%g)", f2->name, h[0],h[1],h[2],h[3]);
323  }
324  xDADV = a;
325  zDADV = c;
326  yDADV = f2->f2(xDADV, zDADV);
328  if (differ(yDADV.val(), yADV.val())) {
329  ++b0;
330  printf("%s(%g,%g): yADV = %g, yDADV = %g, diff = %.2g\n",
331  f2->name, a, c, yADV.val(), yDADV.val(),
332  yADV.val() - yDADV.val());
333  }
334  if (differ(xADV.adj(), xDADV.adj()) || differ(zADV.adj(), zDADV.adj())) {
335  ++b1;
336  printf("%s'(%g,%g): ADV ==> (%g,%g) DADV ==> (%g,%g), diff = (%.2g,%.2g)\n",
337  f2->name, a, c, xADV.adj(), zADV.adj(),
338  xDADV.adj(), zDADV.adj(),
339  xADV.adj() - xDADV.adj(), zADV.adj() - zDADV.adj());
340  }
341  DADVar::Hvprod(2,pDADV,v,h1);
342  DADVar::Hvprod(2,pDADV,v+1,h1+2);
343  for(k = 0; k < 4; k++)
344  if (differ(h[k], h1[k])) {
345  ++b2;
346  printf("%s\"(%g,%g):\n ADV ==> (%8g\t%8g\t%8g\t%g)\n",
347  f2->name, a, c, h[0], h[1], h[2], h[3]);
348  printf("DADV ==> (%8g\t%8g\t%8g\t%g)\n",
349  h1[0], h1[1], h1[2], h1[3]);
350  printf("diff = (%.2g %.2g %.2g %.2g)\n\n",
351  h[0] - h1[0], h[1] - h1[1],
352  h[2] - h1[2], h[3] - h1[3]);
353  break;
354  }
355  tfFDR = 0.;
356  tfFDR.diff(0,1);
357  xFDR = a + tfFDR*v[0];
358  zFDR = c + tfFDR*v[1];
359  yFDR = f2->f3(xFDR, zFDR);
360  tDR = yFDR.dx(0);
361  DReal::Gradcomp();
362 
363  if (differ(yFDR.val().val(), yADV.val())) {
364  ++b0;
365  printf("%s(%g,%g): yFDR = %g, yADV = %g, diff = %.2g\n",
366  f2->name, a, c, yFDR.val().val(), yADV.val(),
367  yFDR.val().val() - yADV.val());
368  }
369  if (differ(xFDR.val().adj(), h[0]) || differ(zFDR.val().adj(), h[1])) {
370  ++b2;
371  printf("%s\"(%g,%g) row 1:\nFDR ==> (%g,%g)\nADV ==> (%g,%g)\n",
372  f2->name, a, c, xFDR.val().adj(), zFDR.val().adj(),
373  h[0], h[1]);
374  printf("diff = %.2g %g\n\n",
375  xFDR.val().adj() - h[0], zFDR.val().adj() - h[1]);
376  }
377  tfFDR.diff(0,1);
378  xFDR = a + tfFDR*v[1];
379  zFDR = c + tfFDR*v[2];
380  yFDR = f2->f3(xFDR, zFDR);
381  tDR = yFDR.dx(0);
382  DReal::Gradcomp();
383  if (differ(xFDR.val().adj(), h[2]) || differ(zFDR.val().adj(), h[3])) {
384  ++b2;
385  printf("%s\"(%g,%g) row 2:\nFDR ==> (%g,%g)\nADV ==> (%g,%g)\n",
386  f2->name, a, c, xFDR.val().adj(), zFDR.val().adj(),
387  h[2], h[3]);
388  printf("diff = %.2g %g\n\n",
389  xFDR.val().adj() - h[2], zFDR.val().adj() - h[3]);
390  }
391 
392  tfFR = 0.;
393  tfFR.diff(0,1);
394  xAFR = a + tfFR*v[0];
395  zAFR = c + tfFR*v[1];
396  yAFR = f2->f4(xAFR, zAFR);
398  if (differ(yAFR.val().val(), yADV.val())) {
399  ++b0;
400  printf("%s(%g,%g): yAFR = %g, yADV = %g, diff = %.2g\n",
401  f2->name, a, c, yAFR.val().val(), yADV.val(),
402  yAFR.val().val() - yADV.val());
403  }
404  if (differ(xAFR.adj().val(), xADV.adj()) || differ(zAFR.adj().val(), zADV.adj())) {
405  ++b1;
406  printf("%s'(%g,%g):\n AFR ==> (%g,%g)\n ADV ==> (%g,%g)\n",
407  f2->name, a, c, xAFR.adj().val(), zAFR.adj().val(),
408  xADV.adj(), zADV.adj());
409  printf(" diff = %.2g %.2g\n\n",
410  xAFR.adj().val() - xADV.adj(),
411  zAFR.adj().val() - zADV.adj());
412  }
413  if (differ(xAFR.adj().dx(0), h[0]) || differ(zAFR.adj().dx(0), h[1])) {
414  ++b2;
415  printf("%s\"(%g,%g) row 1:\n AFR ==> (%g,%g)\n",
416  f2->name, a, c, xAFR.adj().dx(0), zAFR.adj().dx(0));
417  printf(" ADV = (%g,%g)\n diff = %.2g %.2g\n\n", h[0], h[1],
418  xAFR.adj().dx(0) - h[0],
419  zAFR.adj().dx(0) - h[1]);
420  }
421  tfFR.diff(0,1);
422  xAFR = a + tfFR*v[1];
423  zAFR = c + tfFR*v[2];
424  yAFR = f2->f4(xAFR, zAFR);
426  if (differ(xAFR.adj().dx(0), h[2]) || differ(zAFR.adj().dx(0), h[3])) {
427  ++b2;
428  printf("%s\"(%g,%g) row 2:\n AFR ==> (%g,%g)\n",
429  f2->name, a, c, xAFR.adj().dx(0), zAFR.adj().dx(0));
430  printf(" ADV = (%g,%g)\n diff = %.2g %.2g\n\n", h[2], h[3],
431  xAFR.adj().dx(0) - h[2],
432  zAFR.adj().dx(0) - h[3]);
433  }
434  }
435  }
436  k = b0 + b1 + b2;
437  if (k || verbose) {
438  s = progname;
439  if (*s == '.' && s[1] == '/')
440  for(s += 2; *s == '/'; ++s);
441  printf("%s: %d errors seen.\n", s, k);
442  k = 1;
443  }
444  else
445  printf("OK\n");
446  return k;
447  }
asinh(expr.val())
const double * xdom
Definition: hesopcheck.cpp:130
void f()
static void Gradcomp()
asin(expr.val())
cosh(expr.val())
abs(expr.val())
const char * name
Definition: hesopcheck.cpp:129
atanh(expr.val())
Sacado::RadVec::ADvar< double > ADVar
DADVar_uf f2
Definition: hesopcheck.cpp:99
AFReal(* AFReal_uf)(AFReal &)
Definition: hesopcheck.cpp:81
static void Gradcomp()
#define Dcl(T, T1)
Sacado::Fad::SFad< DReal, 1 > FDReal
Definition: hesopcheck.cpp:45
static int usage(int rc)
Definition: hesopcheck.cpp:164
atan(expr.val())
AFReal_uf2 f4
Definition: hesopcheck.cpp:134
KOKKOS_INLINE_FUNCTION mpl::enable_if_c< ExprLevel< Expr< T1 > >::value==ExprLevel< Expr< T2 > >::value, Expr< PowerOp< Expr< T1 >, Expr< T2 > > > >::type pow(const Expr< T1 > &expr1, const Expr< T2 > &expr2)
static int rc
#define U4f(f, d)
Definition: hesopcheck.cpp:104
static double dom_log[]
Definition: hesopcheck.cpp:92
static void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
ADVar(* ADVar_uf)(ADVar &)
Definition: hesopcheck.cpp:78
tanh(expr.val())
FDReal(* FDReal_uf2)(FDReal &, FDReal &)
Definition: hesopcheck.cpp:85
FDReal_uf f3
Definition: hesopcheck.cpp:100
expr expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c *expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr1 c expr2 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 expr2 expr1 expr2 expr1 expr1 expr1 c
#define UF4(f)
Definition: hesopcheck.cpp:51
ADVar(* ADVar_uf2)(ADVar &, ADVar &)
Definition: hesopcheck.cpp:83
Sacado::Fad::SFad< double, 1 > FReal
Definition: hesopcheck.cpp:46
int main()
Definition: ad_example.cpp:191
sqrt(expr.val())
sinh(expr.val())
tan(expr.val())
Sacado::Rad2::ADvar< double > DADVar
Definition: hesopcheck.cpp:43
static int differ(double a, double b)
Definition: hesopcheck.cpp:147
DADVar_uf2 f2
Definition: hesopcheck.cpp:132
atan2(expr1.val(), expr2.val())
static double dom_all[]
Definition: hesopcheck.cpp:90
Sacado::Rad::ADvar< FReal > AFReal
Definition: hesopcheck.cpp:47
static Func4 UT[]
Definition: hesopcheck.cpp:106
ADVar_uf f1
Definition: hesopcheck.cpp:98
sin(expr.val())
static double dom_acosh[]
Definition: hesopcheck.cpp:89
AFReal(* AFReal_uf2)(AFReal &, AFReal &)
Definition: hesopcheck.cpp:86
log(expr.val())
ADVar_uf2 f1
Definition: hesopcheck.cpp:131
const char * name
Definition: hesopcheck.cpp:96
acosh(expr.val())
const double * dom
Definition: hesopcheck.cpp:97
static char * progname
Definition: hesopcheck.cpp:161
acos(expr.val())
#define In(T)
AFReal_uf f4
Definition: hesopcheck.cpp:101
static Func42 UT2[]
Definition: hesopcheck.cpp:137
DADVar(* DADVar_uf2)(DADVar &, DADVar &)
Definition: hesopcheck.cpp:84
DADVar(* DADVar_uf)(DADVar &)
Definition: hesopcheck.cpp:79
exp(expr.val())
FDReal_uf2 f3
Definition: hesopcheck.cpp:133
fabs(expr.val())
log10(expr.val())
cos(expr.val())
FDReal(* FDReal_uf)(FDReal &)
Definition: hesopcheck.cpp:80
Sacado::Rad::ADvar< double > DReal
Definition: hesopcheck.cpp:44
static double dom_invtrig[]
Definition: hesopcheck.cpp:91