50#include "Teuchos_CommandLineProcessor.hpp"
59 Teuchos::CommandLineProcessor
CLP;
61 "This example runs a Gram-Schmidt-based dimension reduction example.\n");
64 CLP.setOption(
"d", &d,
"Stochastic dimension");
67 CLP.setOption(
"p", &p,
"Polynomial order");
72 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
73 for (
int i=0; i<d; i++)
74 bases[i] = Teuchos::rcp(
new basis_type(p,
true));
75 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
78 std::cout <<
"original basis size = " << basis->size() << std::endl;
83 for (
int i=0; i<d; i++) {
88 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
94 std::cout <<
"original quadrature size = " << quad->size() << std::endl;
97 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
98 basis->computeTripleProductTensor();
106 quad_exp.
times(w,v,u);
109 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double> > pces(2);
112 Teuchos::ParameterList params;
113 params.set(
"Reduced Basis Method",
"Monomial Proj Gram-Schmidt");
115 Teuchos::RCP< Stokhos::ReducedPCEBasis<int,double> > gs_basis =
117 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gs_quad =
118 gs_basis->getReducedQuadrature();
121 gs_basis->transformFromOriginalBasis(u.coeff(), u_gs.coeff());
122 gs_basis->transformFromOriginalBasis(v.coeff(), v_gs.coeff());
124 std::cout <<
"reduced basis size = " << gs_basis->size() << std::endl;
125 std::cout <<
"reduced quadrature size = " << gs_quad->size() << std::endl;
128 gs_basis->transformToOriginalBasis(u_gs.coeff(), u2.coeff());
129 gs_basis->transformToOriginalBasis(v_gs.coeff(), v2.
coeff());
132 for (
int i=0; i<basis->size(); i++) {
133 double eu = std::abs(u[i]-u2[i]);
134 double ev = std::abs(v[i]-v2[i]);
135 if (eu > err_u) err_u = eu;
136 if (ev > err_v) err_v = ev;
138 std::cout <<
"error in u transformation = " << err_u << std::endl;
139 std::cout <<
"error in v transformation = " << err_v << std::endl;
142 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk =
146 Teuchos::RCP< Teuchos::ParameterList > gs_exp_params =
147 Teuchos::rcp(
new Teuchos::ParameterList);
148 gs_exp_params->set(
"Use Quadrature for Times",
true);
155 gs_quad_exp.
times(w_gs, u_gs, v_gs);
158 gs_basis->transformToOriginalBasis(w_gs.
coeff(), w2.
coeff());
160 std::cout.precision(12);
161 std::cout <<
"w = " << std::endl << w;
162 std::cout <<
"w2 = " << std::endl << w2;
163 std::cout <<
"w_gs = " << std::endl << w_gs;
166 for (
int i=0; i<basis->size(); i++) {
167 double ew = std::abs(w[i]-w2[i]);
168 if (ew > err_w) err_w = ew;
171 std::cout.setf(std::ios::scientific);
172 std::cout <<
"w.mean() = " << w.mean() << std::endl
173 <<
"w2.mean() = " << w2.
mean() << std::endl
175 << std::abs(w.mean()-w2.
mean()) << std::endl
176 <<
"w.std_dev() = " << w.standard_deviation() << std::endl
178 <<
"std_dev error = "
181 <<
"w coeff error = " << err_w << std::endl;
184 catch (std::exception& e) {
185 std::cout << e.what() << std::endl;
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
value_type mean() const
Compute mean of expansion.
value_type standard_deviation() const
Compute standard deviation of expansion.
pointer coeff()
Return coefficient array.
Orthogonal polynomial expansions based on numerical quadrature.
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void times(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)
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual Teuchos::RCP< Stokhos::ReducedPCEBasis< ordinal_type, value_type > > createReducedBasis(ordinal_type p, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk) const
Get reduced quadrature object.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
int main(int argc, char **argv)
Stokhos::LegendreBasis< int, double > basis_type