hesopcheck.cpp

Go to the documentation of this file.
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   }

Generated on Wed May 12 21:39:33 2010 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.4.7