Sacado Package Browser (Single Doxygen Collection) Version of the Day
fad_expr_depth.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 
00034 #include "Teuchos_Time.hpp"
00035 #include "Teuchos_Array.hpp"
00036 #include <fstream>
00037 
00038 #include "fad_expr_funcs.hpp"
00039 
00040 // A simple performance test that computes the derivative of expressions of
00041 // various depths.
00042 
00043 template <typename T, typename F>
00044 double do_time(const T x[], int nloop, const F& f)
00045 {
00046   T y = 0.0;
00047   Teuchos::Time timer("F", false);
00048 
00049   timer.start(true);
00050   for (int j=0; j<nloop; j++)
00051     f(x, y);
00052   timer.stop();
00053   return timer.totalElapsedTime() / nloop;
00054 }
00055 
00056 template <typename T, template <typename, int> class F>
00057 void do_times(const T x[], int nloop, Teuchos::Array<double>& times)
00058 {
00059   int i = 0;
00060   times[i++] = do_time(x, nloop, F<T,1>());
00061   times[i++] = do_time(x, nloop, F<T,2>());
00062   times[i++] = do_time(x, nloop, F<T,3>());
00063   times[i++] = do_time(x, nloop, F<T,4>());
00064   times[i++] = do_time(x, nloop, F<T,5>());
00065   times[i++] = do_time(x, nloop, F<T,10>());
00066   times[i++] = do_time(x, nloop, F<T,15>());
00067   times[i++] = do_time(x, nloop, F<T,20>());
00068 }
00069 
00070 #ifdef HAVE_ADOLC
00071 template <typename F>
00072 double do_time_adolc(double *x, double **seed, int d, int nloop, 
00073          int tag, const F& f)
00074 {
00075   Teuchos::Time timer("F", false);
00076   int n = F::n;
00077   trace_on(tag);  
00078   adouble *x_ad = new adouble[n];
00079   for (int i=0; i<n; i++)
00080     x_ad[i] <<= x[i];
00081   adouble y_ad = 0.0;
00082   f(x_ad, y_ad);
00083   double y;
00084   y_ad >>= y;
00085   delete [] x_ad;
00086   trace_off();
00087 
00088   double **jac = new double*[1];
00089   jac[0] = new double[d];
00090   timer.start(true);
00091   for (int j=0; j<nloop; j++)
00092     forward(tag, 1, n, d, x, seed, &y, jac);
00093   timer.stop();
00094 
00095   delete [] jac[0];
00096   delete [] jac;
00097 
00098   return timer.totalElapsedTime() / nloop;
00099 }
00100 
00101 template <template <typename,int> class F>
00102 void do_times_adolc(double *x, double **seed, int d, int nloop, 
00103         int& tag, Teuchos::Array<double>& times)
00104 {
00105   int i = 0;
00106   times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,1>());
00107   times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,2>());
00108   times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,3>());
00109   times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,4>());
00110   times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,5>());
00111   times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,10>());
00112   times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,15>());
00113   times[i++] = do_time_adolc(x, seed, d, nloop, tag++, F<adouble,20>());
00114 }
00115 #endif
00116 
00117 template <typename FadType>
00118 void print_times(const std::string& screen_name, const std::string& file_name)
00119 {
00120   const int nderiv = 10;
00121   int deriv_dim[nderiv] = { 0, 1, 3, 5, 10, 15, 20, 30, 40, 50 };
00122   Sacado::Random<double> urand(0.0, 1.0);
00123   int p = 1;
00124   int w = p+7;
00125   std::ofstream file(file_name.c_str(), std::ios::out);
00126 
00127   {
00128   std::cout.setf(std::ios::scientific);
00129   std::cout.precision(p);
00130   std::cout << screen_name << " Relative times (time/(func_time*nderiv)): " 
00131       << std::endl;
00132   std::cout << std::setw(5) << "deriv" << " ";
00133   for (int i=0; i<ExprFuncs::nfunc; i++)
00134     std::cout << std::setw(w) << ExprFuncs::mult_names[i] << " ";
00135   std::cout << std::endl;
00136   std::cout << "===== ";
00137   for (int i=0; i<ExprFuncs::nfunc; i++) {
00138     for (int j=0; j<w; j++)
00139       std::cout << '=';
00140     std::cout << " ";
00141   }
00142   std::cout << std::endl;
00143 
00144   // Get function evaluation times
00145   double x[ExprFuncs::nx_max];
00146   for (int i=0; i<ExprFuncs::nx_max; i++)
00147     x[i] = urand.number();
00148   int nloop_func = 10000000;
00149   Teuchos::Array<double> times_func(ExprFuncs::nfunc);
00150   do_times<double,ExprFuncs::mult>(x, nloop_func, times_func);
00151   
00152   // Get times for each derivative dimension
00153   for (int i=0; i<nderiv; i++) {
00154     FadType fx[ExprFuncs::nx_max];
00155     for (int k=0; k<ExprFuncs::nx_max; k++) {
00156       fx[k] = FadType(deriv_dim[i],  urand.number());
00157       for (int j=0; j<deriv_dim[i]; j++) {
00158   fx[k].fastAccessDx(j) = urand.number();
00159       }
00160     }
00161     
00162     //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
00163     int nloop = static_cast<int>(100000.0);
00164     Teuchos::Array<double> times(ExprFuncs::nfunc);
00165     do_times<FadType,ExprFuncs::mult>(fx, nloop, times);
00166      
00167     // Print times
00168     int d = deriv_dim[i];
00169     if (d == 0)
00170       d = 1;
00171     std::cout << std::setw(5) << deriv_dim[i] << " ";
00172     file << deriv_dim[i] << " ";
00173     for (int j=0; j<times.size(); j++) {
00174       double rel_time = times[j]/(times_func[j]*d);
00175       std::cout << std::setw(w) << rel_time << " ";
00176       file << rel_time << " ";
00177     }
00178     std::cout << std::endl;
00179     file << std::endl;
00180   }
00181   }
00182 
00183   {
00184   std::cout.setf(std::ios::scientific);
00185   std::cout.precision(p);
00186   std::cout << screen_name << " Relative times (time/(func_time*nderiv)): " 
00187       << std::endl;
00188   std::cout << std::setw(5) << "deriv" << " ";
00189   for (int i=0; i<ExprFuncs::nfunc; i++)
00190     std::cout << std::setw(w) << ExprFuncs::add_names[i] << " ";
00191   std::cout << std::endl;
00192   std::cout << "===== ";
00193   for (int i=0; i<ExprFuncs::nfunc; i++) {
00194     for (int j=0; j<w; j++)
00195       std::cout << '=';
00196     std::cout << " ";
00197   }
00198   std::cout << std::endl;
00199 
00200   // Get function evaluation times
00201   double x[ExprFuncs::nx_max];
00202   for (int i=0; i<ExprFuncs::nx_max; i++)
00203     x[i] = urand.number();
00204   int nloop_func = 10000000;
00205   Teuchos::Array<double> times_func(ExprFuncs::nfunc);
00206   do_times<double,ExprFuncs::add>(x, nloop_func, times_func);
00207   
00208   // Get times for each derivative dimension
00209   for (int i=0; i<nderiv; i++) {
00210     FadType fx[ExprFuncs::nx_max];
00211     for (int k=0; k<ExprFuncs::nx_max; k++) {
00212       fx[k] = FadType(deriv_dim[i],  urand.number());
00213       for (int j=0; j<deriv_dim[i]; j++) {
00214   fx[k].fastAccessDx(j) = urand.number();
00215       }
00216     }
00217     
00218     //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
00219     int nloop = static_cast<int>(100000.0);
00220     Teuchos::Array<double> times(ExprFuncs::nfunc);
00221     do_times<FadType,ExprFuncs::add>(fx, nloop, times);
00222      
00223     // Print times
00224     int d = deriv_dim[i];
00225     if (d == 0)
00226       d = 1;
00227     std::cout << std::setw(5) << deriv_dim[i] << " ";
00228     file << deriv_dim[i] << " ";
00229     for (int j=0; j<times.size(); j++) {
00230       double rel_time = times[j]/(times_func[j]*d);
00231       std::cout << std::setw(w) << rel_time << " ";
00232       file << rel_time << " ";
00233     }
00234     std::cout << std::endl;
00235     file << std::endl;
00236   }
00237   }
00238 
00239   {
00240   std::cout.setf(std::ios::scientific);
00241   std::cout.precision(p);
00242   std::cout << screen_name << " Relative times (time/(func_time*nderiv)): " 
00243       << std::endl;
00244   std::cout << std::setw(5) << "deriv" << " ";
00245   for (int i=0; i<ExprFuncs::nfunc; i++)
00246     std::cout << std::setw(w) << ExprFuncs::nest_names[i] << " ";
00247   std::cout << std::endl;
00248   std::cout << "===== ";
00249   for (int i=0; i<ExprFuncs::nfunc; i++) {
00250     for (int j=0; j<w; j++)
00251       std::cout << '=';
00252     std::cout << " ";
00253   }
00254   std::cout << std::endl;
00255 
00256   // Get function evaluation times
00257   double x[ExprFuncs::nx_max];
00258   for (int i=0; i<ExprFuncs::nx_max; i++)
00259     x[i] = urand.number();
00260   int nloop_func = 10000000;
00261   Teuchos::Array<double> times_func(ExprFuncs::nfunc);
00262   do_times<double,ExprFuncs::nest>(x, nloop_func, times_func);
00263   
00264   // Get times for each derivative dimension
00265   for (int i=0; i<nderiv; i++) {
00266     FadType fx[ExprFuncs::nx_max];
00267     for (int k=0; k<ExprFuncs::nx_max; k++) {
00268       fx[k] = FadType(deriv_dim[i],  urand.number());
00269       for (int j=0; j<deriv_dim[i]; j++) {
00270   fx[k].fastAccessDx(j) = urand.number();
00271       }
00272     }
00273     
00274     //int nloop = static_cast<int>(1000000.0/(deriv_dim[i]+1));
00275     int nloop = static_cast<int>(100000.0);
00276     Teuchos::Array<double> times(ExprFuncs::nfunc);
00277     do_times<FadType,ExprFuncs::nest>(fx, nloop, times);
00278      
00279     // Print times
00280     int d = deriv_dim[i];
00281     if (d == 0)
00282       d = 1;
00283     std::cout << std::setw(5) << deriv_dim[i] << " ";
00284     file << deriv_dim[i] << " ";
00285     for (int j=0; j<times.size(); j++) {
00286       double rel_time = times[j]/(times_func[j]*d);
00287       std::cout << std::setw(w) << rel_time << " ";
00288       file << rel_time << " ";
00289     }
00290     std::cout << std::endl;
00291     file << std::endl;
00292   }
00293   }
00294 }
00295 
00296 #ifdef HAVE_ADOLC
00297 void print_times_adolc(const std::string& screen_name, 
00298            const std::string& file_name)
00299 {
00300   const int nderiv = 10;
00301   int deriv_dim[nderiv] = { 0, 1, 3, 5, 10, 15, 20, 30, 40, 50 };
00302   const int deriv_max = 50;
00303   Sacado::Random<double> urand(0.0, 1.0);
00304   int p = 1;
00305   int w = p+7;
00306   std::ofstream file(file_name.c_str(), std::ios::out);
00307   int tag = 0;
00308 
00309   double **seed = new double*[ExprFuncs::nx_max];
00310   for (int i=0; i<ExprFuncs::nx_max; i++) {
00311     seed[i] = new double[deriv_max];
00312     for (int j=0; j<deriv_max; j++)
00313       seed[i][j] = urand.number();
00314   }
00315 
00316   {
00317     std::cout.setf(std::ios::scientific);
00318     std::cout.precision(p);
00319     std::cout << screen_name << " Relative times (time/(func_time*nderiv)): " 
00320         << std::endl;
00321     std::cout << std::setw(5) << "deriv" << " ";
00322     for (int i=0; i<ExprFuncs::nfunc; i++)
00323       std::cout << std::setw(w) << ExprFuncs::mult_names[i] << " ";
00324     std::cout << std::endl;
00325     std::cout << "===== ";
00326     for (int i=0; i<ExprFuncs::nfunc; i++) {
00327       for (int j=0; j<w; j++)
00328   std::cout << '=';
00329       std::cout << " ";
00330     }
00331     std::cout << std::endl;
00332     
00333     // Get function evaluation times
00334     double x[ExprFuncs::nx_max];
00335     for (int i=0; i<ExprFuncs::nx_max; i++)
00336       x[i] = urand.number();
00337     int nloop_func = 10000000;
00338     Teuchos::Array<double> times_func(ExprFuncs::nfunc);
00339     do_times<double,ExprFuncs::mult>(x, nloop_func, times_func);
00340     
00341     // Get times for each derivative dimension
00342     for (int i=0; i<nderiv; i++) {
00343       //int nloop = static_cast<int>(100000.0/(deriv_dim[i]+1));
00344       int nloop = static_cast<int>(10000.0);
00345       Teuchos::Array<double> times(ExprFuncs::nfunc);
00346       do_times_adolc<ExprFuncs::mult>(x, seed, deriv_dim[i], nloop, tag, times);
00347       
00348       // Print times
00349       int d = deriv_dim[i];
00350       if (d == 0)
00351   d = 1;
00352       std::cout << std::setw(5) << deriv_dim[i] << " ";
00353       file << deriv_dim[i] << " ";
00354       for (int j=0; j<times.size(); j++) {
00355   double rel_time = times[j]/(times_func[j]*d);
00356   std::cout << std::setw(w) << rel_time << " ";
00357   file << rel_time << " ";
00358       }
00359       std::cout << std::endl;
00360       file << std::endl;
00361     }
00362   }
00363 
00364   {
00365     std::cout.setf(std::ios::scientific);
00366     std::cout.precision(p);
00367     std::cout << screen_name << " Relative times (time/(func_time*nderiv)): " 
00368         << std::endl;
00369     std::cout << std::setw(5) << "deriv" << " ";
00370     for (int i=0; i<ExprFuncs::nfunc; i++)
00371       std::cout << std::setw(w) << ExprFuncs::add_names[i] << " ";
00372     std::cout << std::endl;
00373     std::cout << "===== ";
00374     for (int i=0; i<ExprFuncs::nfunc; i++) {
00375       for (int j=0; j<w; j++)
00376   std::cout << '=';
00377       std::cout << " ";
00378     }
00379     std::cout << std::endl;
00380     
00381     // Get function evaluation times
00382     double x[ExprFuncs::nx_max];
00383     for (int i=0; i<ExprFuncs::nx_max; i++)
00384       x[i] = urand.number();
00385     int nloop_func = 10000000;
00386     Teuchos::Array<double> times_func(ExprFuncs::nfunc);
00387     do_times<double,ExprFuncs::add>(x, nloop_func, times_func);
00388   
00389     // Get times for each derivative dimension
00390      for (int i=0; i<nderiv; i++) {
00391        //int nloop = static_cast<int>(100000.0/(deriv_dim[i]+1));
00392        int nloop = static_cast<int>(10000.0);
00393        Teuchos::Array<double> times(ExprFuncs::nfunc);
00394        do_times_adolc<ExprFuncs::add>(x, seed, deriv_dim[i], nloop, tag, times);
00395      
00396        // Print times
00397        int d = deriv_dim[i];
00398        if (d == 0)
00399    d = 1;
00400        std::cout << std::setw(5) << deriv_dim[i] << " ";
00401        file << deriv_dim[i] << " ";
00402        for (int j=0; j<times.size(); j++) {
00403    double rel_time = times[j]/(times_func[j]*d);
00404    std::cout << std::setw(w) << rel_time << " ";
00405    file << rel_time << " ";
00406        }
00407        std::cout << std::endl;
00408        file << std::endl;
00409      }
00410   }
00411 
00412   {
00413     std::cout.setf(std::ios::scientific);
00414     std::cout.precision(p);
00415     std::cout << screen_name << " Relative times (time/(func_time*nderiv)): " 
00416         << std::endl;
00417     std::cout << std::setw(5) << "deriv" << " ";
00418     for (int i=0; i<ExprFuncs::nfunc; i++)
00419       std::cout << std::setw(w) << ExprFuncs::nest_names[i] << " ";
00420     std::cout << std::endl;
00421     std::cout << "===== ";
00422     for (int i=0; i<ExprFuncs::nfunc; i++) {
00423       for (int j=0; j<w; j++)
00424   std::cout << '=';
00425       std::cout << " ";
00426     }
00427     std::cout << std::endl;
00428     
00429     // Get function evaluation times
00430     double x[ExprFuncs::nx_max];
00431     for (int i=0; i<ExprFuncs::nx_max; i++)
00432       x[i] = urand.number();
00433     int nloop_func = 10000000;
00434     Teuchos::Array<double> times_func(ExprFuncs::nfunc);
00435     do_times<double,ExprFuncs::nest>(x, nloop_func, times_func);
00436     
00437     // Get times for each derivative dimension
00438     for (int i=0; i<nderiv; i++) {
00439       //int nloop = static_cast<int>(10000.0/(deriv_dim[i]+1));
00440       int nloop = static_cast<int>(1000.0);
00441       Teuchos::Array<double> times(ExprFuncs::nfunc);
00442       do_times_adolc<ExprFuncs::nest>(x, seed, deriv_dim[i], nloop, tag, times);
00443       
00444       // Print times
00445       int d = deriv_dim[i];
00446       if (d == 0)
00447   d = 1;
00448       std::cout << std::setw(5) << deriv_dim[i] << " ";
00449       file << deriv_dim[i] << " ";
00450       for (int j=0; j<times.size(); j++) {
00451   double rel_time = times[j]/(times_func[j]*d);
00452   std::cout << std::setw(w) << rel_time << " ";
00453   file << rel_time << " ";
00454       }
00455       std::cout << std::endl;
00456       file << std::endl;
00457     }
00458   }
00459 
00460   delete [] seed;
00461 }
00462 #endif
00463 
00464 int main() {
00465   print_times< Sacado::Fad::DFad<double> >(
00466     "Sacado::Fad::DFad", "fad_expr_depth_dfad.txt");
00467   print_times< Sacado::ELRFad::DFad<double> >(
00468     "Sacado::ELRFad::DFad", "fad_expr_depth_elr_dfad.txt");
00469   print_times< Sacado::CacheFad::DFad<double> >(
00470     "Sacado::CacheFad::DFad", "fad_expr_depth_cache_dfad.txt");
00471   print_times< Sacado::ELRCacheFad::DFad<double> >(
00472     "Sacado::ELRCacheFad::DFad", "fad_expr_depth_elr_cache_dfad.txt");
00473   print_times< Sacado::Fad::SimpleFad<double> >(
00474     "Sacado::Fad::SimpleFad", "fad_expr_depth_simple_fad.txt");
00475 #ifdef HAVE_ADOLC
00476   print_times_adolc("ADOL-C", "fad_expr_depth_adolc.txt");
00477 #endif
00478 
00479   return 0;
00480 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines