96 const int np = pmax-pmin+1;
98 bool use_pce_quad_points =
false;
99 bool normalize =
false;
100 bool project_integrals =
false;
102 bool sparse_grid =
true;
103#ifndef HAVE_STOKHOS_DAKOTA
106 Teuchos::Array<double> mean(np), mean_st(np), std_dev(np), std_dev_st(np);
107 Teuchos::Array<double> pt(np), pt_st(np);
109 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
111 Teuchos::Array<double> eval_pt(d, 0.56789);
116 for (
int p=pmin; p<=pmax; p++) {
118 std::cout <<
"p = " << p << std::endl;
121 for (
int i=0; i<d; i++)
123 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
128 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double>* > xi(d);
129 for (
int i=0; i<d; i++) {
131 xi[i]->term(i, 1) = 1.0;
135 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
136#ifdef HAVE_STOKHOS_DAKOTA
139 Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(basis, p));
146 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
147 basis->computeTripleProductTensor();
153 s_quad_func<d> s_func(a);
156 for (
int i=0; i<d; i++)
158 quad_exp.nary_op(s_func,s,xip);
159 quad_exp.divide(t1,1.0,s);
163 Teuchos::Array<double> xx(d);
164 for (
int i=0; i<d; i++)
165 xx[i] = xi[i]->evaluate(eval_pt);
166 pt_true = s_func(&(xx[0]));
167 pt_true = 1.0/pt_true;
170 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > st_bases(1);
171 Teuchos::RCP< Stokhos::LanczosProjPCEBasis<int,double> > stp_basis_s;
172 Teuchos::RCP< Stokhos::LanczosPCEBasis<int,double> > st_basis_s;
174 if (project_integrals) {
177 p, Teuchos::rcp(&s,
false), Cijk, normalize,
true));
178 st_bases[0] = stp_basis_s;
183 p, Teuchos::rcp(&s,
false), quad, normalize,
true));
184 st_bases[0] = st_basis_s;
190 p, Teuchos::rcp(&s,
false), quad, use_pce_quad_points,
191 normalize, project_integrals, Cijk));
194 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> >
201 if (project_integrals) {
202 s_st.term(0, 0) = stp_basis_s->getNewCoeffs(0);
203 s_st.term(0, 1) = stp_basis_s->getNewCoeffs(1);
206 s_st.term(0, 0) = st_basis_s->getNewCoeffs(0);
207 s_st.term(0, 1) = st_basis_s->getNewCoeffs(1);
211 s_st.term(0, 0) = s.mean();
212 s_st.term(0, 1) = 1.0;
216 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
217 st_basis->computeTripleProductTensor();
220 Teuchos::RCP<const Stokhos::Quadrature<int,double> > st_quad;
221#ifdef HAVE_STOKHOS_DAKOTA
223 st_quad = Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(st_basis, p));
234 st_quad_exp.divide(t_st, 1.0, s_st);
238 quad_exp.unary_op(st_func, t2, s);
245 mean_st[n] = t2.mean();
246 std_dev[n] = t1.standard_deviation();
247 std_dev_st[n] = t2.standard_deviation();
248 pt[n] = t1.evaluate(eval_pt);
249 pt_st[n] = t2.evaluate(eval_pt);
252 for (
int i=0; i<d; i++)
258 std::cout <<
"Statistical error:" << std::endl;
260 << std::setw(wi) <<
"mean" <<
" "
261 << std::setw(wi) <<
"mean_st" <<
" "
262 << std::setw(wi) <<
"std_dev" <<
" "
263 << std::setw(wi) <<
"std_dev_st" <<
" "
264 << std::setw(wi) <<
"point" <<
" "
265 << std::setw(wi) <<
"point_st" << std::endl;
266 for (
int p=pmin; p<pmax; p++) {
267 std::cout.precision(3);
268 std::cout.setf(std::ios::scientific);
269 std::cout << p <<
" "
270 << std::setw(wi) <<
rel_err(mean[n], mean[np-1]) <<
" "
271 << std::setw(wi) <<
rel_err(mean_st[n], mean[np-1]) <<
" "
272 << std::setw(wi) <<
rel_err(std_dev[n], std_dev[np-1]) <<
" "
273 << std::setw(wi) <<
rel_err(std_dev_st[n], std_dev[np-1])
275 << std::setw(wi) <<
rel_err(pt[n], pt_true) <<
" "
276 << std::setw(wi) <<
rel_err(pt_st[n], pt_true)
282 catch (std::exception& e) {
283 std::cout << e.what() << std::endl;