ad_example.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //                           Sacado Package
00005 //                 Copyright (2006) 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 // Brief demo of Fad and Rad for computing first derivatives
00031 
00032 #include <Sacado.hpp>   // for FAD and RAD
00033 #include <cstdio>   // nicer than streams in some respects
00034 using std::printf;
00035 
00036 using namespace std;
00037 
00038 // Typedefs reduce gobbledygook below
00039 
00040 typedef Sacado::Fad::DFad<double>   F;  // FAD with # of ind. vars given later
00041 typedef Sacado::Fad::SFad<double,2> F2; // FAD with # of ind. vars fixed at 2
00042 typedef Sacado::Rad::ADvar<double>  R;  // for RAD
00043 
00044 template <typename T>
00045 const T func2(T &a, T &b) // sample function of 2 variables
00046 { return sqrt(a*a + 2*b*b); }
00047 
00048 template <typename T>
00049 const T func(int n, T *x) // sample function of n variables
00050         // == func2 when n == 2
00051 {
00052   int i;
00053   T t = 0;
00054   for(i = 1; i < n; i++)
00055     t += i*x[i]*x[i];
00056   return sqrt(t);
00057   }
00058 
00059 // demo of FAD (forward-mode AD), with general number of ind. vars
00060 
00061  void
00062 Fad_demo()
00063 {
00064   F a, b, x[5], y;
00065   int i, n;
00066 
00067   printf("Fad_demo...\n\n");
00068 
00069   // first try n == 2
00070   a = 1.;
00071   b = 2.;
00072   // indicate the independent variables, and initialize their partials to 1:
00073   a.diff(0,2);  // 0 ==> this is the first independent var., of 2
00074   b.diff(1,2);  // 1 ==> this is the second ind. var.
00075   
00076   y = func2(a,b);
00077 
00078   printf("func2(%g,%g) = %g\n", a.val(), b.val(), y.val());
00079 
00080   printf("partials of func2 = %g, %g\n", y.dx(0), y.dx(1));
00081 
00082   // When we know the result was not constant (i.e., did involve ind. vars)
00083   // or when hasFastAccess() is true, we access partials more quickly
00084   // by using member function fastAccessDx rather than dx
00085 
00086   if (y.hasFastAccess())
00087     printf("Repeat with fastAccess: partials of func2 = %g, %g\n",
00088       y.fastAccessDx(0), y.fastAccessDx(1));
00089 
00090   // Similar exercise with general n, in this case n == 5
00091   n = 5;
00092   for(i = 0; i < n; i++) {
00093     x[i] = i;
00094     x[i].diff(i, n);
00095     }
00096   y = func(n, x);
00097   printf("\nfunc(5,x) for x = (0,1,2,3,4) = %g\n", y.val());
00098   for(i = 0; i < n; i++)
00099     printf("d func / d x[%d] = %g == %g\n", i, y.dx(i), y.fastAccessDx(i));
00100   }
00101 
00102 // Fad_demo2 == repeat first part of Fad_Demo with type F2 instead of F
00103 // i.e., with fixed-size allocations
00104 
00105  void
00106 Fad2_demo()
00107 {
00108   F2 a, b, y;
00109 
00110   printf("\n\nFad2_demo...\n\n");
00111 
00112   a = 1.;
00113   b = 2.;
00114   // indicate the independent variables, and initialize their partials to 1:
00115   a.diff(0,2);  // 0 ==> this is the first independent var., of 2
00116   b.diff(1,2);  // 1 ==> this is the second ind. var.
00117   
00118   y = func2(a,b);
00119 
00120   printf("func2(%g,%g) = %g\n", a.val(), b.val(), y.val());
00121 
00122   printf("partials of func2 = %g, %g\n", y.dx(0), y.dx(1));
00123 
00124   if (y.hasFastAccess())
00125     printf("Repeat with fastAccess: partials of func2 = %g, %g\n",
00126       y.fastAccessDx(0), y.fastAccessDx(1));
00127   }
00128 
00129 // Fad_demo3 == repeat of Fad_Demo2 with a different constructor, one that
00130 // indicates the independent variables and their initial values
00131 // and removes the need to invoke .diff()
00132 
00133  void
00134 Fad3_demo()
00135 {
00136   F2 a(2,0,1.), b(2,1,2.), y;
00137 
00138   printf("\n\nFad3_demo...\n\n");
00139 
00140   y = func2(a,b);
00141 
00142   printf("func2(%g,%g) = %g\n", a.val(), b.val(), y.val());
00143 
00144   printf("partials of func2 = %g, %g\n", y.dx(0), y.dx(1));
00145 
00146   if (y.hasFastAccess())
00147     printf("Repeat with fastAccess: partials of func2 = %g, %g\n",
00148       y.fastAccessDx(0), y.fastAccessDx(1));
00149   }
00150 
00151 // Rad_demo == repeat of Fad_Demo with type R instead of F,
00152 // i.e., with reverse-mode rather than forward-mode AD
00153 
00154  void
00155 Rad_demo()
00156 {
00157   R a, b, x[5], y;
00158   int i, n;
00159 
00160   printf("\n\nRad_demo...\n\n");
00161 
00162   // first try n == 2
00163   a = 1.;
00164   b = 2.;
00165   
00166   y = func2(a,b);
00167 
00168   R::Gradcomp();  // do the reverse sweep
00169 
00170   printf("func2(%g,%g) = %g\n", a.val(), b.val(), y.val());
00171 
00172   printf("partials of func2 = %g, %g\n", a.adj(), b.adj());
00173 
00174   // Similar exercise with general n, in this case n == 5
00175   n = 5;
00176   for(i = 0; i < n; i++)
00177     x[i] = i;
00178   y = func(n, x);
00179   printf("\nfunc(5,x) for x = (0,1,2,3,4) = %g\n", y.val());
00180 
00181   // the .val() values are always available; we must call Gradcomp
00182   // before accessing the adjoints, i.e., the .adj() values...
00183 
00184   R::Gradcomp();
00185   
00186   for(i = 0; i < n; i++)
00187     printf("d func / d x[%d] = %g\n", i, x[i].adj());
00188   }
00189 
00190  int
00191 main()
00192 {
00193   Fad_demo();
00194   Fad2_demo();
00195   Fad3_demo();
00196   Rad_demo();
00197   return 0;
00198   }

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