|
Sacado Package Browser (Single Doxygen Collection) Version of the Day
|
|
Sacado Package Browser (Single Doxygen Collection) Version of the Day
|
00001 // @HEADER 00002 // *********************************************************************** 00003 // 00004 // Sacado Package 00005 // Copyright (2007) Sandia Corporation 00006 // 00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation, 00008 // the U.S. Government retains certain rights in this software. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps 00025 // (etphipp@sandia.gov). 00026 // 00027 // *********************************************************************** 00028 // @HEADER 00029 00030 // Individually check concistency among Hv products via rad2.h, trad2.hpp, 00031 // Fad<Rad> and Rad<Fad> for all unary and binary ops. 00032 00033 #include <cstdio> 00034 #include <float.h> // for DBL_MAX 00035 #include "Sacado_MathFunctions.hpp" 00036 #include "Sacado_trad.hpp" 00037 #include "Sacado_trad2.hpp" 00038 #include "Sacado_rad2.hpp" 00039 #include "Sacado_Fad_SFad.hpp" 00040 00041 typedef Sacado::Rad2d::ADvar ADVar; 00042 typedef Sacado::Rad2::ADvar<double> DADVar; 00043 typedef Sacado::Rad::ADvar<double> DReal; 00044 typedef Sacado::Fad::SFad<DReal,1> FDReal; 00045 typedef Sacado::Fad::SFad<double,1> FReal; 00046 typedef Sacado::Rad::ADvar<FReal> AFReal; 00047 00048 #define UF(T,f,fn) T fn(T &x) { return f(x); } 00049 00050 #define UF4(f) UF(ADVar,f,f ## _ADV) UF(DADVar,f,f ## _DADV) UF(FDReal,f,f ## _FDR) UF(AFReal,f,f ## _AFR) 00051 00052 UF4(abs) 00053 UF4(acos) 00054 UF4(acosh) 00055 UF4(asin) 00056 UF4(asinh) 00057 UF4(atan) 00058 UF4(atanh) 00059 UF4(cos) 00060 UF4(cosh) 00061 UF4(exp) 00062 UF4(fabs) 00063 UF4(log) 00064 UF4(log10) 00065 UF4(sin) 00066 UF4(sinh) 00067 UF4(sqrt) 00068 UF4(tan) 00069 UF4(tanh) 00070 00071 #undef UF 00072 #define UF(T,f,fn) T fn(T &x, T &y) { return f(x,y); } 00073 00074 UF4(atan2) 00075 UF4(pow) 00076 00077 typedef ADVar (*ADVar_uf )(ADVar &); 00078 typedef DADVar (*DADVar_uf)(DADVar &); 00079 typedef FDReal (*FDReal_uf)(FDReal &); 00080 typedef AFReal (*AFReal_uf)(AFReal &); 00081 00082 typedef ADVar (*ADVar_uf2 )(ADVar &, ADVar &); 00083 typedef DADVar (*DADVar_uf2)(DADVar &, DADVar &); 00084 typedef FDReal (*FDReal_uf2)(FDReal &, FDReal &); 00085 typedef AFReal (*AFReal_uf2)(AFReal &, AFReal &); 00086 00087 static double 00088 dom_acosh[] = { 1., DBL_MAX }, 00089 dom_all[] = { -DBL_MAX, DBL_MAX }, 00090 dom_invtrig[] = { -1., 1. }, 00091 dom_log[] = { DBL_MIN, DBL_MAX }; 00092 00093 struct 00094 Func4 { 00095 const char *name; 00096 const double *dom; 00097 ADVar_uf f1; 00098 DADVar_uf f2; 00099 FDReal_uf f3; 00100 AFReal_uf f4; 00101 }; 00102 00103 #define U4f(f,d) { #f, d, f ## _ADV, f ## _DADV, f ## _FDR, f ## _AFR } 00104 00105 static Func4 UT[] = { 00106 U4f(abs, dom_all), 00107 U4f(acos, dom_invtrig), 00108 U4f(acosh, dom_acosh), 00109 U4f(asin, dom_invtrig), 00110 U4f(asinh, dom_all), 00111 U4f(atan, dom_all), 00112 U4f(atanh, dom_invtrig), 00113 U4f(cos, dom_all), 00114 U4f(cosh, dom_all), 00115 U4f(exp, dom_all), 00116 U4f(fabs, dom_all), 00117 U4f(log, dom_log), 00118 U4f(log10, dom_log), 00119 U4f(sin, dom_all), 00120 U4f(sinh, dom_all), 00121 U4f(sqrt, dom_log), 00122 U4f(tan, dom_all), 00123 U4f(tanh, dom_all) 00124 }; 00125 00126 struct 00127 Func42 { 00128 const char *name; 00129 const double *xdom; 00130 ADVar_uf2 f1; 00131 DADVar_uf2 f2; 00132 FDReal_uf2 f3; 00133 AFReal_uf2 f4; 00134 }; 00135 00136 static Func42 UT2[] = { 00137 U4f(atan2, dom_all), 00138 U4f(pow, dom_log) 00139 }; 00140 00141 #undef U4f 00142 #undef UF4 00143 #undef UF 00144 00145 static int 00146 differ(double a, double b) 00147 { 00148 double d = a - b;; 00149 if (d == 0.) 00150 return 0; 00151 if (d < 0.) 00152 d = -d; 00153 if (a < 0.) 00154 a = -a; 00155 if (b < 0.) 00156 b = - b; 00157 return d > 5e-15 * (a + b); 00158 } 00159 00160 static char *progname; 00161 00162 static int 00163 usage(int rc) 00164 { 00165 fprintf(rc ? stderr : stdout, "Usage: %s [-v]\n\ 00166 to compare consistency among several schemes for Hessian-vector computations.\n\ 00167 -v ==> verbose (show results)\n\ 00168 No -v ==> just print OK if all goes well; otherwise complain.\n", progname); 00169 return rc; 00170 } 00171 00172 int 00173 main(int argc, char **argv) 00174 { 00175 Func4 *f, *fe; 00176 Func42 *f2, *f2e; 00177 char *s; 00178 double a, *ap, c, h[4], h1[4], v[3]; 00179 int b0, b1, b2, k, verbose; 00180 static double unopargs[] = { .7, -.9, 2.5, -4.2, .1, .3, 1.2, -2.4, 0. }; 00181 static double binargs[] = { .1,.2, -.3,.4, -.5,-.6, .7,-.8, 00182 1.1,2.2, -2.3,3.4, -2.5,-3.6, 2.7,-3.8, 0.}; 00183 00184 #define Dcl(T,T1) T x ## T1, y ## T1, z ## T1 00185 Dcl(ADVar, ADV); 00186 Dcl(DADVar, DADV); 00187 Dcl(FDReal, FDR); 00188 Dcl(AFReal, AFR); 00189 #undef Dcl 00190 ADVar *pADV[2]; 00191 DADVar *pDADV[2]; 00192 DReal tDR; 00193 FDReal tfFDR; 00194 FReal tfFR; 00195 00196 progname = argv[0]; 00197 verbose = 0; 00198 if (argc > 1) { 00199 if (argc > 2 || *(s = argv[1]) != '-') 00200 return usage(1); 00201 switch(s[1]) { 00202 case 'h': 00203 case '?': 00204 return usage(s[2] != 0); 00205 case '-': 00206 if (!s[2]) 00207 break; 00208 return usage(strcmp(s,"--help") ? 1 : 0); 00209 case 'v': 00210 if (!s[2]) { 00211 verbose = 1; 00212 break; 00213 } 00214 default: 00215 return usage(1); 00216 } 00217 } 00218 00219 v[0] = v[2] = 1.; 00220 v[1] = 0.; 00221 #define In(T) p ## T[0] = &x ## T; p ## T[1] = &z ## T 00222 In(ADV); 00223 In(DADV); 00224 #undef In 00225 b0 = b1 = b2 = 0; 00226 fe = UT + sizeof(UT)/sizeof(UT[0]); 00227 for(ap = unopargs; (a = *ap++) != 0.; ) 00228 for(f = UT; f < fe; f++) { 00229 if (a < f->dom[0] || a > f->dom[1]) 00230 continue; 00231 xADV = a; 00232 yADV = f->f1(xADV); 00233 ADVar::Gradcomp(); 00234 ADVar::Hvprod(1,pADV,v,h); 00235 if (verbose) { 00236 printf("%s(%g) = %g\n", f->name, xADV.val(), yADV.val()); 00237 printf("%s' = %g\n", f->name, xADV.adj()); 00238 printf("%s\" = %g\n", f->name, h[0]); 00239 } 00240 xDADV = a; 00241 yDADV = f->f2(xDADV); 00242 DADVar::Gradcomp(); 00243 if (differ(yDADV.val(), yADV.val())) { 00244 ++b0; 00245 printf("%s(%g): yADV = %g, yDADV = %g, diff = %.2g\n", 00246 f->name, a, yADV.val(), yDADV.val(), 00247 yADV.val() - yDADV.val()); 00248 } 00249 if (differ(xADV.adj(), xDADV.adj())) { 00250 ++b1; 00251 printf("%s'(%g): ADV ==> %g, DADV ==> %g, diff = %.2g\n", 00252 f->name, a, xADV.adj(), xDADV.adj(), 00253 xADV.adj() - xDADV.adj()); 00254 } 00255 DADVar::Hvprod(1,pDADV,v,h1); 00256 if (differ(h[0], h1[0])) { 00257 ++b2; 00258 printf("%s\"(%g): ADV ==> %g, DADV ==> %g, diff = %.2g\n", 00259 f->name, a, h[0], h1[0], h[0] - h1[0]); 00260 } 00261 00262 tfFDR = 0.; 00263 tfFDR.diff(0,1); 00264 xFDR = a + tfFDR*v[0]; 00265 yFDR = f->f3(xFDR); 00266 tDR = yFDR.dx(0); 00267 DReal::Gradcomp(); 00268 00269 if (differ(yFDR.val().val(), yADV.val())) { 00270 ++b0; 00271 printf("%s(%g): yFDR = %g, yADV = %g, diff = %.2g\n", 00272 f->name, a, yFDR.val().val(), yADV.val(), 00273 yFDR.val().val() - yADV.val()); 00274 } 00275 if (differ(xFDR.val().adj(), h[0])) { 00276 ++b2; 00277 printf("%s\"(%g): FDR ==> %g, ADV ==> %g, diff = %.2g\n", 00278 f->name, a, xFDR.val().adj(), h[0], 00279 xFDR.val().adj() - h[0]); 00280 } 00281 00282 tfFR = 0.; 00283 tfFR.diff(0,1); 00284 xAFR = a + tfFR*v[0]; 00285 yAFR = f->f4(xAFR); 00286 AFReal::Gradcomp(); 00287 if (differ(yAFR.val().val(), yADV.val())) { 00288 ++b0; 00289 printf("%s(%g): yAFR = %g, yADV = %g, diff = %.2g\n", 00290 f->name, a, yAFR.val().val(), yADV.val(), 00291 yAFR.val().val() - yADV.val()); 00292 } 00293 if (differ(xAFR.adj().val(), xADV.adj())) { 00294 ++b1; 00295 printf("%s'(%g): AFR ==> %g, ADV ==> %g, diff = %.2g\n", 00296 f->name, a, xAFR.adj().val(), xADV.adj(), 00297 xAFR.adj().val() - xADV.adj()); 00298 } 00299 if (differ(xAFR.adj().dx(0), h[0])) { 00300 ++b2; 00301 printf("%s\"(%g): AFR ==> %g, ADV ==> %g, diff = %.2g\n", 00302 f->name, a, xAFR.adj().dx(0), h[0], 00303 xAFR.adj().dx(0) - h[0]); 00304 } 00305 } 00306 f2e = UT2 + sizeof(UT2)/sizeof(UT2[0]); 00307 for(ap = binargs; (a = *ap) != 0.; ap += 2) { 00308 c = ap[1]; 00309 for(f2 = UT2; f2 < f2e; f2++) { 00310 if (a < f2->xdom[0] || a > f2->xdom[1]) 00311 continue; 00312 xADV = a; 00313 zADV = c; 00314 yADV = f2->f1(xADV, zADV); 00315 ADVar::Gradcomp(); 00316 ADVar::Hvprod(2,pADV,v,h); 00317 ADVar::Hvprod(2,pADV,v+1,h+2); 00318 if (verbose) { 00319 printf("%s(%g,%g) = %g\n", f2->name, xADV.val(), zADV.val(), yADV.val()); 00320 printf("%s' = (%g, %g)\n", f2->name, xADV.adj(), zADV.adj()); 00321 printf("%s\" = (%8g\t%g)\n(%8g\t%g)", f2->name, h[0],h[1],h[2],h[3]); 00322 } 00323 xDADV = a; 00324 zDADV = c; 00325 yDADV = f2->f2(xDADV, zDADV); 00326 DADVar::Gradcomp(); 00327 if (differ(yDADV.val(), yADV.val())) { 00328 ++b0; 00329 printf("%s(%g,%g): yADV = %g, yDADV = %g, diff = %.2g\n", 00330 f2->name, a, c, yADV.val(), yDADV.val(), 00331 yADV.val() - yDADV.val()); 00332 } 00333 if (differ(xADV.adj(), xDADV.adj()) || differ(zADV.adj(), zDADV.adj())) { 00334 ++b1; 00335 printf("%s'(%g,%g): ADV ==> (%g,%g) DADV ==> (%g,%g), diff = (%.2g,%.2g)\n", 00336 f2->name, a, c, xADV.adj(), zADV.adj(), 00337 xDADV.adj(), zDADV.adj(), 00338 xADV.adj() - xDADV.adj(), zADV.adj() - zDADV.adj()); 00339 } 00340 DADVar::Hvprod(2,pDADV,v,h1); 00341 DADVar::Hvprod(2,pDADV,v+1,h1+2); 00342 for(k = 0; k < 4; k++) 00343 if (differ(h[k], h1[k])) { 00344 ++b2; 00345 printf("%s\"(%g,%g):\n ADV ==> (%8g\t%8g\t%8g\t%g)\n", 00346 f2->name, a, c, h[0], h[1], h[2], h[3]); 00347 printf("DADV ==> (%8g\t%8g\t%8g\t%g)\n", 00348 h1[0], h1[1], h1[2], h1[3]); 00349 printf("diff = (%.2g %.2g %.2g %.2g)\n\n", 00350 h[0] - h1[0], h[1] - h1[1], 00351 h[2] - h1[2], h[3] - h1[3]); 00352 break; 00353 } 00354 tfFDR = 0.; 00355 tfFDR.diff(0,1); 00356 xFDR = a + tfFDR*v[0]; 00357 zFDR = c + tfFDR*v[1]; 00358 yFDR = f2->f3(xFDR, zFDR); 00359 tDR = yFDR.dx(0); 00360 DReal::Gradcomp(); 00361 00362 if (differ(yFDR.val().val(), yADV.val())) { 00363 ++b0; 00364 printf("%s(%g,%g): yFDR = %g, yADV = %g, diff = %.2g\n", 00365 f2->name, a, c, yFDR.val().val(), yADV.val(), 00366 yFDR.val().val() - yADV.val()); 00367 } 00368 if (differ(xFDR.val().adj(), h[0]) || differ(zFDR.val().adj(), h[1])) { 00369 ++b2; 00370 printf("%s\"(%g,%g) row 1:\nFDR ==> (%g,%g)\nADV ==> (%g,%g)\n", 00371 f2->name, a, c, xFDR.val().adj(), zFDR.val().adj(), 00372 h[0], h[1]); 00373 printf("diff = %.2g %g\n\n", 00374 xFDR.val().adj() - h[0], zFDR.val().adj() - h[1]); 00375 } 00376 tfFDR.diff(0,1); 00377 xFDR = a + tfFDR*v[1]; 00378 zFDR = c + tfFDR*v[2]; 00379 yFDR = f2->f3(xFDR, zFDR); 00380 tDR = yFDR.dx(0); 00381 DReal::Gradcomp(); 00382 if (differ(xFDR.val().adj(), h[2]) || differ(zFDR.val().adj(), h[3])) { 00383 ++b2; 00384 printf("%s\"(%g,%g) row 2:\nFDR ==> (%g,%g)\nADV ==> (%g,%g)\n", 00385 f2->name, a, c, xFDR.val().adj(), zFDR.val().adj(), 00386 h[2], h[3]); 00387 printf("diff = %.2g %g\n\n", 00388 xFDR.val().adj() - h[2], zFDR.val().adj() - h[3]); 00389 } 00390 00391 tfFR = 0.; 00392 tfFR.diff(0,1); 00393 xAFR = a + tfFR*v[0]; 00394 zAFR = c + tfFR*v[1]; 00395 yAFR = f2->f4(xAFR, zAFR); 00396 AFReal::Gradcomp(); 00397 if (differ(yAFR.val().val(), yADV.val())) { 00398 ++b0; 00399 printf("%s(%g,%g): yAFR = %g, yADV = %g, diff = %.2g\n", 00400 f2->name, a, c, yAFR.val().val(), yADV.val(), 00401 yAFR.val().val() - yADV.val()); 00402 } 00403 if (differ(xAFR.adj().val(), xADV.adj()) || differ(zAFR.adj().val(), zADV.adj())) { 00404 ++b1; 00405 printf("%s'(%g,%g):\n AFR ==> (%g,%g)\n ADV ==> (%g,%g)\n", 00406 f2->name, a, c, xAFR.adj().val(), zAFR.adj().val(), 00407 xADV.adj(), zADV.adj()); 00408 printf(" diff = %.2g %.2g\n\n", 00409 xAFR.adj().val() - xADV.adj(), 00410 zAFR.adj().val() - zADV.adj()); 00411 } 00412 if (differ(xAFR.adj().dx(0), h[0]) || differ(zAFR.adj().dx(0), h[1])) { 00413 ++b2; 00414 printf("%s\"(%g,%g) row 1:\n AFR ==> (%g,%g)\n", 00415 f2->name, a, c, xAFR.adj().dx(0), zAFR.adj().dx(0)); 00416 printf(" ADV = (%g,%g)\n diff = %.2g %.2g\n\n", h[0], h[1], 00417 xAFR.adj().dx(0) - h[0], 00418 zAFR.adj().dx(0) - h[1]); 00419 } 00420 tfFR.diff(0,1); 00421 xAFR = a + tfFR*v[1]; 00422 zAFR = c + tfFR*v[2]; 00423 yAFR = f2->f4(xAFR, zAFR); 00424 AFReal::Gradcomp(); 00425 if (differ(xAFR.adj().dx(0), h[2]) || differ(zAFR.adj().dx(0), h[3])) { 00426 ++b2; 00427 printf("%s\"(%g,%g) row 2:\n AFR ==> (%g,%g)\n", 00428 f2->name, a, c, xAFR.adj().dx(0), zAFR.adj().dx(0)); 00429 printf(" ADV = (%g,%g)\n diff = %.2g %.2g\n\n", h[2], h[3], 00430 xAFR.adj().dx(0) - h[2], 00431 zAFR.adj().dx(0) - h[3]); 00432 } 00433 } 00434 } 00435 k = b0 + b1 + b2; 00436 if (k || verbose) { 00437 s = progname; 00438 if (*s == '.' && s[1] == '/') 00439 for(s += 2; *s == '/'; ++s); 00440 printf("%s: %d errors seen.\n", s, k); 00441 k = 1; 00442 } 00443 else 00444 printf("OK\n"); 00445 return k; 00446 }
1.7.4
