54#include "Teuchos_oblackholestream.hpp"
55#include "Teuchos_RCP.hpp"
56#include "Teuchos_RefCountPtr.hpp"
57#include "Teuchos_GlobalMPISession.hpp"
59using namespace Intrepid;
60std::vector<long double> alpha(1,0);
61std::vector<long double> beta(1,0);
65 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
69 StdVector(
const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec )
70 : std_vec_(std_vec) {}
72 Teuchos::RefCountPtr<StdVector<Scalar> > Create()
const {
74 Teuchos::rcp(
new std::vector<Scalar>(std_vec_->size(),0))));
78 int dimension = (int)(std_vec_->size());
79 for (
int i=0; i<dimension; i++)
80 (*std_vec_)[i] += s[i];
84 int dimension = (int)(std_vec_->size());
85 for (
int i=0; i<dimension; i++)
86 (*std_vec_)[i] += alpha*s[i];
89 Scalar operator[](
int i) {
90 return (*std_vec_)[i];
97 void resize(
int n, Scalar p) {
98 std_vec_->resize(n,p);
102 return (
int)std_vec_->size();
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;
112template<
class Scalar,
class UserVector>
118 ASGdata(
int dimension,std::vector<EIntrepidBurkardt> rule1D,
119 std::vector<EIntrepidGrowth> growth1D,
int maxLevel,
124 int dimension = (int)alpha.size();
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);
131 output.clear(); output.resize(1,std::exp(-total));
135 int dimension = (int)input.size();
137 for (
int i=0; i<dimension; i++)
138 norm2 += input[i]*input[i];
142 norm2 = std::sqrt(norm2)/ID;
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);
171 expntl = exp(-zabs*zabs/2.0);
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);
178 p = expntl/(zabs+1.0/(zabs+2.0/(zabs+3.0/(zabs+4.0/(zabs+0.65)))))/rootpi;
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));
201 problem_data,
long double TOL) {
204 int dimension = problem_data.getDimension();
205 std::vector<int> index(dimension,1);
208 long double eta = 1.0;
211 std::multimap<long double,std::vector<int> > activeIndex;
212 activeIndex.insert(std::pair<
long double,std::vector<int> >(eta,index));
215 std::set<std::vector<int> > oldIndex;
220 activeIndex,oldIndex,
221 iv,eta,problem_data);
226int main(
int argc,
char *argv[]) {
228 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
232 int iprint = argc - 1;
233 Teuchos::RCP<std::ostream> outStream;
234 Teuchos::oblackholestream bhs;
236 outStream = Teuchos::rcp(&std::cout,
false);
238 outStream = Teuchos::rcp(&bhs,
false);
241 Teuchos::oblackholestream oldFormatState;
242 oldFormatState.copyfmt(std::cout);
245 <<
"===============================================================================\n" \
247 <<
"| Unit Test (AdaptiveSparseGrid) |\n" \
249 <<
"| 1) Integrate product Gaussians in 5 dimensions (Genz integration test). |\n" \
251 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
252 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
254 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
255 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
257 <<
"===============================================================================\n"\
258 <<
"| TEST 19: Integrate a product of Gaussians in 5D |\n"\
259 <<
"===============================================================================\n";
264 long double TOL = INTREPID_TOL;
267 bool isNormalized =
true;
269 std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_PATTERSON);
270 std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
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;
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));
283 problem_data.init(sol);
285 long double eta = adaptSG(sol,problem_data,TOL);
287 long double analyticInt = compExactInt();
288 long double abstol = std::sqrt(INTREPID_TOL);
289 long double absdiff = std::abs(analyticInt-sol[0]);
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) {
299 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
302 catch (
const std::logic_error & err) {
303 *outStream << err.what() <<
"\n";
308 std::cout <<
"End Result: TEST FAILED\n";
310 std::cout <<
"End Result: TEST PASSED\n";
313 std::cout.copyfmt(oldFormatState);
Header file for the Intrepid::AdaptiveSparseGrid class.
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.
Scalar error_indicator(UserVector &input)
User defined error indicator function.
Scalar getInitialDiff()
Return initial error indicator.
bool isNormalized()
Return whether or not cubature weights are normalized.
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...