50 template <
typename ordinal_type,
typename value_type>
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();
67 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coord_bases;
68 if (prod_basis != Teuchos::null)
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);
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();
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++) {
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();
180 nrm += A(i,
j)*A(i,
j)*basis_norms[i];
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);
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);
223 template <
typename ordinal_type,
typename value_type>
229 template <
typename ordinal_type,
typename value_type>
237 template <
typename ordinal_type,
typename value_type>
245 template <
typename ordinal_type,
typename value_type>
253 template <
typename ordinal_type,
typename value_type>
254 const Teuchos::Array<value_type>&
261 template <
typename ordinal_type,
typename value_type>
269 template <
typename ordinal_type,
typename value_type>
270 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
275 return Teuchos::null;
278 template <
typename ordinal_type,
typename value_type>
279 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
284 return Teuchos::null;
287 template <
typename ordinal_type,
typename value_type>
292 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented!");
295 template <
typename ordinal_type,
typename value_type>
299 Teuchos::Array<value_type>& basis_vals)
const 301 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented!");
304 template <
typename ordinal_type,
typename value_type>
309 os <<
"Product Lanczos Gram-Schmidt basis of order " << p <<
", dimension " 311 <<
", and size " << sz <<
". Matrix coefficients:\n";
313 os <<
"Basis vector norms (squared):\n\t";
315 os << norms[i] <<
" ";
319 template <
typename ordinal_type,
typename value_type>
327 template <
typename ordinal_type,
typename value_type>
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);
349 template <
typename ordinal_type,
typename value_type>
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);
371 template <
typename ordinal_type,
typename value_type>
372 Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
379 template <
typename ordinal_type,
typename value_type>
384 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis =
391 if (prod_basis == Teuchos::null)
397 Teuchos::Array<ordinal_type> dependent_dims;
398 tmp_pce.reset(basis);
400 ordinal_type p = prod_basis->getCoordinateBases()[i]->order();
403 tmp_pce.term(i,
j) = pce.
term(i,
j);
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)
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
ordinal_type p
Total order of basis.
virtual const std::string & getName() const
Return string name of basis.
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > reduced_quad
Reduced quadrature object.
ordinal_type pce_sz
Size of original pce basis.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
SDM Qp
Coefficients of transformed basis in original basis.
ordinal_type order() const
Return order of basis.
SDM Q
Values of transformed basis at quadrature points.
virtual ordinal_type size() const
Return total size of basis.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
ordinal_type isInvariant(const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce) const
ordinal_type sz
Total size of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
Teuchos::ParameterList params
Algorithm parameters.
virtual void print(std::ostream &os) const
Print basis to stream os.
Abstract base class for quadrature methods.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
virtual ~ProductLanczosGramSchmidtPCEBasis()
Destructor.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
ordinal_type d
Total dimension of basis.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const =0
Return array of coordinate bases.
virtual void transformFromOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients from original basis to this basis.
ordinal_type dimension() const
Return dimension of basis.
Teuchos::Array< value_type > norms
Norms.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
Teuchos::RCP< Stokhos::CompletePolynomialBasis< ordinal_type, value_type > > tensor_lanczos_basis
Product Lanczos basis.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
ProductLanczosGramSchmidtPCEBasis(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 Teuchos::ParameterList ¶ms=Teuchos::ParameterList())
Constructor.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual void transformToOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients to original basis from this basis.
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Encapsulate various orthogonalization (ie QR) methods.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.