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
59using namespace Intrepid;
60std::vector<long double> alpha(1,0);
61std::vector<long double> beta(1,0);
62template<class Scalar>
63class StdVector {
64private:
65 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
66
67public:
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
112template<class Scalar, class UserVector>
113class ASGdata :
114 public Intrepid::AdaptiveSparseGridInterface<Scalar,UserVector> {
115public:
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
147long 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
187long 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
199long double adaptSG(StdVector<long double> & iv,
200 AdaptiveSparseGridInterface<long double,StdVector<long double> > &
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) {
219 eta = AdaptiveSparseGrid<long double,StdVector<long double> >::refine_grid(
220 activeIndex,oldIndex,
221 iv,eta,problem_data);
222 }
223 return eta;
224}
225
226int 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
278 ASGdata<long double,StdVector<long double> > problem_data(
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 (const 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}
Header file for the Intrepid::AdaptiveSparseGrid class.
Intrepid utilities.
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.
Definition test_19.cpp:123
Scalar error_indicator(UserVector &input)
User defined error indicator function.
Definition test_19.cpp:134
bool isNormalized()
Return whether or not cubature weights are normalized.