Sacado Package Browser (Single Doxygen Collection) Version of the Day
fad_fe_jac_fill.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 #include "Fad/fad.h"
00035 #include "TinyFadET/tfad.h"
00036 
00037 // A performance test that computes a finite-element-like Jacobian using
00038 // several Fad variants
00039 
00040 template <>
00041 Sacado::Fad::MemPool* Sacado::Fad::MemPoolStorage<double>::defaultPool_ = NULL;
00042 
00043 void FAD::error(const char *msg) {
00044   std::cout << msg << std::endl;
00045 }
00046 
00047 int main(int argc, char* argv[]) {
00048   int ierr = 0;
00049 
00050   try {
00051     double t, ta, tr;
00052     int p = 2;
00053     int w = p+7;
00054 
00055     // Maximum number of derivative components for SLFad
00056     const int slfad_max = 130;
00057 
00058     // Set up command line options
00059     Teuchos::CommandLineProcessor clp;
00060     clp.setDocString("This program tests the speed of various forward mode AD implementations for a finite-element-like Jacobian fill");
00061     int num_nodes = 100000;
00062     int num_eqns = 2;
00063     int rt = 0;
00064     clp.setOption("n", &num_nodes, "Number of nodes");
00065     clp.setOption("p", &num_eqns, "Number of equations");
00066     clp.setOption("rt", &rt, "Include ADOL-C retaping test");
00067 
00068     // Parse options
00069     Teuchos::CommandLineProcessor::EParseCommandLineReturn
00070       parseReturn= clp.parse(argc, argv);
00071     if(parseReturn != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
00072       return 1;
00073 
00074     double mesh_spacing = 1.0 / static_cast<double>(num_nodes - 1);
00075 
00076     // Memory pool & manager
00077     Sacado::Fad::MemPoolManager<double> poolManager(num_nodes*num_eqns);
00078     Sacado::Fad::MemPool* pool = poolManager.getMemoryPool(num_nodes*num_eqns);
00079     Sacado::Fad::DMFad<double>::setDefaultPool(pool);
00080 
00081     std::cout.setf(std::ios::scientific);
00082     std::cout.precision(p);
00083     std::cout << "num_nodes =  " << num_nodes 
00084         << ", num_eqns = " << num_eqns << ":  " << std::endl
00085         << "               " << "   Time   " << "\t"<< "Time/Analytic" << "\t"
00086         << "Time/(2*p*Residual)" << std::endl;
00087 
00088     ta = 1.0;
00089     tr = 1.0;
00090 
00091     tr = residual_fill(num_nodes, num_eqns, mesh_spacing);
00092 
00093     ta = analytic_jac_fill(num_nodes, num_eqns, mesh_spacing);
00094     std::cout << "Analytic:      " << std::setw(w) << ta << "\t" << std::setw(w) << ta/ta << "\t" << std::setw(w) << ta/(2.0*num_eqns*tr) << std::endl;
00095 
00096 #ifdef HAVE_ADOLC
00097 #ifndef ADOLC_TAPELESS
00098     t = adolc_jac_fill(num_nodes, num_eqns, mesh_spacing);
00099     std::cout << "ADOL-C:        " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00100 
00101     if (rt != 0) {
00102       t = adolc_retape_jac_fill(num_nodes, num_eqns, mesh_spacing);
00103       std::cout << "ADOL-C(rt):  " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00104     }
00105 
00106 #else
00107     t = adolc_tapeless_jac_fill(num_nodes, num_eqns, mesh_spacing);
00108     std::cout << "ADOL-C(tl):    " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00109 #endif
00110 #endif
00111 
00112 #ifdef HAVE_ADIC
00113     t = adic_jac_fill(num_nodes, num_eqns, mesh_spacing);
00114     std::cout << "ADIC:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00115 #endif
00116 
00117     if (num_eqns*2 == 4) {
00118       t = fad_jac_fill< FAD::TFad<16,double> >(num_nodes, num_eqns, mesh_spacing);
00119       std::cout << "TFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00120     }
00121     else if (num_eqns*2 == 16) {
00122       t = fad_jac_fill< FAD::TFad<16,double> >(num_nodes, num_eqns, mesh_spacing);
00123       std::cout << "TFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00124     }
00125     else if (num_eqns*2 == 32) {
00126       t = fad_jac_fill< FAD::TFad<32,double> >(num_nodes, num_eqns, mesh_spacing);
00127       std::cout << "TFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00128     }
00129     else if (num_eqns*2 == 64) {
00130       t = fad_jac_fill< FAD::TFad<64,double> >(num_nodes, num_eqns, mesh_spacing);
00131       std::cout << "TFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00132     }
00133 
00134     t = fad_jac_fill< FAD::Fad<double> >(num_nodes, num_eqns, mesh_spacing);
00135     std::cout << "Fad:           " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00136 
00137     if (num_eqns*2 == 4) {
00138       t = fad_jac_fill< Sacado::Fad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
00139       std::cout << "SFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00140     }
00141     else if (num_eqns*2 == 16) {
00142       t = fad_jac_fill< Sacado::Fad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
00143       std::cout << "SFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00144     }
00145     else if (num_eqns*2 == 32) {
00146       t = fad_jac_fill< Sacado::Fad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
00147       std::cout << "SFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00148     }
00149     else if (num_eqns*2 == 64) {
00150       t = fad_jac_fill< Sacado::Fad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
00151       std::cout << "SFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00152     }
00153 
00154     if (num_eqns*2 < slfad_max) {
00155       t = fad_jac_fill< Sacado::Fad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
00156       std::cout << "SLFad:         " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00157     }
00158     
00159     t = fad_jac_fill< Sacado::Fad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
00160     std::cout << "DFad:          " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00161 
00162     t = fad_jac_fill< Sacado::Fad::SimpleFad<double> >(num_nodes, num_eqns, mesh_spacing);
00163     std::cout << "SimpleFad:     " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00164 
00165     t = fad_jac_fill< Sacado::Fad::DMFad<double> >(num_nodes, num_eqns, mesh_spacing);
00166     std::cout << "DMFad:         " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl; 
00167 
00168     if (num_eqns*2 == 4) {
00169       t = fad_jac_fill< Sacado::ELRFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
00170       std::cout << "ELRSFad:       " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00171     }
00172     else if (num_eqns*2 == 16) {
00173       t = fad_jac_fill< Sacado::ELRFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
00174       std::cout << "ELRSFad:       " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00175     }
00176     else if (num_eqns*2 == 32) {
00177       t = fad_jac_fill< Sacado::ELRFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
00178       std::cout << "ELRSFad:       " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00179     }
00180     else if (num_eqns*2 == 64) {
00181       t = fad_jac_fill< Sacado::ELRFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
00182       std::cout << "ELRSFad:       " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00183     }
00184 
00185     if (num_eqns*2 < slfad_max) {
00186       t = fad_jac_fill< Sacado::ELRFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
00187       std::cout << "ELRSLFad:      " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00188     }
00189 
00190     t = fad_jac_fill< Sacado::ELRFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
00191     std::cout << "ELRDFad:       " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00192 
00193     if (num_eqns*2 == 4) {
00194       t = fad_jac_fill< Sacado::CacheFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
00195       std::cout << "CacheSFad:     " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00196     }
00197     else if (num_eqns*2 == 16) {
00198       t = fad_jac_fill< Sacado::CacheFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
00199       std::cout << "CacheSFad:     " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00200     }
00201     else if (num_eqns*2 == 32) {
00202       t = fad_jac_fill< Sacado::CacheFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
00203       std::cout << "CacheSFad:     " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00204     }
00205     else if (num_eqns*2 == 64) {
00206       t = fad_jac_fill< Sacado::CacheFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
00207       std::cout << "CacheSFad:     " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00208     }
00209 
00210     if (num_eqns*2 < slfad_max) {
00211       t = fad_jac_fill< Sacado::CacheFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
00212       std::cout << "CacheSLFad:    " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00213     }
00214     
00215     t = fad_jac_fill< Sacado::CacheFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
00216     std::cout << "CacheFad:      " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00217 
00218     if (num_eqns*2 == 4) {
00219       t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,4> >(num_nodes, num_eqns, mesh_spacing);
00220       std::cout << "ELRCacheSFad:  " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00221     }
00222     else if (num_eqns*2 == 16) {
00223       t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,16> >(num_nodes, num_eqns, mesh_spacing);
00224       std::cout << "ELRCacheSFad:  " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00225     }
00226     else if (num_eqns*2 == 32) {
00227       t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,32> >(num_nodes, num_eqns, mesh_spacing);
00228       std::cout << "ELRCacheSFad:  " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00229     }
00230     else if (num_eqns*2 == 64) {
00231       t = fad_jac_fill< Sacado::ELRCacheFad::SFad<double,64> >(num_nodes, num_eqns, mesh_spacing);
00232       std::cout << "ELRCacheSFad:  " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00233     }
00234 
00235     if (num_eqns*2 < slfad_max) {
00236       t = fad_jac_fill< Sacado::ELRCacheFad::SLFad<double,slfad_max> >(num_nodes, num_eqns, mesh_spacing);
00237       std::cout << "ELRCacheSLFad: " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00238     }
00239 
00240     t = fad_jac_fill< Sacado::ELRCacheFad::DFad<double> >(num_nodes, num_eqns, mesh_spacing);
00241     std::cout << "ELRCacheFad:   " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00242 
00243     t = fad_jac_fill< Sacado::Fad::DVFad<double> >(num_nodes, num_eqns, mesh_spacing);
00244     std::cout << "DVFad:         " << std::setw(w) << t << "\t" << std::setw(w) << t/ta << "\t" << std::setw(w) << t/(2.0*num_eqns*tr) << std::endl;
00245     
00246   }
00247   catch (std::exception& e) {
00248     cout << e.what() << endl;
00249     ierr = 1;
00250   }
00251   catch (const char *s) {
00252     cout << s << endl;
00253     ierr = 1;
00254   }
00255   catch (...) {
00256     cout << "Caught unknown exception!" << endl;
00257     ierr = 1;
00258   }
00259 
00260   return ierr;
00261 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines