57 static const int N = Num;
64 for (
int i=0; i<Num; i++)
66 return 1.0 / (prod -
a);
82 mutable Teuchos::Array<double>
vec;
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();
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);
245 mean_st[n] = t2.
mean();
246 std_dev[n] = t1.standard_deviation();
248 pt[n] = t1.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;
int main(int argc, char **argv)
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
double operator()(double x[Num]) const
pce_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
Teuchos::Array< double > vec
void divide(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
value_type standard_deviation() const
Compute standard deviation of expansion.
const Stokhos::OrthogPolyBasis< int, double > & basis
void nary_op(const FuncT &func, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > **a)
double rel_err(double a, double b)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
const Stokhos::OrthogPolyApprox< int, double > & pce
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
value_type mean() const
Compute mean of expansion.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Orthogonal polynomial expansions based on numerical quadrature.
double operator()(const double &a, const double &b) const
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
void unary_op(const FuncT &func, OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Nonlinear unary function.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
Stokhos::LegendreBasis< int, double > basis_type