Sacado Package Browser (Single Doxygen Collection)  Version of the Day
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  }
FDReal(* FDReal_uf)(FDReal &)
Definition: hesopcheck.cpp:80
const double * xdom
Definition: hesopcheck.cpp:130
void f()
static void Gradcomp()
ADVar(* ADVar_uf2)(ADVar &, ADVar &)
Definition: hesopcheck.cpp:83
const char * name
Definition: hesopcheck.cpp:129
DADVar_uf f2
Definition: hesopcheck.cpp:99
SimpleFad< ValueT > atan2(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
Sacado::Rad2d::ADvar ADVar
Definition: hesopcheck.cpp:42
AFReal(* AFReal_uf)(AFReal &)
Definition: hesopcheck.cpp:81
static void Gradcomp()
SimpleFad< ValueT > sqrt(const SimpleFad< ValueT > &a)
#define Dcl(T, T1)
SimpleFad< ValueT > asin(const SimpleFad< ValueT > &a)
Sacado::Fad::SFad< DReal, 1 > FDReal
Definition: hesopcheck.cpp:45
static int usage(int rc)
Definition: hesopcheck.cpp:164
AFReal_uf2 f4
Definition: hesopcheck.cpp:134
FDReal(* FDReal_uf2)(FDReal &, FDReal &)
Definition: hesopcheck.cpp:85
SimpleFad< ValueT > log(const SimpleFad< ValueT > &a)
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)
SimpleFad< ValueT > exp(const SimpleFad< ValueT > &a)
SimpleFad< ValueT > atan(const SimpleFad< ValueT > &a)
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
AFReal(* AFReal_uf2)(AFReal &, AFReal &)
Definition: hesopcheck.cpp:86
SimpleFad< ValueT > log10(const SimpleFad< ValueT > &a)
Sacado::Fad::SFad< double, 1 > FReal
Definition: hesopcheck.cpp:46
DADVar(* DADVar_uf)(DADVar &)
Definition: hesopcheck.cpp:79
SimpleFad< ValueT > sin(const SimpleFad< ValueT > &a)
Sacado::Rad2::ADvar< double > DADVar
Definition: hesopcheck.cpp:43
static int differ(double a, double b)
Definition: hesopcheck.cpp:147
SimpleFad< ValueT > sinh(const SimpleFad< ValueT > &a)
DADVar_uf2 f2
Definition: hesopcheck.cpp:132
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
SimpleFad< ValueT > pow(const SimpleFad< ValueT > &a, const SimpleFad< ValueT > &b)
static double dom_acosh[]
Definition: hesopcheck.cpp:89
SimpleFad< ValueT > tan(const SimpleFad< ValueT > &a)
expr acosh(expr.val())) FAD_UNARYOP_MACRO(asinh
expr asinh(expr.val())) FAD_UNARYOP_MACRO(atanh
ADVar_uf2 f1
Definition: hesopcheck.cpp:131
SimpleFad< ValueT > cos(const SimpleFad< ValueT > &a)
const char * name
Definition: hesopcheck.cpp:96
SimpleFad< ValueT > cosh(const SimpleFad< ValueT > &a)
const double * dom
Definition: hesopcheck.cpp:97
static char * progname
Definition: hesopcheck.cpp:161
#define In(T)
AFReal_uf f4
Definition: hesopcheck.cpp:101
static Func42 UT2[]
Definition: hesopcheck.cpp:137
KOKKOS_INLINE_FUNCTION Expr< AbsOp< Expr< T > > > abs(const Expr< T > &expr)
KOKKOS_INLINE_FUNCTION Expr< FAbsOp< Expr< T > > > fabs(const Expr< T > &expr)
FDReal_uf2 f3
Definition: hesopcheck.cpp:133
ADVar(* ADVar_uf)(ADVar &)
Definition: hesopcheck.cpp:78
expr atanh(expr.val())) FAD_UNARYOP_MACRO(abs
SimpleFad< ValueT > acos(const SimpleFad< ValueT > &a)
SimpleFad< ValueT > tanh(const SimpleFad< ValueT > &a)
int main(int argc, char **argv)
Definition: hesopcheck.cpp:174
DADVar(* DADVar_uf2)(DADVar &, DADVar &)
Definition: hesopcheck.cpp:84
Sacado::Rad::ADvar< double > DReal
Definition: hesopcheck.cpp:44
static double dom_invtrig[]
Definition: hesopcheck.cpp:91