Sacado Package Browser (Single Doxygen Collection) Version of the Day
taylor_expr.cpp
Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 // @HEADER
00031 
00032 #include "Sacado_Random.hpp"
00033 #include "Sacado.hpp"
00034 #include "Sacado_Tay_CacheTaylor.hpp"
00035 
00036 #include "Teuchos_Time.hpp"
00037 #include "Teuchos_CommandLineProcessor.hpp"
00038 
00039 // ADOL-C includes
00040 #ifdef HAVE_ADOLC
00041 #ifdef PACKAGE
00042 #undef PACKAGE
00043 #endif
00044 #ifdef PACKAGE_NAME
00045 #undef PACKAGE_NAME
00046 #endif
00047 #ifdef PACKAGE_BUGREPORT
00048 #undef PACKAGE_BUGREPORT
00049 #endif
00050 #ifdef PACKAGE_STRING
00051 #undef PACKAGE_STRING
00052 #endif
00053 #ifdef PACKAGE_TARNAME
00054 #undef PACKAGE_TARNAME
00055 #endif
00056 #ifdef PACKAGE_VERSION
00057 #undef PACKAGE_VERSION
00058 #endif
00059 #ifdef VERSION
00060 #undef VERSION
00061 #endif
00062 #include "adolc/adouble.h"
00063 #include "adolc/interfaces.h"
00064 #endif
00065 
00066 using std::cout;
00067 using std::endl;
00068 
00069 // A simple performance test that computes a Taylor expansion of a simple
00070 // expression using many variants of Taylor polynomial classes.
00071 
00072 template <typename T>
00073 inline void
00074 func(const T& x1, const T& x2, T& y) {
00075   y = x1*x2 + sin(x1)/x2;
00076 }
00077 
00078 template <typename TaylorType>
00079 double
00080 do_time(int degree, int nloop)
00081 {
00082   TaylorType x1, x2, y;
00083   Sacado::Random<double> urand(0.0, 1.0);
00084 
00085   x1 = TaylorType(degree, urand.number());
00086   x2 = TaylorType(degree, urand.number());
00087   y = 0.0;
00088   for (int j=0; j<=degree; j++) {
00089     x1.fastAccessCoeff(j) = urand.number();
00090     x2.fastAccessCoeff(j) = urand.number();
00091   }
00092   
00093   Teuchos::Time timer("mult", false);
00094   timer.start(true);
00095   for (int j=0; j<nloop; j++) {
00096     func(x1, x2, y);
00097   }
00098   timer.stop();
00099 
00100   return timer.totalElapsedTime() / nloop;
00101 }
00102 
00103 #ifdef HAVE_ADOLC
00104 double
00105 do_time_adolc(int degree, int nloop)
00106 {
00107   Sacado::Random<double> urand(0.0, 1.0);
00108   double **X, **Y;
00109 
00110   X = new double*[2];
00111   X[0] = new double[degree+1];
00112   X[1] = new double[degree+1];
00113   Y = new double*[1];
00114   Y[0] = new double[degree+1]; 
00115   for (int j=0; j<=degree; j++) {
00116     X[0][j] = urand.number();
00117     X[1][j] = urand.number();
00118     Y[0][j] = 0.0;
00119   }
00120 
00121   trace_on(0);
00122   adouble x1, x2, y;
00123   x1 <<= X[0][0];
00124   x2 <<= X[1][0];
00125   func(x1, x2, y);
00126   y >>= Y[0][0];
00127   trace_off();
00128   
00129   Teuchos::Time timer("mult", false);
00130   timer.start(true);
00131   for (int j=0; j<nloop; j++) {
00132     forward(0,1,2,degree,0,X,Y);
00133   }
00134   timer.stop();
00135 
00136   delete [] X[1];
00137   delete [] X[0];
00138   delete [] X;
00139 
00140   delete [] Y[0];
00141   delete [] Y;
00142 
00143   return timer.totalElapsedTime() / nloop;
00144 }
00145 
00146 double
00147 do_time_adolc_retape(int degree, int nloop)
00148 {
00149   Sacado::Random<double> urand(0.0, 1.0);
00150   double **X, **Y;
00151 
00152   X = new double*[2];
00153   X[0] = new double[degree+1];
00154   X[1] = new double[degree+1];
00155   Y = new double*[1];
00156   Y[0] = new double[degree+1]; 
00157   for (int j=0; j<=degree; j++) {
00158     X[0][j] = urand.number();
00159     X[1][j] = urand.number();
00160     Y[0][j] = 0.0;
00161   }
00162   adouble x1, x2, y;
00163   
00164   Teuchos::Time timer("mult", false);
00165   timer.start(true);
00166   for (int j=0; j<nloop; j++) {
00167     trace_on(0);
00168     x1 <<= X[0][0];
00169     x2 <<= X[1][0];
00170     func(x1, x2, y);
00171     y >>= Y[0][0];
00172     trace_off();
00173     forward(0,1,2,degree,0,X,Y);
00174   }
00175   timer.stop();
00176 
00177   delete [] X[1];
00178   delete [] X[0];
00179   delete [] X;
00180 
00181   delete [] Y[0];
00182   delete [] Y;
00183 
00184   return timer.totalElapsedTime() / nloop;
00185 }
00186 #endif
00187 
00188 int main(int argc, char* argv[]) {
00189   int ierr = 0;
00190 
00191   try {
00192     double t;
00193     int p = 2;
00194     int w = p+7;
00195 
00196     // Set up command line options
00197     Teuchos::CommandLineProcessor clp;
00198     clp.setDocString("This program tests the speed of various forward mode AD implementations for a single multiplication operation");
00199     int degree = 10;
00200     clp.setOption("degree", &degree, "Polynomial degree");
00201     int nloop = 1000000;
00202     clp.setOption("nloop", &nloop, "Number of loops");
00203 
00204     // Parse options
00205     Teuchos::CommandLineProcessor::EParseCommandLineReturn
00206       parseReturn= clp.parse(argc, argv);
00207     if(parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
00208       return 1;
00209 
00210     std::cout.setf(std::ios::scientific);
00211     std::cout.precision(p);
00212     std::cout << "Times (sec) for degree = " << degree
00213         << " nloop =  " << nloop << ":  " << std::endl;
00214 
00215     t = do_time< Sacado::Tay::Taylor<double> >(degree, nloop);
00216     std::cout << "Taylor:          " << std::setw(w) << t << std::endl;
00217 
00218     t = do_time< Sacado::Tay::CacheTaylor<double> >(degree, nloop);
00219     std::cout << "CacheTaylor:     " << std::setw(w) << t << std::endl;
00220 
00221 #ifdef HAVE_ADOLC
00222     t = do_time_adolc(degree, nloop);
00223     std::cout << "ADOL-C:          " << std::setw(w) << t << std::endl;
00224     t = do_time_adolc_retape(degree, nloop);
00225     std::cout << "ADOL-C (retape): " << std::setw(w) << t << std::endl;
00226 #endif
00227     
00228   }
00229   catch (std::exception& e) {
00230     cout << e.what() << endl;
00231     ierr = 1;
00232   }
00233   catch (const char *s) {
00234     cout << s << endl;
00235     ierr = 1;
00236   }
00237   catch (...) {
00238     cout << "Caught unknown exception!" << endl;
00239     ierr = 1;
00240   }
00241 
00242   return ierr;
00243 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines