44 #include "Teuchos_Assert.hpp" 45 #include "Teuchos_BLAS.hpp" 46 #include "Teuchos_TimeMonitor.hpp" 48 template <
typename ordinal_type,
typename value_type>
54 bool use_pce_quad_points_,
56 bool project_integrals_,
61 pce_weights(quad->getQuadWeights()),
62 basis_values(quad->getBasisAtQuadPoints()),
65 use_pce_quad_points(use_pce_quad_points_),
66 fromStieltjesMat(p+1,pce->size()),
67 project_integrals(project_integrals_),
76 template <
typename ordinal_type,
typename value_type>
82 template <
typename ordinal_type,
typename value_type>
86 Teuchos::Array<value_type>& quad_points,
87 Teuchos::Array<value_type>& quad_weights,
88 Teuchos::Array< Teuchos::Array<value_type> >& quad_values)
const 90 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 91 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::StieltjesPCEBasis -- compute Gauss points");
95 if (use_pce_quad_points) {
96 quad_points = pce_vals;
97 quad_weights = pce_weights;
98 quad_values = phi_vals;
108 if (quad_order > 2*this->p)
109 quad_order = 2*this->p;
116 if (quad_weights.size() < num_points) {
118 quad_weights.resize(num_points);
119 quad_points.resize(num_points);
120 quad_values.resize(num_points);
123 quad_points[i] = quad_points[0];
124 quad_values[i].resize(this->p+1);
125 this->evaluateBases(quad_points[i], quad_values[i]);
130 template <
typename ordinal_type,
typename value_type>
134 Teuchos::Array<value_type>& alpha,
135 Teuchos::Array<value_type>& beta,
136 Teuchos::Array<value_type>& delta,
137 Teuchos::Array<value_type>& gamma)
const 140 Teuchos::Array<value_type> nrm(n);
141 Teuchos::Array< Teuchos::Array<value_type> > vals(nqp);
144 stieltjes(0, n, pce_weights, pce_vals, alpha, beta, nrm, vals);
157 template <
typename ordinal_type,
typename value_type>
163 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
164 quad->getQuadPoints();
166 pce_vals.resize(nqp);
167 phi_vals.resize(nqp);
169 pce_vals[i] = pce->evaluate(quad_points[i], basis_values[i]);
170 phi_vals[i].resize(this->p+1);
173 if (project_integrals)
174 phi_pce_coeffs.resize(basis->size());
179 fromStieltjesMat.putScalar(0.0);
183 fromStieltjesMat(i,
j) +=
184 pce_weights[k]*phi_vals[k][i]*basis_values[k][
j];
185 fromStieltjesMat(i,
j) /= basis->norm_squared(
j);
190 template <
typename ordinal_type,
typename value_type>
195 const Teuchos::Array<value_type>& weights,
196 const Teuchos::Array<value_type>& points,
197 Teuchos::Array<value_type>& a,
198 Teuchos::Array<value_type>& b,
199 Teuchos::Array<value_type>& nrm,
200 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals)
const 202 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 203 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::StieltjesPCEBasis -- Discretized Stieltjes Procedure");
209 if (project_integrals)
210 integrateBasisSquaredProj(0, a, b, weights, points, phi_vals, val1, val2);
212 integrateBasisSquared(0, a, b, weights, points, phi_vals, val1, val2);
219 if (project_integrals)
220 integrateBasisSquaredProj(i, a, b, weights, points, phi_vals, val1, val2);
222 integrateBasisSquared(i, a, b, weights, points, phi_vals, val1, val2);
225 TEUCHOS_TEST_FOR_EXCEPTION(val1 < 0.0, std::logic_error,
226 "Stokhos::StieltjesPCEBasis::stieltjes(): " 227 <<
" Polynomial " << i <<
" out of " << nfinish
228 <<
" has norm " << val1
229 <<
"! Try increasing number of quadrature points");
232 b[i] = nrm[i]/nrm[i-1];
239 template <
typename ordinal_type,
typename value_type>
243 const Teuchos::Array<value_type>& a,
244 const Teuchos::Array<value_type>& b,
245 const Teuchos::Array<value_type>& weights,
246 const Teuchos::Array<value_type>& points,
247 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals,
250 evaluateRecurrence(k, a, b, points, phi_vals);
255 val1 += weights[i]*phi_vals[i][k]*phi_vals[i][k];
256 val2 += weights[i]*phi_vals[i][k]*phi_vals[i][k]*points[i];
260 template <
typename ordinal_type,
typename value_type>
264 const Teuchos::Array<value_type>& a,
265 const Teuchos::Array<value_type>& b,
266 const Teuchos::Array<value_type>& points,
267 Teuchos::Array< Teuchos::Array<value_type> >& values)
const 275 values[i][k] = points[i] - a[k-1];
279 (points[i] - a[k-1])*values[i][k-1] - b[k-1]*values[i][k-2];
282 template <
typename ordinal_type,
typename value_type>
287 const Teuchos::Array<value_type>& a,
288 const Teuchos::Array<value_type>& b,
289 const Teuchos::Array<value_type>& weights,
290 const Teuchos::Array<value_type>& points,
291 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals,
296 const Teuchos::Array<value_type>& norms = basis->norm_squared();
299 evaluateRecurrence(k, a, b, points, phi_vals);
303 c += weights[i]*phi_vals[i][k]*basis_values[i][
j];
305 phi_pce_coeffs[
j] = c;
311 val1 += phi_pce_coeffs[
j]*phi_pce_coeffs[
j]*norms[
j];
316 k_it != Cijk->k_end(); ++k_it) {
319 j_it != Cijk->j_end(k_it); ++j_it) {
322 i_it != Cijk->i_end(j_it); ++i_it) {
325 val2 += phi_pce_coeffs[i]*phi_pce_coeffs[
j]*(*pce)[l]*c;
331 template <
typename ordinal_type,
typename value_type>
336 Teuchos::BLAS<ordinal_type, value_type> blas;
337 blas.GEMV(Teuchos::TRANS, fromStieltjesMat.numRows(),
338 fromStieltjesMat.numCols(), 1.0, fromStieltjesMat.values(),
339 fromStieltjesMat.numRows(), in, 1, 0.0, out, 1);
342 template <
typename ordinal_type,
typename value_type>
343 Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
350 template <
typename ordinal_type,
typename value_type>
356 pce_weights(quad->getQuadWeights()),
357 basis_values(quad->getBasisAtQuadPoints()),
358 pce_vals(sbasis.pce_vals),
360 use_pce_quad_points(sbasis.use_pce_quad_points),
361 fromStieltjesMat(p+1,pce->size()),
362 project_integrals(sbasis.project_integrals),
void integrateBasisSquared(ordinal_type k, const Teuchos::Array< value_type > &a, const Teuchos::Array< value_type > &b, const Teuchos::Array< value_type > &weights, const Teuchos::Array< value_type > &points, Teuchos::Array< Teuchos::Array< value_type > > &phi_vals, value_type &val1, value_type &val2) const
Compute and .
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship: for ...
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
virtual bool computeRecurrenceCoefficients(ordinal_type n, Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &delta, Teuchos::Array< value_type > &gamma) const
Compute recurrence coefficients.
void evaluateRecurrence(ordinal_type k, const Teuchos::Array< value_type > &a, const Teuchos::Array< value_type > &b, const Teuchos::Array< value_type > &points, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Evaluate polynomials via 3-term recurrence.
virtual Teuchos::RCP< OneDOrthogPolyBasis< ordinal_type, value_type > > cloneWithOrder(ordinal_type p) const
Clone this object with the option of building a higher order basis.
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Compute quadrature points, weights, and values of basis polynomials at given set of points points...
StieltjesPCEBasis(ordinal_type p, const Teuchos::RCP< const Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, bool use_pce_quad_points, bool normalize=false, bool project_integrals=false, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk=Teuchos::null)
Constructor.
Bi-directional iterator for traversing a sparse array.
virtual void getQuadPoints(ordinal_type quad_order, Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &weights, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Get Gauss quadrature points, weights, and values of basis at points.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
Abstract base class for quadrature methods.
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
void transformCoeffsFromStieltjes(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
~StieltjesPCEBasis()
Destructor.
virtual void setup()
Setup basis after computing recurrence coefficients.
virtual void setup()
Setup basis after computing recurrence coefficients.
void integrateBasisSquaredProj(ordinal_type k, const Teuchos::Array< value_type > &a, const Teuchos::Array< value_type > &b, const Teuchos::Array< value_type > &weights, const Teuchos::Array< value_type > &points, Teuchos::Array< Teuchos::Array< value_type > > &phi_vals, value_type &val1, value_type &val2) const
Compute and by projecting onto original PCE basis.
void stieltjes(ordinal_type nstart, ordinal_type nfinish, const Teuchos::Array< value_type > &weights, const Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &a, Teuchos::Array< value_type > &b, Teuchos::Array< value_type > &nrm, Teuchos::Array< Teuchos::Array< value_type > > &phi_vals) const
Compute 3-term recurrence using Stieljtes procedure.