Sacado Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
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
53UF4(abs)
54UF4(acos)
55UF4(acosh)
56UF4(asin)
57UF4(asinh)
58UF4(atan)
59UF4(atanh)
60UF4(cos)
61UF4(cosh)
62UF4(exp)
63UF4(fabs)
64UF4(log)
65UF4(log10)
66UF4(sin)
67UF4(sinh)
68UF4(sqrt)
69UF4(tan)
70UF4(tanh)
71
72#undef UF
73#define UF(T,f,fn) T fn(T &x, T &y) { return f(x,y); }
74
75UF4(atan2)
76UF4(pow)
77
78typedef ADVar (*ADVar_uf )(ADVar &);
79typedef DADVar (*DADVar_uf)(DADVar &);
80typedef FDReal (*FDReal_uf)(FDReal &);
81typedef AFReal (*AFReal_uf)(AFReal &);
82
83typedef ADVar (*ADVar_uf2 )(ADVar &, ADVar &);
84typedef DADVar (*DADVar_uf2)(DADVar &, DADVar &);
85typedef FDReal (*FDReal_uf2)(FDReal &, FDReal &);
86typedef AFReal (*AFReal_uf2)(AFReal &, AFReal &);
87
88static 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
95Func4 {
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
106static 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),
113 U4f(atanh, dom_invtrig),
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
128Func42 {
129 const char *name;
130 const double *xdom;
135 };
136
137static 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
147differ(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
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
174main(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);
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);
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);
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);
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);
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 }
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
int main()
Fad specializations for Teuchos::BLAS wrappers.
static void Gradcomp()
static void Hvprod(int n, ADVar **vp, Double *v, Double *hv)
static void Hvprod(int n, ADvar **vp, double *v, double *hv)
static void Gradcomp()
Sacado::Rad::ADvar< FReal > AFReal
#define UF4(f)
static double dom_log[]
ADVar(* ADVar_uf)(ADVar &)
static Func42 UT2[]
DADVar(* DADVar_uf)(DADVar &)
AFReal(* AFReal_uf)(AFReal &)
Sacado::Rad::ADvar< double > DReal
static int differ(double a, double b)
#define In(T)
static double dom_acosh[]
AFReal(* AFReal_uf2)(AFReal &, AFReal &)
Sacado::Rad2::ADvar< double > DADVar
static double dom_all[]
Sacado::Rad2d::ADvar ADVar
Sacado::Fad::SFad< DReal, 1 > FDReal
static Func4 UT[]
ADVar(* ADVar_uf2)(ADVar &, ADVar &)
FDReal(* FDReal_uf)(FDReal &)
static char * progname
DADVar(* DADVar_uf2)(DADVar &, DADVar &)
static int usage(int rc)
#define U4f(f, d)
#define Dcl(T, T1)
FDReal(* FDReal_uf2)(FDReal &, FDReal &)
Sacado::Fad::SFad< double, 1 > FReal
static double dom_invtrig[]
AFReal_uf2 f4
const char * name
DADVar_uf2 f2
FDReal_uf2 f3
ADVar_uf2 f1
const double * xdom
const char * name
ADVar_uf f1
const double * dom
DADVar_uf f2
AFReal_uf f4
FDReal_uf f3
static int rc