57 const Teuchos::ParameterList& params_) :
58 name(
"Product Lanczos Gram-Schmidt PCE Basis"),
62 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > pce_basis = pce[0].basis();
63 pce_sz = pce_basis->size();
66 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(pce_basis);
67 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coord_bases;
68 if (prod_basis != Teuchos::null)
69 coord_bases = prod_basis->getCoordinateBases();
72 bool project =
params.get(
"Project",
true);
73 bool normalize =
params.get(
"Normalize",
true);
74 bool limit_integration_order =
params.get(
"Limit Integration Order",
false);
75 bool use_stieltjes =
params.get(
"Use Old Stieltjes Method",
false);
76 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double > > > coordinate_bases;
77 Teuchos::Array<int> is_invariant(pce.size(),-2);
78 for (ordinal_type i=0; i<pce.size(); i++) {
84 if (prod_basis != Teuchos::null)
86 if (is_invariant[i] >= 0) {
87 coordinate_bases.push_back(coord_bases[is_invariant[i]]);
92 else if (is_invariant[i] != -1) {
94 coordinate_bases.push_back(
97 p, Teuchos::rcp(&(pce[i]),
false), quad,
false,
98 normalize, project, Cijk)));
102 coordinate_bases.push_back(
105 p, Teuchos::rcp(&(pce[i]),
false), Cijk,
106 normalize, limit_integration_order)));
108 coordinate_bases.push_back(
111 p, Teuchos::rcp(&(pce[i]),
false), quad,
112 normalize, limit_integration_order)));
116 d = coordinate_bases.size();
123 params.get(
"Cijk Drop Tolerance", 1.0e-15),
124 params.get(
"Use Old Cijk Algorithm",
false)));
127 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
128 const Teuchos::Array< Teuchos::Array<value_type> >& points =
129 quad->getQuadPoints();
130 const Teuchos::Array< Teuchos::Array<value_type> >& basis_vals =
131 quad->getBasisAtQuadPoints();
132 ordinal_type nqp = weights.size();
134 for (ordinal_type i=0; i<
pce_sz; i++)
135 for (ordinal_type k=0; k<nqp; k++)
136 Psi(i,k) = basis_vals[k][i]*weights[k]/pce_basis->norm_squared(i);
140 Teuchos::Array<value_type> red_basis_vals(
sz);
141 Teuchos::Array<value_type> pce_vals(
d);
144 for (
int k=0; k<nqp; k++) {
145 ordinal_type jdx = 0;
146 for (
int j=0;
j<pce.size();
j++) {
149 if (is_invariant[
j] != -1) {
152 if (is_invariant[
j] >= 0)
153 pce_vals[jdx] = points[k][is_invariant[
j]];
155 pce_vals[jdx] = pce[
j].evaluate(points[k], basis_vals[k]);
156 F(k,jdx) = pce_vals[jdx];
163 for (
int i=0; i<
sz; i++)
164 Phi(i,k) = red_basis_vals[i];
167 bool verbose =
params.get(
"Verbose",
false);
172 A.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, Psi, Phi, 0.0);
173 TEUCHOS_ASSERT(ret == 0);
176 const Teuchos::Array<value_type>& basis_norms = pce_basis->norm_squared();
177 for (ordinal_type
j=0;
j<
sz;
j++) {
178 value_type nrm = 0.0;
179 for (ordinal_type i=0; i<
pce_sz; i++)
180 nrm += A(i,
j)*A(i,
j)*basis_norms[i];
181 nrm = std::sqrt(nrm);
182 for (ordinal_type i=0; i<
pce_sz; i++)
189 value_type rank_threshold =
params.get(
"Rank Threshold", 1.0e-12);
190 std::string orthogonalization_method =
191 params.get(
"Orthogonalization Method",
"Householder");
192 Teuchos::Array<value_type> w(
pce_sz, 1.0);
194 Teuchos::Array<ordinal_type> piv(
sz);
195 for (
int i=0; i<
d+1; i++)
198 sz = SOF::createOrthogonalBasis(
199 orthogonalization_method, rank_threshold, verbose, A, w,
Qp, R, piv);
204 for (ordinal_type i=0; i<nqp; i++)
206 B(i,
j) = basis_vals[i][
j];
210 ret =
Q.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, B,
Qp, 0.0);
211 TEUCHOS_ASSERT(ret == 0);
215 params.sublist(
"Reduced Quadrature"));
216 Teuchos::SerialDenseMatrix<ordinal_type, value_type>
Q2;
217 reduced_quad = quad_factory.createReducedQuadrature(
Q,
Q2, F, weights);
331 ordinal_type ncol,
bool transpose)
const
334 SDM zbar(Teuchos::View,
const_cast<value_type*
>(in), ncol, ncol, sz);
335 SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
337 z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, Qp, 0.0);
338 TEUCHOS_ASSERT(ret == 0);
341 SDM zbar(Teuchos::View,
const_cast<value_type*
>(in), sz, sz, ncol);
342 SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
344 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Qp, zbar, 0.0);
345 TEUCHOS_ASSERT(ret == 0);
353 ordinal_type ncol,
bool transpose)
const
356 SDM z(Teuchos::View,
const_cast<value_type*
>(in), ncol, ncol, pce_sz);
357 SDM zbar(Teuchos::View, out, ncol, ncol, sz);
359 zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, z, Qp, 0.0);
360 TEUCHOS_ASSERT(ret == 0);
363 SDM z(Teuchos::View,
const_cast<value_type*
>(in), pce_sz, pce_sz, ncol);
364 SDM zbar(Teuchos::View, out, sz, sz, ncol);
366 zbar.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Qp, z, 0.0);
367 TEUCHOS_ASSERT(ret == 0);
384 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis =
386 ordinal_type dim = basis->dimension();
387 value_type tol = 1.0e-15;
390 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(basis);
391 if (prod_basis == Teuchos::null)
397 Teuchos::Array<ordinal_type> dependent_dims;
398 tmp_pce.reset(basis);
399 for (ordinal_type i=0; i<dim; i++) {
400 ordinal_type p = prod_basis->getCoordinateBases()[i]->order();
402 for (ordinal_type
j=1;
j<=p;
j++)
403 tmp_pce.term(i,
j) = pce.
term(i,
j);
404 value_type nrm = tmp_pce.two_norm();
405 if (nrm > tol) dependent_dims.push_back(i);
410 if (dependent_dims.size() == 1)
411 return dependent_dims[0];
414 else if (dependent_dims.size() == 0)
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)