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