48template <
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;
69 for (ordinal_type i=0; i<pce.size(); i++) {
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;
106 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(
pce_basis);
107 ordinal_type dim = prod_basis->dimension();
108 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
109 for (ordinal_type i=0; i<dim; i++)
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()
128 for (ordinal_type i=0; i<
d; i++) {
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();
144 ordinal_type nqp = weights.size();
148 ordinal_type pce_sz2 = basis2->size();
149 SDM AA(nqp, pce_sz2);
150 Teuchos::Array<value_type> basis_vals(pce_sz2);
151 for (ordinal_type i=0; i<nqp; i++) {
152 basis2->evaluateBases(points[i], basis_vals);
153 for (ordinal_type
j=0;
j<pce_sz2;
j++)
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();
178template <
typename ordinal_type,
typename value_type>
183 value_type threshold,
187 Teuchos::Array<ordinal_type>& num_terms_,
190 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > pce_basis_ = pce[0].basis();
191 ordinal_type pce_sz_ = pce_basis_->size();
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();
199 ordinal_type nqp = weights.size();
203 A_.reshape(nqp, pce_sz_);
204 for (ordinal_type i=0; i<nqp; i++)
205 for (ordinal_type
j=0;
j<pce_sz_;
j++)
206 A_(i,
j) = basis_values[i][
j];
209 Teuchos::Array<value_type> pce_norms(d, 0.0);
210 for (ordinal_type
j=0;
j<d;
j++) {
211 for (ordinal_type i=0; i<pce_sz_; i++)
212 pce_norms[
j] += (pce[
j])[i]*(pce[
j])[i]*pce_basis_->norm_squared(i);
213 pce_norms[
j] = std::sqrt(pce_norms[
j]);
223 Teuchos::Array< Teuchos::Array<value_type> > values(nqp);
224 for (ordinal_type i=0; i<nqp; i++)
225 for (ordinal_type
j=0;
j<d;
j++)
226 F_(i,
j) = pce[
j].evaluate(points[i], basis_values[i]);
231 CPBUtils::compute_terms(max_p, d, max_sz, terms_, num_terms_);
240 for (ordinal_type i=0; i<nqp; i++) {
241 for (ordinal_type
j=0;
j<max_sz;
j++) {
243 for (ordinal_type k=0; k<d; k++)
244 B(i,
j) *= std::pow(F_(i,k), terms_[
j][k]);
249 SDM Bp(pce_sz_, max_sz);
250 const Teuchos::Array<value_type>& basis_norms =
251 pce_basis_->norm_squared();
252 for (ordinal_type i=0; i<pce_sz_; i++) {
253 for (ordinal_type
j=0;
j<max_sz;
j++) {
255 for (ordinal_type k=0; k<nqp; k++)
256 Bp(i,
j) += weights[k]*B(k,
j)*A_(k,i);
257 Bp(i,
j) /= basis_norms[i];
262 for (ordinal_type
j=0;
j<max_sz;
j++) {
263 value_type nrm = 0.0;
264 for (ordinal_type i=0; i<pce_sz_; i++)
265 nrm += Bp(i,
j)*Bp(i,
j)*basis_norms[i];
266 nrm = std::sqrt(nrm);
267 for (ordinal_type i=0; i<pce_sz_; 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++)
280 ordinal_type sz_ = SOF::createOrthogonalBasis(
281 orthogonalization_method, threshold, verbose, Bp, w, Qp_, R, piv);
286template <
typename ordinal_type,
typename value_type>
292template <
typename ordinal_type,
typename value_type>
300template <
typename ordinal_type,
typename value_type>
308template <
typename ordinal_type,
typename value_type>
316template <
typename ordinal_type,
typename value_type>
317const Teuchos::Array<value_type>&
324template <
typename ordinal_type,
typename value_type>
332template <
typename ordinal_type,
typename value_type>
333Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
338 return Teuchos::null;
341template <
typename ordinal_type,
typename value_type>
342Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
347 return Teuchos::null;
350template <
typename ordinal_type,
typename value_type>
355 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented!");
358template <
typename ordinal_type,
typename value_type>
361evaluateBases(
const Teuchos::ArrayView<const value_type>& point,
362 Teuchos::Array<value_type>& basis_vals)
const
364 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented!");
367template <
typename ordinal_type,
typename value_type>
375template <
typename ordinal_type,
typename value_type>
378print(std::ostream& os)
const
380 os <<
"Gram-Schmidt basis of order " << p <<
", dimension " << d
381 <<
", and size " << sz <<
". Matrix coefficients:\n";
383 os <<
"Basis vector norms (squared):\n\t";
384 for (ordinal_type i=0; i<sz; i++)
385 os << norms[i] <<
" ";
389template <
typename ordinal_type,
typename value_type>
393 ordinal_type ncol,
bool transpose)
const
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);
411template <
typename ordinal_type,
typename value_type>
415 ordinal_type ncol,
bool transpose)
const
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);
433template <
typename ordinal_type,
typename value_type>
434Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
int quad2(const Epetra_Map &map, bool verbose)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
SDM Q
Values of transformed basis at quadrature points.
virtual ~MonomialProjGramSchmidtPCEBasis2()
Destructor.
Teuchos::ParameterList params
Algorithm parameters.
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
2-D array of basis terms
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.
virtual const std::string & getName() const
Return string name of basis.
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 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.
Teuchos::Array< ordinal_type > num_terms
Number of terms up to each order.
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.
ordinal_type order() const
Return order of basis.
ordinal_type dimension() const
Return dimension of basis.
Teuchos::Array< value_type > norms
Norms.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
value_type rank_threshold
Rank threshold.
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.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
ordinal_type d
Total dimension of basis.
ordinal_type sz
Total size of basis.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > pce_basis
Original pce basis.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
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.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
SDM Qp
Coefficients of transformed basis in original basis.
A multidimensional index.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
Encapsulate various orthogonalization (ie QR) methods.
Abstract base class for quadrature methods.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
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.
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)