55 const Teuchos::ParameterList& params_) :
56 name(
"Product Lanczos PCE Basis"),
59 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > pce_basis = pce[0].basis();
60 ordinal_type pce_sz = pce_basis->size();
63 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(pce_basis);
64 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coord_bases;
65 if (prod_basis != Teuchos::null)
66 coord_bases = prod_basis->getCoordinateBases();
69 bool project =
params.get(
"Project",
true);
70 bool normalize =
params.get(
"Normalize",
true);
71 bool limit_integration_order =
params.get(
"Limit Integration Order",
false);
72 bool use_stieltjes =
params.get(
"Use Old Stieltjes Method",
false);
73 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double > > > coordinate_bases;
74 Teuchos::Array<int> is_invariant(pce.size(),-2);
75 for (ordinal_type i=0; i<pce.size(); i++) {
81 if (prod_basis != Teuchos::null)
83 if (is_invariant[i] >= 0) {
84 coordinate_bases.push_back(coord_bases[is_invariant[i]]);
89 else if (is_invariant[i] != -1) {
91 coordinate_bases.push_back(
94 p, Teuchos::rcp(&(pce[i]),
false), quad,
false,
95 normalize, project, Cijk)));
99 coordinate_bases.push_back(
102 p, Teuchos::rcp(&(pce[i]),
false), Cijk,
103 normalize, limit_integration_order)));
105 coordinate_bases.push_back(
108 p, Teuchos::rcp(&(pce[i]),
false), quad,
109 normalize, limit_integration_order)));
113 ordinal_type d = coordinate_bases.size();
120 params.get(
"Cijk Drop Tolerance", 1.0e-15),
121 params.get(
"Use Old Cijk Algorithm",
false)));
124 Teuchos::ParameterList sg_params;
125 sg_params.sublist(
"Basis").set< Teuchos::RCP< const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > >(
"Stochastic Galerkin Basis",
tensor_lanczos_basis);
126 sg_params.sublist(
"Quadrature") =
params.sublist(
"Reduced Quadrature");
131 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
132 const Teuchos::Array< Teuchos::Array<value_type> >& points =
133 quad->getQuadPoints();
134 const Teuchos::Array< Teuchos::Array<value_type> >& basis_vals =
135 quad->getBasisAtQuadPoints();
136 ordinal_type nqp = weights.size();
137 SDM Psi(pce_sz, nqp);
138 for (ordinal_type i=0; i<pce_sz; i++)
139 for (ordinal_type k=0; k<nqp; k++)
140 Psi(i,k) = basis_vals[k][i]*weights[k]/pce_basis->norm_squared(i);
144 Teuchos::Array<value_type> red_basis_vals(sz);
145 Teuchos::Array<value_type> pce_vals(d);
147 for (
int k=0; k<nqp; k++) {
148 ordinal_type jdx = 0;
149 for (
int j=0;
j<pce.size();
j++) {
152 if (is_invariant[
j] != -1) {
155 if (is_invariant[
j] >= 0)
156 pce_vals[jdx] = points[k][is_invariant[
j]];
158 pce_vals[jdx] = pce[
j].evaluate(points[k], basis_vals[k]);
165 for (
int i=0; i<sz; i++)
166 Phi(i,k) = red_basis_vals[i];
169 bool verbose =
params.get(
"Verbose",
false);
174 A.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, Psi,
Phi, 0.0);
175 TEUCHOS_ASSERT(ret == 0);
182 Teuchos::Array<value_type> sigma;
184 value_type rank_threshold =
params.get(
"Rank Threshold", 1.0e-12);
185 ordinal_type rank =
svd_threshold(rank_threshold,
A, sigma, U, Vt);
186 Ainv.shape(sz, pce_sz);
187 TEUCHOS_ASSERT(rank == Vt.numRows());
188 for (ordinal_type i=0; i<Vt.numRows(); i++)
189 for (ordinal_type
j=0;
j<Vt.numCols();
j++)
190 Vt(i,
j) = Vt(i,
j) / sigma[i];
191 ret =
Ainv.multiply(Teuchos::TRANS, Teuchos::TRANS, 1.0, Vt, U, 0.0);
192 TEUCHOS_ASSERT(ret == 0);
196 std::cout <<
"rank = " << rank << std::endl;
198 std::cout <<
"diag(S) = [";
199 for (ordinal_type i=0; i<rank; i++)
200 std::cout << sigma[i] <<
" ";
201 std::cout <<
"]" << std::endl;
204 SDM SVt(rank, Vt.numCols());
205 for (ordinal_type i=0; i<Vt.numRows(); i++)
206 for (ordinal_type
j=0;
j<Vt.numCols();
j++)
207 SVt(i,
j) = Vt(i,
j) * sigma[i] * sigma[i];
209 SDM err_A(pce_sz,sz);
211 ret = err_A.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -1.0, U, SVt,
213 TEUCHOS_ASSERT(ret == 0);
214 std::cout <<
"||A - U*S*V^T||_infty = " << err_A.normInf() << std::endl;
220 for (ordinal_type i=0; i<sz; i++)
222 ret = err.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0,
Ainv,
A,
224 TEUCHOS_ASSERT(ret == 0);
225 std::cout <<
"||Ainv*A - I||_infty = " << err.normInf() << std::endl;
229 SDM err2(pce_sz,pce_sz);
231 for (ordinal_type i=0; i<pce_sz; i++)
233 ret = err2.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0,
A,
Ainv, -1.0);
234 TEUCHOS_ASSERT(ret == 0);
235 std::cout <<
"||A*Ainv - I||_infty = " << err2.normInf() << std::endl;
379 ordinal_type ncol,
bool transpose)
const
381 ordinal_type pce_sz = A.numRows();
382 ordinal_type sz = A.numCols();
384 SDM zbar(Teuchos::View,
const_cast<value_type*
>(in), ncol, ncol, sz);
385 SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
387 z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, A, 0.0);
388 TEUCHOS_ASSERT(ret == 0);
391 SDM zbar(Teuchos::View,
const_cast<value_type*
>(in), sz, sz, ncol);
392 SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
394 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A, zbar, 0.0);
395 TEUCHOS_ASSERT(ret == 0);
403 ordinal_type ncol,
bool transpose)
const
405 ordinal_type pce_sz = A.numRows();
406 ordinal_type sz = A.numCols();
408 SDM z(Teuchos::View,
const_cast<value_type*
>(in), ncol, ncol, pce_sz);
409 SDM zbar(Teuchos::View, out, ncol, ncol, sz);
411 zbar.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, z, Ainv, 0.0);
412 TEUCHOS_ASSERT(ret == 0);
415 SDM z(Teuchos::View,
const_cast<value_type*
>(in), pce_sz, pce_sz, ncol);
416 SDM zbar(Teuchos::View, out, sz, sz, ncol);
418 zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Ainv, z, 0.0);
419 TEUCHOS_ASSERT(ret == 0);
436 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis =
438 ordinal_type dim = basis->dimension();
439 value_type tol = 1.0e-15;
442 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(basis);
443 if (prod_basis == Teuchos::null)
449 Teuchos::Array<ordinal_type> dependent_dims;
450 tmp_pce.reset(basis);
451 for (ordinal_type i=0; i<dim; i++) {
452 ordinal_type p = prod_basis->getCoordinateBases()[i]->order();
454 for (ordinal_type
j=1;
j<=p;
j++)
455 tmp_pce.term(i,
j) = pce.
term(i,
j);
456 value_type nrm = tmp_pce.two_norm();
457 if (nrm > tol) dependent_dims.push_back(i);
462 if (dependent_dims.size() == 1)
463 return dependent_dims[0];
466 else if (dependent_dims.size() == 0)