99 Teuchos::GlobalMPISession mpiSession (&argc, &
argv);
104 Teuchos::CommandLineProcessor
CLP;
105 CLP.setDocString(
"This example tests the / operator.\n");
107 CLP.setOption(
"dim", &d,
"Stochastic dimension");
109 CLP.setOption(
"order", &p,
"Polynomial order");
111 CLP.setOption(
"samples", &n,
"Number of samples");
113 CLP.setOption(
"shift", &shift,
"Shift point");
114 double tolerance = 1e-6;
115 CLP.setOption(
"tolerance", &tolerance,
"Tolerance in Iterative Solver");
117 CLP.setOption(
"prec_level", &prec_level,
"Level in Schur Complement Prec 0->Solve A0u0=g0 with division; 1->Form 1x1 Schur Complement");
119 CLP.setOption(
"max_it_div", &max_it_div,
"Maximum # of iterations for division iterative solver");
120 bool equilibrate =
true;
121 CLP.setOption(
"equilibrate",
"noequilibrate", &equilibrate,
122 "Equilibrate the linear system");
125 CLP.setOption(
"solver", &solve_method,
129 CLP.setOption(
"prec", &prec_method,
131 "Preconditioner Method");
133 CLP.setOption(
"schur_option", &schur_option,
137 CLP.setOption(
"prec_option", &prec_option,
143 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
144 for (
int i=0; i<d; i++) {
147 RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
151 RCP<const Stokhos::Quadrature<int,double> > quad =
155 RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
156 basis->computeTripleProductTensor();
159 Teuchos::RCP<Stokhos::QuadOrthogPolyExpansion<int,double> > quad_expn =
162 Teuchos::RCP<Teuchos::ParameterList> alg_params =
163 Teuchos::rcp(
new Teuchos::ParameterList);
165 alg_params->set(
"Division Tolerance", tolerance);
166 alg_params->set(
"prec_iter", prec_level);
167 alg_params->set(
"max_it_div", max_it_div);
169 alg_params->set(
"Division Strategy",
"Dense Direct");
170 else if (solve_method ==
GMRES)
171 alg_params->set(
"Division Strategy",
"GMRES");
172 else if (solve_method ==
CG)
173 alg_params->set(
"Division Strategy",
"CG");
176 if (prec_method ==
None)
177 alg_params->set(
"Prec Strategy",
"None");
178 else if (prec_method ==
Diag)
179 alg_params->set(
"Prec Strategy",
"Diag");
180 else if (prec_method ==
Jacobi)
181 alg_params->set(
"Prec Strategy",
"Jacobi");
182 else if (prec_method ==
GS)
183 alg_params->set(
"Prec Strategy",
"GS");
184 else if (prec_method ==
Schur)
185 alg_params->set(
"Prec Strategy",
"Schur");
187 if (schur_option ==
diag)
188 alg_params->set(
"Schur option",
"diag");
190 alg_params->set(
"Schur option",
"full");
191 if (prec_option ==
linear)
192 alg_params->set(
"Prec option",
"linear");
194 alg_params->set(
"Use Quadrature for Division",
false);
197 alg_params->set(
"Equilibrate", 1);
199 alg_params->set(
"Equilibrate", 0);
203 Teuchos::RCP<Stokhos::QuadOrthogPolyExpansion<int,double> > alg_expn =
205 basis, Cijk, quad, alg_params));
208 pce_type u_quad(quad_expn), v_quad(quad_expn);
209 u_quad.term(0,0) = 0.0;
210 for (
int i=0; i<d; i++) {
211 u_quad.term(i,1) = 1.0;
213 pce_type u_alg(alg_expn), v_alg(alg_expn);
214 u_alg.term(0,0) = 0.0;
215 for (
int i=0; i<d; i++) {
216 u_alg.term(i,1) = 1.0;
220 double scale = std::exp(shift);
221 pce_type b_alg = std::exp(shift + u_alg)/scale;
222 pce_type b_quad = std::exp(shift + u_quad)/scale;
225 v_alg = 1.0 / std::exp(shift + u_alg);
226 v_quad = 1.0 /std::exp(shift + u_quad);
236 double h = 2.0 / (n-1);
237 double err_quad = 0.0;
238 double err_alg = 0.0;
239 for (
int i=0; i<n; i++) {
241 double x = -1.0 + h*i;
243 for (
int j=0;
j<d;
j++)
245 double up = u_quad.evaluate(pt);
246 double vp = 1.0/(shift+up);
247 double vp_quad = v_quad.evaluate(pt);
248 double vp_alg = v_alg.evaluate(pt);
252 double point_err_quad = std::abs(vp-vp_quad);
253 double point_err_alg = std::abs(vp-vp_alg);
254 if (point_err_quad > err_quad) err_quad = point_err_quad;
255 if (point_err_alg > err_alg) err_alg = point_err_alg;
257 std::cout <<
"\tL_infty norm of quadrature error = " << err_quad
259 std::cout <<
"\tL_infty norm of solver error = " << err_alg
264 std::cout <<
"\nExample Passed!" << std::endl;
266 Teuchos::TimeMonitor::summarize(std::cout);
268 catch (std::exception& e) {
269 std::cout << e.what() << std::endl;