Intrepid
http://trilinos.sandia.gov/packages/docs/r11.2/packages/intrepid/test/Discretization/Integration/test_19.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00044 
00051 #include "Intrepid_AdaptiveSparseGrid.hpp"
00052 //#include "Intrepid_CubatureLineSorted.hpp"
00053 #include "Intrepid_Utils.hpp"
00054 #include "Teuchos_oblackholestream.hpp"
00055 #include "Teuchos_RCP.hpp"
00056 #include "Teuchos_RefCountPtr.hpp"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 
00059 using namespace Intrepid;
00060 std::vector<long double> alpha(1,0);
00061 std::vector<long double> beta(1,0);
00062 template<class Scalar>
00063 class StdVector {
00064 private:
00065   Teuchos::RefCountPtr<std::vector<Scalar> >  std_vec_;
00066 
00067 public:
00068 
00069   StdVector( const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec ) 
00070   : std_vec_(std_vec) {}
00071 
00072   Teuchos::RefCountPtr<StdVector<Scalar> > Create() const {
00073     return Teuchos::rcp( new StdVector<Scalar>(
00074            Teuchos::rcp(new std::vector<Scalar>(std_vec_->size(),0))));
00075   }
00076 
00077   void Update( StdVector<Scalar> & s ) {
00078     int dimension  = (int)(std_vec_->size());
00079     for (int i=0; i<dimension; i++)
00080       (*std_vec_)[i] += s[i]; 
00081   }
00082 
00083   void Update( Scalar alpha, StdVector<Scalar> & s )  {
00084     int dimension  = (int)(std_vec_->size());
00085     for (int i=0; i<dimension; i++)
00086       (*std_vec_)[i] += alpha*s[i];
00087   }
00088 
00089   Scalar operator[](int i) {
00090     return (*std_vec_)[i];
00091   }
00092 
00093   void clear() {
00094     std_vec_->clear();
00095   }
00096 
00097   void resize(int n, Scalar p) {
00098     std_vec_->resize(n,p);
00099   }
00100 
00101   int size() {
00102     return (int)std_vec_->size();
00103   }
00104 
00105   void Set( Scalar alpha )  {
00106     int dimension  = (int)(std_vec_->size());
00107     for (int i=0; i<dimension; i++)
00108       (*std_vec_)[i] = alpha;
00109   }
00110 };
00111 
00112 template<class Scalar, class UserVector>
00113 class ASGdata : 
00114   public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> {  
00115 public:  
00116   ~ASGdata() {}
00117 
00118   ASGdata(int dimension,std::vector<EIntrepidBurkardt> rule1D,
00119           std::vector<EIntrepidGrowth> growth1D, int maxLevel,
00120           bool isNormalized) : AdaptiveSparseGridInterface<Scalar,UserVector>(
00121           dimension,rule1D,growth1D,maxLevel,isNormalized) {}
00122 
00123   void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
00124     int    dimension = (int)alpha.size();
00125     Scalar total     = 0.0;
00126     Scalar point     = 0.0;
00127     for (int i=0; i<dimension; i++) {
00128       point     = 0.5*input[i]+0.5;
00129       total    += powl(alpha[i]*(point-beta[i]),(long double)2.0);
00130     }
00131     output.clear(); output.resize(1,std::exp(-total));
00132   }  
00133 
00134   Scalar error_indicator(UserVector & input) {
00135     int dimension = (int)input.size();
00136     Scalar norm2  = 0.0;
00137     for (int i=0; i<dimension; i++)
00138       norm2 += input[i]*input[i];
00139     
00140     Scalar ID = AdaptiveSparseGridInterface<Scalar,UserVector>::
00141       getInitialDiff();
00142     norm2 = std::sqrt(norm2)/ID;
00143     return norm2;
00144   }
00145 };
00146 
00147 long double nCDF(long double z) {
00148   long double p      = 0.0, expntl = 0.0;
00149   long double p0     = 220.2068679123761;
00150   long double p1     = 221.2135961699311;
00151   long double p2     = 112.0792914978709;
00152   long double p3     = 33.91286607838300;
00153   long double p4     = 6.373962203531650;
00154   long double p5     = 0.7003830644436881;
00155   long double p6     = 0.03526249659989109;
00156   long double q0     = 440.4137358247522;
00157   long double q1     = 793.8265125199484;
00158   long double q2     = 637.3336333788311;
00159   long double q3     = 296.5642487796737;
00160   long double q4     = 86.78073220294608;
00161   long double q5     = 16.06417757920695;
00162   long double q6     = 1.755667163182642;
00163   long double q7     = 0.08838834764831844;
00164   long double rootpi = std::sqrt(M_PI);
00165   long double zabs   = fabs(z);
00166 
00167   if (12.0 < zabs) {
00168     p = 0.0;
00169   }
00170   else {
00171     expntl = exp(-zabs*zabs/2.0);
00172     if (zabs < 7.0) {
00173       p = expntl*
00174         ((((((p6*zabs+p5)*zabs+p4)*zabs+p3)*zabs+p2)*zabs+p1)*zabs+p0)/ 
00175         (((((((q7*zabs+q6)*zabs+q5)*zabs+q4)*zabs+q3)*zabs+q2)*zabs+q1)*zabs+q0);
00176     }
00177     else {
00178       p = expntl/(zabs+1.0/(zabs+2.0/(zabs+3.0/(zabs+4.0/(zabs+0.65)))))/rootpi;
00179     }                                
00180   }
00181   if(0.0 < z){
00182     p = 1.0-p;
00183   }
00184   return p;
00185 }
00186 
00187 long double compExactInt(void) {
00188   long double val       = 1.0;
00189   int         dimension = alpha.size();    
00190   long double s2        = std::sqrt(2.0);
00191   long double sp        = std::sqrt(M_PI);
00192   for (int i=0; i<dimension; i++) {
00193     long double s2a = s2*alpha[i];
00194     val *= (sp/alpha[i])*(nCDF((1.0-beta[i])*s2a)-nCDF(-beta[i]*s2a));
00195   }
00196   return val;
00197 }
00198 
00199 long double adaptSG(StdVector<long double> & iv,
00200    AdaptiveSparseGridInterface<long double,StdVector<long double> > & 
00201    problem_data, long double TOL) {
00202 
00203   // Construct a Container for the adapted rule
00204   int dimension = problem_data.getDimension();
00205   std::vector<int> index(dimension,1);
00206   
00207   // Initialize global error indicator
00208   long double eta = 1.0;
00209   
00210   // Initialize the Active index set
00211   std::multimap<long double,std::vector<int> > activeIndex;  
00212   activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index));
00213 
00214   // Initialize the old index set
00215   std::set<std::vector<int> > oldIndex;
00216 
00217   // Perform Adaptation
00218   while (eta > TOL) {
00219     eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(
00220                                                        activeIndex,oldIndex,
00221                                                        iv,eta,problem_data);
00222   }
00223   return eta;
00224 }
00225 
00226 int main(int argc, char *argv[]) {
00227   
00228   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00229 
00230   // This little trick lets us print to std::cout only if
00231   // a (dummy) command-line argument is provided.
00232   int iprint     = argc - 1;
00233   Teuchos::RCP<std::ostream> outStream;
00234   Teuchos::oblackholestream bhs; // outputs nothing
00235   if (iprint > 0)
00236     outStream = Teuchos::rcp(&std::cout, false);
00237   else
00238     outStream = Teuchos::rcp(&bhs, false);
00239 
00240   // Save the format state of the original std::cout.
00241   Teuchos::oblackholestream oldFormatState;
00242   oldFormatState.copyfmt(std::cout);
00243  
00244   *outStream \
00245   << "===============================================================================\n" \
00246   << "|                                                                             |\n" \
00247   << "|                         Unit Test (AdaptiveSparseGrid)                      |\n" \
00248   << "|                                                                             |\n" \
00249   << "|     1) Integrate product Gaussians in 5 dimensions (Genz integration test). |\n" \
00250   << "|                                                                             |\n" \
00251   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00252   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00253   << "|                                                                             |\n" \
00254   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00255   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00256   << "|                                                                             |\n" \
00257   << "===============================================================================\n"\
00258   << "| TEST 19: Integrate a product of Gaussians in 5D                             |\n"\
00259   << "===============================================================================\n";
00260 
00261 
00262   // internal variables:
00263   int         errorFlag    = 0;
00264   long double TOL          = INTREPID_TOL;
00265   int         dimension    = 5;
00266   int         maxLevel     = 7;
00267   bool        isNormalized = true;                  
00268 
00269   std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_PATTERSON);
00270   std::vector<EIntrepidGrowth>   growth1D(dimension,GROWTH_FULLEXP);
00271  
00272   alpha.resize(dimension,0); beta.resize(dimension,0);
00273   for (int i=0; i<dimension; i++) {
00274     alpha[i] = (long double)std::rand()/(long double)RAND_MAX;
00275     beta[i]  = (long double)std::rand()/(long double)RAND_MAX;
00276   }
00277 
00278   ASGdata<long double,StdVector<long double> > problem_data(
00279         dimension,rule1D,growth1D,maxLevel,isNormalized);  
00280   Teuchos::RCP<std::vector<long double> > integralValue = 
00281     Teuchos::rcp(new std::vector<long double>(1,0.0));
00282   StdVector<long double> sol(integralValue); sol.Set(0.0);
00283   problem_data.init(sol);
00284 
00285   long double eta = adaptSG(sol,problem_data,TOL); 
00286 
00287   long double analyticInt = compExactInt();
00288   long double abstol      = std::sqrt(INTREPID_TOL);
00289   long double absdiff     = fabs(analyticInt-sol[0]);
00290   try { 
00291     *outStream << "Adaptive Sparse Grid exited with global error " 
00292                << std::scientific << std::setprecision(16) << eta << "\n"
00293                << "Approx = " << std::scientific << std::setprecision(16) << sol[0] 
00294                << ",  Exact = " << std::scientific << std::setprecision(16) << analyticInt << "\n"
00295                << "Error = " << std::scientific << std::setprecision(16) << absdiff << "   " 
00296                << "<?" << "   " << abstol << "\n"; 
00297     if (absdiff > abstol) {
00298       errorFlag++;
00299       *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00300     }
00301   }
00302   catch (std::logic_error err) {    
00303     *outStream << err.what() << "\n";
00304     errorFlag = -1;
00305   };  
00306 
00307   if (errorFlag != 0)
00308     std::cout << "End Result: TEST FAILED\n";
00309   else
00310     std::cout << "End Result: TEST PASSED\n";
00311 
00312   // reset format state of std::cout
00313   std::cout.copyfmt(oldFormatState);
00314 
00315   return errorFlag;
00316 }