209 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(options.d);
210 for (
int i=0; i<options.d; i++)
211 bases[i] = Teuchos::rcp(
new basis_type(p,
true));
212 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
215 results.basis_size = basis->size();
218 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad;
219 if (options.quad_method ==
TENSOR)
222 else if (options.quad_method ==
SPARSE) {
223#ifdef HAVE_STOKHOS_DAKOTA
224 if (options.level == -1)
226 Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(basis));
229 Teuchos::rcp(
new Stokhos::SparseGridQuadrature<int,double>(
230 basis, p, 1e-12, Pecos::SLOW_RESTRICTED_GROWTH));
232 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Sparse grid quadrature only supported when compiled with Dakota!");
236 results.quad_size = quad->size();
239 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
240 basis->computeTripleProductTensor();
243 Teuchos::RCP<Stokhos::QuadOrthogPolyExpansion<int,double> > quad_exp =
249 Teuchos::Array<double> point(options.d, 1.0);
250 Teuchos::Array<double> basis_vals(basis->size());
251 basis->evaluateBases(point, basis_vals);
252 Teuchos::Array<pce_type> x(options.d);
253 for (
int i=0; i<options.d; i++) {
255 x[i].reset(quad_exp);
256 x[i].term(i,1) = 1.0 / basis_vals[i+1];
258 Teuchos::Array<pce_type> x2(options.d2);
259 for (
int i=0; i<options.d2; i++) {
264 pce_type y =
f(x, options.pole, options.shift);
265 results.z =
g(x2, y);
271 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double> > pces(options.d2+1);
272 for (
int i=0; i<options.d2; i++)
273 pces[i] = x2[i].getOrthogPolyApprox();
274 pces[options.d2] = y.getOrthogPolyApprox();
275 Teuchos::ParameterList params;
277 params.set(
"Reduced Basis Method",
"Monomial Proj Gram-Schmidt");
279 params.set(
"Reduced Basis Method",
"Monomial Proj Gram-Schmidt2");
280 else if (options.reduced_basis_method ==
MONOMIAL_GS)
281 params.set(
"Reduced Basis Method",
"Monomial Gram-Schmidt");
282 else if (options.reduced_basis_method ==
LANCZOS)
283 params.set(
"Reduced Basis Method",
"Product Lanczos");
284 else if (options.reduced_basis_method ==
LANCZOS_GS)
285 params.set(
"Reduced Basis Method",
"Product Lanczos Gram-Schmidt");
286 params.set(
"Verbose", options.verbose);
287 params.set(
"Project", options.project);
289 params.set(
"Use Old Stieltjes Method", options.use_stieltjes);
290 params.set(
"Orthogonalization Method",
292 params.set(
"Rank Threshold", options.rank_threshold);
293 Teuchos::ParameterList& red_quad_params =
294 params.sublist(
"Reduced Quadrature");
296 "Reduced Quadrature Method",
301 "Eliminate Dependent Rows", options.eliminate_dependent_rows);
302 red_quad_params.set(
"Write MPS File",
false);
303 red_quad_params.set(
"Reduction Tolerance", options.reduction_tolerance);
304 red_quad_params.set(
"Verbose", options.verbose);
305 red_quad_params.set(
"Objective Value", options.objective_value);
306 red_quad_params.set(
"Q2 Rank Threshold", options.rank_threshold2);
308 "Orthogonalization Method",
310 red_quad_params.set(
"Use Q in LP", options.use_Q);
311 red_quad_params.set(
"Restrict Rank", options.restrict_r);
312 red_quad_params.set(
"Order Restriction", 2*p);
314 Teuchos::RCP< Stokhos::ReducedPCEBasis<int,double> > gs_basis =
316 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gs_quad =
317 gs_basis->getReducedQuadrature();
319 results.reduced_basis_size = gs_basis->size();
320 results.reduced_quad_size = gs_quad->size();
323 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk;
324 Teuchos::RCP< Teuchos::ParameterList > gs_exp_params =
325 Teuchos::rcp(
new Teuchos::ParameterList);
326 if (options.reduced_basis_method ==
LANCZOS)
327 gs_Cijk = gs_basis->computeTripleProductTensor();
329 gs_Cijk = Teuchos::null;
330 gs_exp_params->set(
"Use Quadrature for Times",
true);
332 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<int,double> > gs_quad_exp =
334 gs_basis, gs_Cijk, gs_quad, gs_exp_params));
337 Teuchos::Array<pce_type> x2_gs(options.d2);
338 for (
int i=0; i<options.d2; i++) {
339 x2_gs[i].copyForWrite();
340 x2_gs[i].reset(gs_quad_exp);
341 gs_basis->transformFromOriginalBasis(x2[i].coeff(), x2_gs[i].coeff());
344 gs_basis->transformFromOriginalBasis(y.coeff(), y_gs.coeff());
351 gs_basis->transformToOriginalBasis(z_gs.coeff(), z2.coeff());
363 Teuchos::CommandLineProcessor
CLP;
365 "This example runs a Gram-Schmidt-based dimension reduction example.\n");
369 CLP.setOption(
"d", &options.d,
"Stochastic dimension");
372 CLP.setOption(
"d2", &options.d2,
"Intermediate stochastic dimension");
375 CLP.setOption(
"p_begin", &options.p_begin,
"Starting polynomial order");
378 CLP.setOption(
"p_end", &options.p_end,
"Ending polynomial order");
381 CLP.setOption(
"p_truth", &options.p_truth,
"Truth polynomial order");
384 CLP.setOption(
"pole", &options.pole,
"Pole location");
387 CLP.setOption(
"shift", &options.shift,
"Shift location");
389 options.rank_threshold = 1.0e-120;
390 CLP.setOption(
"rank_threshold", &options.rank_threshold,
"Rank threshold");
392 options.rank_threshold2 = 1.0e-120;
393 CLP.setOption(
"rank_threshold2", &options.rank_threshold2,
"Rank threshold for Q2");
395 options.reduction_tolerance = 1.0e-12;
396 CLP.setOption(
"reduction_tolerance", &options.reduction_tolerance,
"Quadrature reduction tolerance");
398 options.verbose =
false;
399 CLP.setOption(
"verbose",
"quiet", &options.verbose,
"Verbose output");
401 options.quad_method =
TENSOR;
402 CLP.setOption(
"quadrature_method", &options.quad_method,
407 CLP.setOption(
"level", &options.level,
408 "Sparse grid level (set to -1 to use default)");
411 CLP.setOption(
"reduced_basis_method", &options.reduced_basis_method,
415 options.orthogonalization_method =
MGSRO;
416 CLP.setOption(
"orthogonalization_method",
417 &options.orthogonalization_method,
421 "Orthogonalization method");
423 options.quad_reduction_method =
QSQUARED;
424 CLP.setOption(
"reduced_quadrature_method", &options.quad_reduction_method,
428 options.solver_method =
TRSM;
431 "Reduced quadrature solver method");
433 options.quad_orthogonalization_method =
MGSRO;
434 CLP.setOption(
"quad_orthogonalization_method",
435 &options.quad_orthogonalization_method,
439 "Quadrature Orthogonalization method");
441 options.eliminate_dependent_rows =
true;
442 CLP.setOption(
"cpqr",
"no-cpqr", &options.eliminate_dependent_rows,
443 "Eliminate dependent rows in quadrature constraint matrix");
445 options.use_Q =
true;
446 CLP.setOption(
"use-Q",
"no-use-Q", &options.use_Q,
"Use Q in LP");
448 options.restrict_r =
false;
449 CLP.setOption(
"restrict-rank",
"no-restrict-rank", &options.restrict_r,
450 "Restrict rank in LP");
452 options.objective_value = 0.0;
453 CLP.setOption(
"objective_value", &options.objective_value,
454 "Value for LP objective function");
456 options.project =
true;
457 CLP.setOption(
"project",
"no-project", &options.project,
458 "Use Projected Lanczos Method");
460 options.use_stieltjes =
false;
461 CLP.setOption(
"stieltjes",
"no-stieltjes", &options.use_stieltjes,
462 "Use Old Stieltjes Method");
466 std::cout <<
"Summary of command line options:" << std::endl
467 <<
"\tquadrature_method = "
470 <<
"\tlevel = " << options.level
472 <<
"\treduced_basis_method = "
475 <<
"\torthogonalization_method = "
478 <<
"\tquadrature_reduction_method = "
481 <<
"\tsolver_method = "
483 <<
"\tquad_orthogonalization_method = "
486 <<
"\tcpqr = " << options.eliminate_dependent_rows
488 <<
"\tuse-Q = " << options.use_Q << std::endl
489 <<
"\trestrict-rank = " << options.restrict_r << std::endl
490 <<
"\tobjective_value = " << options.objective_value << std::endl
491 <<
"\tproject = " << options.project << std::endl
492 <<
"\tstieljtes = " << options.use_stieltjes << std::endl
493 <<
"\tp_begin = " << options.p_begin << std::endl
494 <<
"\tp_end = " << options.p_end << std::endl
495 <<
"\tp_truth = " << options.p_truth << std::endl
496 <<
"\td = " << options.d << std::endl
497 <<
"\td2 = " << options.d2 << std::endl
498 <<
"\tpole = " << options.pole << std::endl
499 <<
"\tshift = " << options.shift << std::endl
500 <<
"\trank_threshold = " << options.rank_threshold << std::endl
501 <<
"\trank_threshold2 = " << options.rank_threshold2 << std::endl
502 <<
"\treduction_tolerance = " << options.reduction_tolerance
504 <<
"\tverbose = " << options.verbose << std::endl
505 << std::endl << std::endl;
507 std::stringstream ss;
508 typedef Teuchos::TabularOutputter TO;
510 out.setFieldTypePrecision(TO::DOUBLE, 5);
511 out.pushFieldSpec(
"\"Order\"", TO::INT);
512 out.pushFieldSpec(
"\"Basis Size\"", TO::INT);
513 out.pushFieldSpec(
"\"Quad. Size\"", TO::INT);
514 out.pushFieldSpec(
"\"Red. Basis Size\"", TO::INT);
515 out.pushFieldSpec(
"\"Red. Quad. Size\"", TO::INT);
516 out.pushFieldSpec(
"\"Coeff. Error\"", TO::DOUBLE);
517 out.pushFieldSpec(
"\"Red. Coeff. Error\"", TO::DOUBLE);
518 out.pushFieldSpec(
"\"Disc. Orthog. Error\"", TO::DOUBLE);
522 compute_pces(
false, options.p_truth, options, results_truth);
524 int n = options.p_end - options.p_begin + 1;
525 Teuchos::Array<MyResults> results(n);
526 for (
int i=0; i<n; ++i) {
527 int p = options.p_begin + i;
529 results[i].mean_error =
530 std::abs(results_truth.z.mean()-results[i].z.mean()) /
531 std::abs(results_truth.z.mean());
532 results[i].std_dev_error =
533 std::abs(results_truth.z.standard_deviation()-results[i].z.standard_deviation()) / std::abs(results_truth.z.standard_deviation());
534 results[i].coeff_error =
coeff_error(results_truth.z,
537 results[i].reduced_mean_error =
538 std::abs(results_truth.z.mean()-results[i].z_red.mean()) /
539 std::abs(results_truth.z.mean());
540 results[i].reduced_std_dev_error =
541 std::abs(results_truth.z.standard_deviation()-results[i].z_red.standard_deviation()) / std::abs(results_truth.z.standard_deviation());
542 results[i].reduced_coeff_error =
coeff_error(results_truth.z,
546 out.outputField(results[i].basis_size);
547 out.outputField(results[i].quad_size);
548 out.outputField(results[i].reduced_basis_size);
549 out.outputField(results[i].reduced_quad_size);
551 out.outputField(results[i].reduced_coeff_error);
555 std::cout << std::endl << ss.str() << std::endl;
558 catch (std::exception& e) {
559 std::cout << e.what() << std::endl;