53 #include "Teuchos_CommandLineProcessor.hpp" 55 template <
typename ordinal_type,
typename scalar_type>
57 f_func(
const Teuchos::SerialDenseMatrix<ordinal_type,scalar_type>& A,
58 const Teuchos::SerialDenseVector<ordinal_type,scalar_type>&
x,
60 Teuchos::SerialDenseVector<ordinal_type,scalar_type>&
f)
63 Teuchos::SerialDenseVector<ordinal_type,scalar_type>
y(n), s(n);
67 f.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A,
y, 0.0);
70 template <
typename ordinal_type,
typename scalar_type>
72 g_func(
const Teuchos::SerialDenseVector<ordinal_type,scalar_type>&
f)
80 typedef Teuchos::SerialDenseMatrix<int,pce_type>
SDM;
81 typedef Teuchos::SerialDenseVector<int,pce_type>
SDV;
82 typedef Teuchos::SerialDenseMatrix<int,double>
SDM0;
83 typedef Teuchos::SerialDenseVector<int,double>
SDV0;
101 Teuchos::CommandLineProcessor
CLP;
103 "This example runs a Stieltjes-based dimension reduction example.\n");
106 CLP.setOption(
"n", &n,
"Number of random variables");
109 CLP.setOption(
"m", &m,
"Number of intermediate variables");
112 CLP.setOption(
"pmin", &pmin,
"Starting expansion order");
115 CLP.setOption(
"pmax", &pmax,
"Final expansion order");
118 CLP.setOption(
"mt_method", &mt_method,
120 "Measure transformation method");
122 bool normalize =
false;
123 CLP.setOption(
"normalize",
"unnormalize", &normalize,
124 "Normalize PC basis");
126 bool project_integrals =
false;
127 CLP.setOption(
"project",
"no_project", &project_integrals,
128 "Project integrals");
130 #ifdef HAVE_STOKHOS_DAKOTA 131 bool sparse_grid =
true;
132 CLP.setOption(
"sparse",
"tensor", &sparse_grid,
133 "Use Sparse or Tensor grid");
135 bool sparse_grid =
false;
138 double rank_threshold = 1.0e-12;
139 CLP.setOption(
"rank_threshold", &rank_threshold,
"Rank threshold");
141 double reduction_tolerance = 1.0e-12;
142 CLP.setOption(
"reduction_tolerance", &reduction_tolerance,
"Quadrature reduction tolerance");
144 bool verbose =
false;
145 CLP.setOption(
"verbose",
"quient", &verbose,
150 std::cout <<
"Summary of command line options:" << std::endl
151 <<
"\tm = " << m << std::endl
152 <<
"\tn = " << n << std::endl
154 <<
"\tnormalize = " << normalize << std::endl
155 <<
"\tproject = " << project_integrals << std::endl
156 <<
"\tsparse = " << sparse_grid << std::endl
157 <<
"\trank_threshold = " << rank_threshold << std::endl
158 <<
"\treduction_tolerance = " << reduction_tolerance << std::endl;
160 int np = pmax-pmin+1;
162 Teuchos::Array<double> mean(np), mean_st(np), std_dev(np), std_dev_st(np);
163 Teuchos::Array<double> pt(np), pt_st(np);
165 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(n);
167 Teuchos::Array<double> eval_pt(n, 0.56789);
168 double pt_true = 0.0;
170 SDM0 B0(n,n), Q0, R0;
172 Teuchos::Array<double> w(n, 1.0);
174 SDM0 A0(Teuchos::View, B0, m, n);
178 for (
int p=pmin; p<=pmax; p++) {
180 std::cout <<
"p = " << p;
183 for (
int i=0; i<n; i++)
184 bases[i] = Teuchos::rcp(
new basis_type(p, normalize));
185 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
187 std::cout <<
", basis sz = " << basis->size();
190 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
191 #ifdef HAVE_STOKHOS_DAKOTA 194 Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(basis, p));
199 std::cout <<
", quad sz = " << quad->size();
202 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
203 basis->computeTripleProductTensor();
206 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<int,double> > quad_exp =
211 for (
int i=0; i<n; i++) {
213 x[i].reset(quad_exp);
217 x[i].term(i, 1) = 1.0;
222 for (
int i=0; i<m; i++)
223 for (
int j=0;
j<n;
j++)
237 for (
int i=0; i<n; i++)
238 x0[i] =
x[i].evaluate(eval_pt);
243 Teuchos::ParameterList params;
244 params.set(
"Verbose", verbose);
246 params.set(
"Reduced Basis Method",
"Monomial Proj Gram-Schmidt");
248 params.set(
"Reduced Basis Method",
"Product Lanczos");
250 params.set(
"Reduced Basis Method",
"Product Lanczos");
251 params.set(
"Use Old Stieltjes Method",
true);
253 params.set(
"Project", project_integrals);
254 params.set(
"Normalize", normalize);
255 params.set(
"Rank Threshold", rank_threshold);
256 Teuchos::ParameterList& red_quad_params =
257 params.sublist(
"Reduced Quadrature");
258 red_quad_params.set(
"Reduction Tolerance", reduction_tolerance);
259 red_quad_params.set(
"Verbose", verbose);
260 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double> > pces(m);
261 for (
int i=0; i<m; i++)
262 pces[i] =
f[i].getOrthogPolyApprox();
264 Teuchos::RCP< Stokhos::ReducedPCEBasis<int,double> > gs_basis =
266 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gs_quad =
267 gs_basis->getReducedQuadrature();
268 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk =
270 Teuchos::RCP< Teuchos::ParameterList > gs_exp_params =
271 Teuchos::rcp(
new Teuchos::ParameterList);
272 gs_exp_params->set(
"Use Quadrature for Times",
true);
273 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<int,double> > st_quad_exp =
279 std::cout <<
", red. basis sz = " <<gs_basis->size();
280 std::cout <<
", red. quad sz = " << gs_quad->size();
283 for (
int i=0; i<m; i++) {
284 f_st(i).copyForWrite();
285 f_st(i).reset(st_quad_exp);
286 gs_basis->transformFromOriginalBasis(
f(i).coeff(), f_st(i).coeff());
294 gs_basis->transformToOriginalBasis(g_st.coeff(), g2.coeff());
300 mean[idx] =
g.mean();
301 mean_st[idx] = g2.mean();
302 std_dev[idx] =
g.standard_deviation();
303 std_dev_st[idx] = g2.standard_deviation();
304 pt[idx] =
g.evaluate(eval_pt);
305 pt_st[idx] = g2.evaluate(eval_pt);
308 std::cout << std::endl;
313 std::cout <<
"Statistical error:" << std::endl;
315 << std::setw(wi) <<
"mean" <<
" " 316 << std::setw(wi) <<
"mean_st" <<
" " 317 << std::setw(wi) <<
"std_dev" <<
" " 318 << std::setw(wi) <<
"std_dev_st" <<
" " 319 << std::setw(wi) <<
"point" <<
" " 320 << std::setw(wi) <<
"point_st" << std::endl;
321 for (
int p=pmin; p<pmax; p++) {
322 std::cout.precision(3);
323 std::cout.setf(std::ios::scientific);
324 std::cout << p <<
" " 325 << std::setw(wi) <<
rel_err(mean[idx], mean[np-1]) <<
" " 326 << std::setw(wi) <<
rel_err(mean_st[idx], mean[np-1]) <<
" " 327 << std::setw(wi) <<
rel_err(std_dev[idx], std_dev[np-1]) <<
" " 328 << std::setw(wi) <<
rel_err(std_dev_st[idx], std_dev[np-1])
330 << std::setw(wi) <<
rel_err(pt[idx], pt_true) <<
" " 331 << std::setw(wi) <<
rel_err(pt_st[idx], pt_true)
337 catch (std::exception& e) {
338 std::cout << e.what() << std::endl;
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
int main(int argc, char **argv)
void QR_MGS(ordinal_type k, const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::Array< scalar_type > &w, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &Q, Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &R)
Compute thin QR using modified Gram-Schmidt.
Teuchos::SerialDenseVector< int, double > SDV0
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
Teuchos::SerialDenseMatrix< int, pce_type > SDM
Teuchos::SerialDenseVector< int, pce_type > SDV
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Sacado::PCE::OrthogPoly< double, Stokhos::StandardStorage< int, double > > pce_type
void f_func(const Teuchos::SerialDenseMatrix< ordinal_type, scalar_type > &A, const Teuchos::SerialDenseVector< ordinal_type, scalar_type > &x, double shift, Teuchos::SerialDenseVector< ordinal_type, scalar_type > &f)
double rel_err(double a, double b)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
const MT_METHOD mt_method_values[]
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Teuchos::SerialDenseMatrix< int, double > SDM0
Legendre polynomial basis.
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.
ScalarType f(const Teuchos::Array< ScalarType > &x, double a, double b)
scalar_type g_func(const Teuchos::SerialDenseVector< ordinal_type, scalar_type > &f)
const char * mt_method_names[]
Orthogonal polynomial expansions based on numerical quadrature.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
Stokhos::LegendreBasis< int, double > basis_type