54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_RefCountPtr.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
59 using namespace Intrepid;
60 std::vector<long double> alpha(1,0);
61 std::vector<long double> beta(1,0);
62 template<
class Scalar>
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;
112 template<
class Scalar,
class UserVector>
118 ASGdata(
int dimension,std::vector<EIntrepidBurkardt> rule1D,
119 std::vector<EIntrepidGrowth> growth1D,
int maxLevel,
121 dimension,rule1D,growth1D,maxLevel,isNormalized) {}
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;
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);
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;
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));
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);
226 int 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 (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);
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.
Scalar error_indicator(UserVector &input)
User defined error indicator function.
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.