Sacado Package Browser (Single Doxygen Collection) Version of the Day
fad_fe_jac_fill_range.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 "fe_jac_fill_funcs.hpp"
00033 
00034 std::vector<double> 
00035 do_times(int work_count, int num_eqns_begin, int num_eqns_end,
00036    int num_eqns_delta,
00037    double (*func)(unsigned int,unsigned int,double)) {
00038   std::vector<double> times;
00039   for (int num_eqns = num_eqns_begin; num_eqns <= num_eqns_end; 
00040        num_eqns += num_eqns_delta) {
00041     int num_nodes = work_count / num_eqns;
00042     double mesh_spacing = 1.0 / (num_nodes - 1);
00043     times.push_back(func(num_nodes, num_eqns, mesh_spacing));
00044   }
00045   return times;
00046 }
00047 
00048 template <template <typename> class FadType>
00049 std::vector<double> 
00050 do_times_fad(int work_count, int num_eqns_begin, int num_eqns_end,
00051        int num_eqns_delta) {
00052   std::vector<double> times;
00053   for (int num_eqns = num_eqns_begin; num_eqns <= num_eqns_end; 
00054        num_eqns += num_eqns_delta) {
00055     int num_nodes = work_count / num_eqns;
00056     double mesh_spacing = 1.0 / (num_nodes - 1);
00057     times.push_back(fad_jac_fill<FadType<double> >(num_nodes, num_eqns, mesh_spacing));
00058   }
00059   return times;
00060 }
00061 
00062 template <template <typename,int> class FadType>
00063 std::vector<double> 
00064 do_times_slfad(int work_count, int num_eqns_begin, int num_eqns_end,
00065        int num_eqns_delta) {
00066   const int slfad_max = 130; // Maximum number of derivative components for SLFad
00067   std::vector<double> times;
00068   for (int num_eqns = num_eqns_begin; num_eqns <= num_eqns_end; 
00069        num_eqns += num_eqns_delta) {
00070     int num_nodes = work_count / num_eqns;
00071     double mesh_spacing = 1.0 / (num_nodes - 1);
00072     if (num_eqns*2 < slfad_max)
00073       times.push_back(fad_jac_fill<FadType<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing));
00074   }
00075   return times;
00076 }
00077 
00078 template <template <typename,int> class FadType>
00079 std::vector<double> 
00080 do_times_sfad(int work_count, int num_eqns_begin, int num_eqns_end,
00081        int num_eqns_delta) {
00082   std::vector<double> times;
00083   for (int num_eqns = num_eqns_begin; num_eqns <= num_eqns_end; 
00084        num_eqns += num_eqns_delta) {
00085     int num_nodes = work_count / num_eqns;
00086     double mesh_spacing = 1.0 / (num_nodes - 1);
00087     if (num_eqns*2 == 10)
00088       times.push_back(fad_jac_fill<FadType<double,10> >(num_nodes, num_eqns, mesh_spacing));
00089     else if (num_eqns*2 == 30)
00090       times.push_back(fad_jac_fill<FadType<double,30> >(num_nodes, num_eqns, mesh_spacing));
00091     else if (num_eqns*2 == 50)
00092       times.push_back(fad_jac_fill<FadType<double,50> >(num_nodes, num_eqns, mesh_spacing));
00093     else if (num_eqns*2 == 70)
00094       times.push_back(fad_jac_fill<FadType<double,70> >(num_nodes, num_eqns, mesh_spacing));
00095     else if (num_eqns*2 == 90)
00096       times.push_back(fad_jac_fill<FadType<double,90> >(num_nodes, num_eqns, mesh_spacing));
00097     else if (num_eqns*2 == 110)
00098       times.push_back(fad_jac_fill<FadType<double,110> >(num_nodes, num_eqns, mesh_spacing));
00099     else if (num_eqns*2 == 130)
00100       times.push_back(fad_jac_fill<FadType<double,130> >(num_nodes, num_eqns, mesh_spacing));
00101   }
00102   return times;
00103 }
00104 
00105 void print_times(const std::vector<double>& times,
00106      const std::vector<double>& base,
00107      const std::string& name, int p, int w, int w_name) {
00108   std::cout.setf(std::ios::scientific);
00109   std::cout.precision(p);
00110   std::cout.setf(std::ios::right);
00111   std::cout << std::setw(w_name) << name << " ";
00112   std::cout.setf(std::ios::right);
00113   for (unsigned int i=0; i<times.size(); i++)
00114     std::cout << std::setw(w) << times[i]/base[i] << " ";
00115   std::cout << std::endl;
00116 }
00117 
00118 int main(int argc, char* argv[]) {
00119   int ierr = 0;
00120   int p = 1;
00121   int w = p+7;
00122   int w_name = 13;
00123 
00124   try {
00125  
00126     // Set up command line options
00127     Teuchos::CommandLineProcessor clp;
00128     clp.setDocString("This program tests the speed of various forward mode AD implementations for a finite-element-like Jacobian fill");
00129     int work_count = 200000;
00130     int num_eqns_begin = 5;
00131     int num_eqns_end = 65;
00132     int num_eqns_delta = 10;
00133     int rt = 0;
00134     clp.setOption("wc", &work_count, "Work count = num_nodes*num_eqns");
00135     clp.setOption("p_begin", &num_eqns_begin, "Intitial number of equations");
00136     clp.setOption("p_end", &num_eqns_end, "Final number of equations");
00137     clp.setOption("p_delta", &num_eqns_delta, "Step in number of equations");
00138     clp.setOption("rt", &rt, "Include ADOL-C retaping test");
00139 
00140     // Parse options
00141     Teuchos::CommandLineProcessor::EParseCommandLineReturn
00142       parseReturn= clp.parse(argc, argv);
00143     if(parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
00144       return 1;
00145 
00146     // Print header
00147     std::cout.setf(std::ios::right);
00148     std::cout << std::setw(w_name) << "Name" << " ";
00149     for (int num_eqns = num_eqns_begin; num_eqns <= num_eqns_end; 
00150    num_eqns += num_eqns_delta)
00151       std::cout << std::setw(w) << num_eqns << " ";
00152     std::cout << std::endl;
00153     for (int j=0; j<w_name; j++)
00154       std::cout << '=';
00155     std::cout << " ";
00156     for (int num_eqns = num_eqns_begin; num_eqns <= num_eqns_end; 
00157    num_eqns += num_eqns_delta) {
00158       for (int j=0; j<w; j++)
00159   std::cout << '=';
00160       std::cout << " ";
00161     }
00162     std::cout << std::endl;
00163 
00164     // Analytic
00165     std::vector<double> times_analytic =
00166       do_times(work_count, num_eqns_begin, num_eqns_end, num_eqns_delta, 
00167          analytic_jac_fill);
00168     print_times(times_analytic, times_analytic, "Analytic", p, w, w_name);
00169 
00170 #ifdef HAVE_ADIC
00171     // Note there seems to be a bug in ADIC where doing more than one num_eqns
00172     // value results in incorrect timings after the first.  Doing one value
00173     // at a time seems to give correct values though.
00174     std::vector<double> times_adic =
00175       do_times(work_count, num_eqns_begin, num_eqns_end, num_eqns_delta, 
00176          adic_jac_fill);
00177     print_times(times_adic, times_analytic, "ADIC", p, w, w_name);
00178 #endif
00179 
00180     // Original Fad
00181     std::vector<double> times_sfad =
00182       do_times_sfad<Sacado::Fad::SFad>(
00183   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00184     print_times(times_sfad, times_analytic, "SFAD", p, w, w_name);
00185 
00186     std::vector<double> times_slfad =
00187       do_times_sfad<Sacado::Fad::SLFad>(
00188   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00189     print_times(times_slfad, times_analytic, "SLFAD", p, w, w_name);
00190 
00191     std::vector<double> times_dfad =
00192       do_times_fad<Sacado::Fad::DFad>(
00193   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00194     print_times(times_dfad, times_analytic, "DFAD", p, w, w_name);
00195     
00196     
00197     // ELR Fad
00198     std::vector<double> times_elr_sfad =
00199       do_times_sfad<Sacado::ELRFad::SFad>(
00200   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00201     print_times(times_elr_sfad, times_analytic, "ELRSFAD", p, w, w_name);
00202 
00203     std::vector<double> times_elr_slfad =
00204       do_times_sfad<Sacado::ELRFad::SLFad>(
00205   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00206     print_times(times_elr_slfad, times_analytic, "ELRSLFAD", p, w, w_name);
00207 
00208     std::vector<double> times_elr_dfad =
00209       do_times_fad<Sacado::ELRFad::DFad>(
00210   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00211     print_times(times_elr_dfad, times_analytic, "ELRDFAD", p, w, w_name);
00212    
00213 
00214     // Cache Fad
00215     std::vector<double> times_cache_sfad =
00216       do_times_sfad<Sacado::CacheFad::SFad>(
00217   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00218     print_times(times_cache_sfad, times_analytic, "CacheSFAD", p, w, w_name);
00219 
00220     std::vector<double> times_cache_slfad =
00221       do_times_sfad<Sacado::CacheFad::SLFad>(
00222   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00223     print_times(times_cache_slfad, times_analytic, "CacheSLFAD", p, w, w_name);
00224 
00225     std::vector<double> times_cache_dfad =
00226       do_times_fad<Sacado::CacheFad::DFad>(
00227   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00228     print_times(times_cache_dfad, times_analytic, "CacheDFAD", p, w, w_name);
00229 
00230     // ELR Cache Fad
00231     std::vector<double> times_cache_elr_sfad =
00232       do_times_sfad<Sacado::ELRCacheFad::SFad>(
00233   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00234     print_times(times_cache_elr_sfad, times_analytic, "ELRCacheSFAD", p, w, w_name);
00235 
00236     std::vector<double> times_cache_elr_slfad =
00237       do_times_sfad<Sacado::ELRCacheFad::SLFad>(
00238   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00239     print_times(times_cache_elr_slfad, times_analytic, "ELRCacheSLFAD", p, w, w_name);
00240 
00241     std::vector<double> times_cache_elr_dfad =
00242       do_times_fad<Sacado::ELRCacheFad::DFad>(
00243   work_count, num_eqns_begin, num_eqns_end, num_eqns_delta);
00244     print_times(times_cache_elr_dfad, times_analytic, "ELRCacheDFAD", p, w, w_name);
00245     
00246   }
00247   catch (std::exception& e) {
00248     std::cout << e.what() << std::endl;
00249     ierr = 1;
00250   }
00251   catch (const char *s) {
00252     std::cout << s << std::endl;
00253     ierr = 1;
00254   }
00255   catch (...) {
00256     std::cout << "Caught unknown exception!" << std::endl;
00257     ierr = 1;
00258   }
00259 
00260   return ierr;
00261 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines