Intrepid
test_19.cpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Intrepid Package
5 // Copyright (2007) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38 // Denis Ridzal (dridzal@sandia.gov), or
39 // Kara Peterson (kjpeter@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 
52 //#include "Intrepid_CubatureLineSorted.hpp"
53 #include "Intrepid_Utils.hpp"
54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_RefCountPtr.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
58 
59 using namespace Intrepid;
60 std::vector<long double> alpha(1,0);
61 std::vector<long double> beta(1,0);
62 template<class Scalar>
63 class StdVector {
64 private:
65  Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
66 
67 public:
68 
69  StdVector( const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec )
70  : std_vec_(std_vec) {}
71 
72  Teuchos::RefCountPtr<StdVector<Scalar> > Create() const {
73  return Teuchos::rcp( new StdVector<Scalar>(
74  Teuchos::rcp(new std::vector<Scalar>(std_vec_->size(),0))));
75  }
76 
77  void Update( StdVector<Scalar> & s ) {
78  int dimension = (int)(std_vec_->size());
79  for (int i=0; i<dimension; i++)
80  (*std_vec_)[i] += s[i];
81  }
82 
83  void Update( Scalar alpha, StdVector<Scalar> & s ) {
84  int dimension = (int)(std_vec_->size());
85  for (int i=0; i<dimension; i++)
86  (*std_vec_)[i] += alpha*s[i];
87  }
88 
89  Scalar operator[](int i) {
90  return (*std_vec_)[i];
91  }
92 
93  void clear() {
94  std_vec_->clear();
95  }
96 
97  void resize(int n, Scalar p) {
98  std_vec_->resize(n,p);
99  }
100 
101  int size() {
102  return (int)std_vec_->size();
103  }
104 
105  void Set( Scalar alpha ) {
106  int dimension = (int)(std_vec_->size());
107  for (int i=0; i<dimension; i++)
108  (*std_vec_)[i] = alpha;
109  }
110 };
111 
112 template<class Scalar, class UserVector>
113 class ASGdata :
114  public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> {
115 public:
116  ~ASGdata() {}
117 
118  ASGdata(int dimension,std::vector<EIntrepidBurkardt> rule1D,
119  std::vector<EIntrepidGrowth> growth1D, int maxLevel,
120  bool isNormalized) : AdaptiveSparseGridInterface<Scalar,UserVector>(
121  dimension,rule1D,growth1D,maxLevel,isNormalized) {}
122 
123  void eval_integrand(UserVector & output, std::vector<Scalar> & input) {
124  int dimension = (int)alpha.size();
125  Scalar total = 0.0;
126  Scalar point = 0.0;
127  for (int i=0; i<dimension; i++) {
128  point = 0.5*input[i]+0.5;
129  total += powl(alpha[i]*(point-beta[i]),(long double)2.0);
130  }
131  output.clear(); output.resize(1,std::exp(-total));
132  }
133 
134  Scalar error_indicator(UserVector & input) {
135  int dimension = (int)input.size();
136  Scalar norm2 = 0.0;
137  for (int i=0; i<dimension; i++)
138  norm2 += input[i]*input[i];
139 
142  norm2 = std::sqrt(norm2)/ID;
143  return norm2;
144  }
145 };
146 
147 long double nCDF(long double z) {
148  long double p = 0.0, expntl = 0.0;
149  long double p0 = 220.2068679123761;
150  long double p1 = 221.2135961699311;
151  long double p2 = 112.0792914978709;
152  long double p3 = 33.91286607838300;
153  long double p4 = 6.373962203531650;
154  long double p5 = 0.7003830644436881;
155  long double p6 = 0.03526249659989109;
156  long double q0 = 440.4137358247522;
157  long double q1 = 793.8265125199484;
158  long double q2 = 637.3336333788311;
159  long double q3 = 296.5642487796737;
160  long double q4 = 86.78073220294608;
161  long double q5 = 16.06417757920695;
162  long double q6 = 1.755667163182642;
163  long double q7 = 0.08838834764831844;
164  long double rootpi = std::sqrt(M_PI);
165  long double zabs = std::abs(z);
166 
167  if (12.0 < zabs) {
168  p = 0.0;
169  }
170  else {
171  expntl = exp(-zabs*zabs/2.0);
172  if (zabs < 7.0) {
173  p = expntl*
174  ((((((p6*zabs+p5)*zabs+p4)*zabs+p3)*zabs+p2)*zabs+p1)*zabs+p0)/
175  (((((((q7*zabs+q6)*zabs+q5)*zabs+q4)*zabs+q3)*zabs+q2)*zabs+q1)*zabs+q0);
176  }
177  else {
178  p = expntl/(zabs+1.0/(zabs+2.0/(zabs+3.0/(zabs+4.0/(zabs+0.65)))))/rootpi;
179  }
180  }
181  if(0.0 < z){
182  p = 1.0-p;
183  }
184  return p;
185 }
186 
187 long double compExactInt(void) {
188  long double val = 1.0;
189  int dimension = alpha.size();
190  long double s2 = std::sqrt(2.0);
191  long double sp = std::sqrt(M_PI);
192  for (int i=0; i<dimension; i++) {
193  long double s2a = s2*alpha[i];
194  val *= (sp/alpha[i])*(nCDF((1.0-beta[i])*s2a)-nCDF(-beta[i]*s2a));
195  }
196  return val;
197 }
198 
199 long double adaptSG(StdVector<long double> & iv,
201  problem_data, long double TOL) {
202 
203  // Construct a Container for the adapted rule
204  int dimension = problem_data.getDimension();
205  std::vector<int> index(dimension,1);
206 
207  // Initialize global error indicator
208  long double eta = 1.0;
209 
210  // Initialize the Active index set
211  std::multimap<long double,std::vector<int> > activeIndex;
212  activeIndex.insert(std::pair<long double,std::vector<int> >(eta,index));
213 
214  // Initialize the old index set
215  std::set<std::vector<int> > oldIndex;
216 
217  // Perform Adaptation
218  while (eta > TOL) {
220  activeIndex,oldIndex,
221  iv,eta,problem_data);
222  }
223  return eta;
224 }
225 
226 int main(int argc, char *argv[]) {
227 
228  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
229 
230  // This little trick lets us print to std::cout only if
231  // a (dummy) command-line argument is provided.
232  int iprint = argc - 1;
233  Teuchos::RCP<std::ostream> outStream;
234  Teuchos::oblackholestream bhs; // outputs nothing
235  if (iprint > 0)
236  outStream = Teuchos::rcp(&std::cout, false);
237  else
238  outStream = Teuchos::rcp(&bhs, false);
239 
240  // Save the format state of the original std::cout.
241  Teuchos::oblackholestream oldFormatState;
242  oldFormatState.copyfmt(std::cout);
243 
244  *outStream \
245  << "===============================================================================\n" \
246  << "| |\n" \
247  << "| Unit Test (AdaptiveSparseGrid) |\n" \
248  << "| |\n" \
249  << "| 1) Integrate product Gaussians in 5 dimensions (Genz integration test). |\n" \
250  << "| |\n" \
251  << "| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
252  << "| Denis Ridzal (dridzal@sandia.gov). |\n" \
253  << "| |\n" \
254  << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
255  << "| Trilinos website: http://trilinos.sandia.gov |\n" \
256  << "| |\n" \
257  << "===============================================================================\n"\
258  << "| TEST 19: Integrate a product of Gaussians in 5D |\n"\
259  << "===============================================================================\n";
260 
261 
262  // internal variables:
263  int errorFlag = 0;
264  long double TOL = INTREPID_TOL;
265  int dimension = 5;
266  int maxLevel = 7;
267  bool isNormalized = true;
268 
269  std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_PATTERSON);
270  std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
271 
272  alpha.resize(dimension,0); beta.resize(dimension,0);
273  for (int i=0; i<dimension; i++) {
274  alpha[i] = (long double)std::rand()/(long double)RAND_MAX;
275  beta[i] = (long double)std::rand()/(long double)RAND_MAX;
276  }
277 
279  dimension,rule1D,growth1D,maxLevel,isNormalized);
280  Teuchos::RCP<std::vector<long double> > integralValue =
281  Teuchos::rcp(new std::vector<long double>(1,0.0));
282  StdVector<long double> sol(integralValue); sol.Set(0.0);
283  problem_data.init(sol);
284 
285  long double eta = adaptSG(sol,problem_data,TOL);
286 
287  long double analyticInt = compExactInt();
288  long double abstol = std::sqrt(INTREPID_TOL);
289  long double absdiff = std::abs(analyticInt-sol[0]);
290  try {
291  *outStream << "Adaptive Sparse Grid exited with global error "
292  << std::scientific << std::setprecision(16) << eta << "\n"
293  << "Approx = " << std::scientific << std::setprecision(16) << sol[0]
294  << ", Exact = " << std::scientific << std::setprecision(16) << analyticInt << "\n"
295  << "Error = " << std::scientific << std::setprecision(16) << absdiff << " "
296  << "<?" << " " << abstol << "\n";
297  if (absdiff > abstol) {
298  errorFlag++;
299  *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
300  }
301  }
302  catch (std::logic_error err) {
303  *outStream << err.what() << "\n";
304  errorFlag = -1;
305  };
306 
307  if (errorFlag != 0)
308  std::cout << "End Result: TEST FAILED\n";
309  else
310  std::cout << "End Result: TEST PASSED\n";
311 
312  // reset format state of std::cout
313  std::cout.copyfmt(oldFormatState);
314 
315  return errorFlag;
316 }
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...
Scalar getInitialDiff()
Return initial error indicator.
Header file for the Intrepid::AdaptiveSparseGrid class.
Intrepid utilities.
Scalar error_indicator(UserVector &input)
User defined error indicator function.
Definition: test_19.cpp:134
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.
Definition: test_19.cpp:123