48 template <
typename ordinal_type,
typename value_type>
54 const Teuchos::ParameterList& params_) :
55 name(
"Monomial Proj Gram Schmidt PCE Basis"),
57 pce_basis(pce[0].basis()),
58 pce_sz(pce_basis->size()),
61 verbose(params.get(
"Verbose", false)),
62 rank_threshold(params.get(
"Rank Threshold", 1.0e-12)),
63 orthogonalization_method(params.get(
"Orthogonalization Method",
68 Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> > non_const_pce;
70 if (pce[i].standard_deviation() > 1.0e-15)
71 non_const_pce.push_back(pce[i]);
73 d = non_const_pce.size();
79 Q.reshape(A.numRows(),
sz);
81 Q.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A,
Qp, 0.0);
82 TEUCHOS_ASSERT(ret == 0);
87 Teuchos::ParameterList quad_params =
params.sublist(
"Reduced Quadrature");
91 if (quad_params.isParameter(
"Reduced Quadrature Method") &&
92 quad_params.get<std::string>(
"Reduced Quadrature Method") ==
"Q2") {
93 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> > terms2;
94 Teuchos::Array<ordinal_type> num_terms2;
95 value_type rank_threshold2 = quad_params.get(
"Q2 Rank Threshold",
100 Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> > pce2(non_const_pce);
101 Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
quad2 = quad;
102 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis2 =
pce_basis;
108 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
110 bases[i] = prod_basis->getCoordinateBases()[i]->cloneWithOrder(2*max_p);
114 std::cout <<
"built new basis of dimension " << basis2->dimension()
115 <<
" and order " << basis2->order()
116 <<
" with total size " << basis2->size() << std::endl;
124 std::cout <<
"built new quadrature with total size " <<
quad2->size()
129 pce2[i].reset(basis2);
136 buildQ(2*max_p, rank_threshold2, pce2,
quad2, terms2, num_terms2,
141 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
142 const Teuchos::Array< Teuchos::Array<value_type> >& points =
143 quad->getQuadPoints();
149 SDM AA(nqp, pce_sz2);
150 Teuchos::Array<value_type> basis_vals(pce_sz2);
152 basis2->evaluateBases(points[i], basis_vals);
154 AA(i,
j) = basis_vals[
j];
156 Q2.reshape(nqp, sz2);
157 ret =
Q2.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, AA, Qp2, 0.0);
158 TEUCHOS_ASSERT(ret == 0);
170 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
178 template <
typename ordinal_type,
typename value_type>
187 Teuchos::Array<ordinal_type>& num_terms_,
190 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > pce_basis_ = pce[0].basis();
194 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
195 const Teuchos::Array< Teuchos::Array<value_type> >& points =
196 quad->getQuadPoints();
197 const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
198 quad->getBasisAtQuadPoints();
203 A_.reshape(nqp, pce_sz_);
206 A_(i,
j) = basis_values[i][
j];
209 Teuchos::Array<value_type> pce_norms(d, 0.0);
212 pce_norms[
j] += (pce[
j])[i]*(pce[
j])[i]*pce_basis_->norm_squared(i);
223 Teuchos::Array< Teuchos::Array<value_type> > values(nqp);
226 F_(i,
j) = pce[
j].evaluate(points[i], basis_values[i]);
231 CPBUtils::compute_terms(max_p, d, max_sz, terms_, num_terms_);
249 SDM Bp(pce_sz_, max_sz);
250 const Teuchos::Array<value_type>& basis_norms =
251 pce_basis_->norm_squared();
256 Bp(i,
j) += weights[k]*B(k,
j)*A_(k,i);
257 Bp(i,
j) /= basis_norms[i];
265 nrm += Bp(i,
j)*Bp(i,
j)*basis_norms[i];
274 Teuchos::Array<value_type> w(pce_sz_, 1.0);
276 Teuchos::Array<ordinal_type> piv(max_sz);
277 for (
int i=0; i<d+1; i++)
281 orthogonalization_method, threshold, verbose, Bp, w, Qp_, R, piv);
286 template <
typename ordinal_type,
typename value_type>
292 template <
typename ordinal_type,
typename value_type>
300 template <
typename ordinal_type,
typename value_type>
308 template <
typename ordinal_type,
typename value_type>
316 template <
typename ordinal_type,
typename value_type>
317 const Teuchos::Array<value_type>&
324 template <
typename ordinal_type,
typename value_type>
332 template <
typename ordinal_type,
typename value_type>
333 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
338 return Teuchos::null;
341 template <
typename ordinal_type,
typename value_type>
342 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
347 return Teuchos::null;
350 template <
typename ordinal_type,
typename value_type>
355 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented!");
358 template <
typename ordinal_type,
typename value_type>
362 Teuchos::Array<value_type>& basis_vals)
const 364 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented!");
367 template <
typename ordinal_type,
typename value_type>
375 template <
typename ordinal_type,
typename value_type>
380 os <<
"Gram-Schmidt basis of order " << p <<
", dimension " << d
381 <<
", and size " << sz <<
". Matrix coefficients:\n";
382 os << Qp << std::endl;
383 os <<
"Basis vector norms (squared):\n\t";
385 os << norms[i] <<
" ";
389 template <
typename ordinal_type,
typename value_type>
396 SDM zbar(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, sz);
397 SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
399 z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, Qp, 0.0);
400 TEUCHOS_ASSERT(ret == 0);
403 SDM zbar(Teuchos::View, const_cast<value_type*>(in), sz, sz, ncol);
404 SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
406 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Qp, zbar, 0.0);
407 TEUCHOS_ASSERT(ret == 0);
411 template <
typename ordinal_type,
typename value_type>
418 SDM z(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, pce_sz);
419 SDM zbar(Teuchos::View, out, ncol, ncol, sz);
421 zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, z, Qp, 0.0);
422 TEUCHOS_ASSERT(ret == 0);
425 SDM z(Teuchos::View, const_cast<value_type*>(in), pce_sz, pce_sz, ncol);
426 SDM zbar(Teuchos::View, out, sz, sz, ncol);
428 zbar.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Qp, z, 0.0);
429 TEUCHOS_ASSERT(ret == 0);
433 template <
typename ordinal_type,
typename value_type>
434 Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
virtual void print(std::ostream &os) const
Print basis to stream os.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > reduced_quad
Reduced quadrature object.
MonomialProjGramSchmidtPCEBasis2(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::ParameterList ¶ms=Teuchos::ParameterList())
Constructor.
ordinal_type d
Total dimension of basis.
KOKKOS_INLINE_FUNCTION PCE< Storage > pow(const PCE< Storage > &a, const PCE< Storage > &b)
ordinal_type order() const
Return order of basis.
ordinal_type buildQ(ordinal_type max_p, value_type threshold, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > &terms_, Teuchos::Array< ordinal_type > &num_terms_, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Qp_, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &A_, Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F_)
Build the reduced basis, parameterized by total order max_p.
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.
int quad2(const Epetra_Map &map, bool verbose)
SDM Q
Values of transformed basis at quadrature points.
value_type rank_threshold
Rank threshold.
Abstract base class for quadrature methods.
virtual ordinal_type size() const
Return total size of basis.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
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.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
virtual const std::string & getName() const
Return string name of basis.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
ordinal_type sz
Total size of basis.
Teuchos::ParameterList params
Algorithm parameters.
virtual ~MonomialProjGramSchmidtPCEBasis2()
Destructor.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
SDM Qp
Coefficients of transformed basis in original basis.
ordinal_type dimension() const
Return dimension of basis.
Teuchos::Array< value_type > norms
Norms.
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
2-D array of basis terms
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > pce_basis
Original pce basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
Teuchos::Array< ordinal_type > num_terms
Number of terms up to each order.
virtual Teuchos::RCP< const Stokhos::UserDefinedQuadrature< ordinal_type, value_type > > createReducedQuadrature(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q2, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights) const
Get reduced quadrature object.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
Encapsulate various orthogonalization (ie QR) methods.