44template <
typename ordinal_type,
typename value_type>
50 const Teuchos::ParameterList& params_) :
52 pce_basis(pce[0].basis()),
53 pce_sz(pce_basis->size()),
56 verbose(params.get(
"Verbose", false)),
57 rank_threshold(params.get(
"Rank Threshold", 1.0e-12)),
58 orthogonalization_method(params.get(
"Orthogonalization Method",
63template <
typename ordinal_type,
typename value_type>
73 Teuchos::Array< const Stokhos::OrthogPolyApprox<ordinal_type, value_type>* > pce2;
74 for (ordinal_type i=0; i<pce.size(); i++) {
75 if (pce[i].standard_deviation() > 1.0e-15)
76 pce2.push_back(&pce[i]);
81 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
82 const Teuchos::Array< Teuchos::Array<value_type> >& points =
83 quad->getQuadPoints();
84 const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
85 quad->getBasisAtQuadPoints();
86 ordinal_type nqp = weights.size();
91 for (ordinal_type i=0; i<nqp; i++)
92 for (ordinal_type
j=0;
j<pce_sz;
j++)
93 A(i,
j) = basis_values[i][
j];
96 Teuchos::Array<value_type> pce_norms(d, 0.0);
97 for (ordinal_type
j=0;
j<d;
j++) {
98 for (ordinal_type i=0; i<pce_sz; i++)
99 pce_norms[
j] += (*pce2[
j])[i]*(*pce2[
j])[i]*pce_basis->norm_squared(i);
100 pce_norms[
j] = std::sqrt(pce_norms[
j]);
110 Teuchos::Array< Teuchos::Array<value_type> > values(nqp);
111 for (ordinal_type i=0; i<nqp; i++)
112 for (ordinal_type
j=0;
j<d;
j++)
113 F(i,
j) = pce2[
j]->evaluate(points[i], basis_values[i]);
116 sz = buildReducedBasis(max_p, rank_threshold, A, F, weights, terms, num_terms,
120 Teuchos::ParameterList quad_params = params.sublist(
"Reduced Quadrature");
124 if (quad_params.isParameter(
"Reduced Quadrature Method") &&
125 quad_params.get<std::string>(
"Reduced Quadrature Method") ==
"Q2") {
126 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> > terms2;
127 Teuchos::Array<ordinal_type> num_terms2;
128 value_type rank_threshold2 = quad_params.get(
"Q2 Rank Threshold",
132 buildReducedBasis(2*max_p, rank_threshold2, A, F, weights, terms2,
133 num_terms2, Qp2,
Q2);
138 norms.resize(sz, 1.0);
141template <
typename ordinal_type,
typename value_type>
147template <
typename ordinal_type,
typename value_type>
155template <
typename ordinal_type,
typename value_type>
163template <
typename ordinal_type,
typename value_type>
171template <
typename ordinal_type,
typename value_type>
172const Teuchos::Array<value_type>&
179template <
typename ordinal_type,
typename value_type>
187template <
typename ordinal_type,
typename value_type>
188Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
193 return Teuchos::null;
196template <
typename ordinal_type,
typename value_type>
197Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
202 return Teuchos::null;
205template <
typename ordinal_type,
typename value_type>
210 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented!");
213template <
typename ordinal_type,
typename value_type>
216evaluateBases(
const Teuchos::ArrayView<const value_type>& point,
217 Teuchos::Array<value_type>& basis_vals)
const
219 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::logic_error,
"Not implemented!");
222template <
typename ordinal_type,
typename value_type>
225print(std::ostream& os)
const
227 os <<
"Gram-Schmidt basis of order " << p <<
", dimension " << d
228 <<
", and size " << sz <<
". Matrix coefficients:\n";
230 os <<
"Basis vector norms (squared):\n\t";
231 for (ordinal_type i=0; i<sz; i++)
232 os << norms[i] <<
" ";
236template <
typename ordinal_type,
typename value_type>
240 ordinal_type ncol,
bool transpose)
const
243 SDM zbar(Teuchos::View,
const_cast<value_type*
>(in), ncol, ncol, sz);
244 SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
246 z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, Qp, 0.0);
247 TEUCHOS_ASSERT(ret == 0);
250 SDM zbar(Teuchos::View,
const_cast<value_type*
>(in), sz, sz, ncol);
251 SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
253 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Qp, zbar, 0.0);
254 TEUCHOS_ASSERT(ret == 0);
258template <
typename ordinal_type,
typename value_type>
262 ordinal_type ncol,
bool transpose)
const
265 SDM z(Teuchos::View,
const_cast<value_type*
>(in), ncol, ncol, pce_sz);
266 SDM zbar(Teuchos::View, out, ncol, ncol, sz);
268 zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, z, Qp, 0.0);
269 TEUCHOS_ASSERT(ret == 0);
272 SDM z(Teuchos::View,
const_cast<value_type*
>(in), pce_sz, pce_sz, ncol);
273 SDM zbar(Teuchos::View, out, sz, sz, ncol);
275 zbar.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Qp, z, 0.0);
276 TEUCHOS_ASSERT(ret == 0);
280template <
typename ordinal_type,
typename value_type>
281Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
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 dimension() const
Return dimension of basis.
GSReducedPCEBasisBase(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 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 value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
void setup(ordinal_type p, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad)
virtual ordinal_type size() const
Return total size of basis.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual void print(std::ostream &os) const
Print basis to stream os.
ordinal_type order() const
Return order of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
virtual ~GSReducedPCEBasisBase()
Destructor.
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.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
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.
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)