42#include "Teuchos_Assert.hpp"
43#include "Teuchos_TimeMonitor.hpp"
45template <
typename ordinal_type,
typename value_type,
typename func_type>
49 const Teuchos::RCP<const func_type>& func_,
51 bool use_pce_quad_points_,
53 RecurrenceBasis<ordinal_type, value_type>(
"Stieltjes PCE", p, normalize),
56 pce_weights(quad->getQuadWeights()),
57 basis_values(quad->getBasisAtQuadPoints()),
60 use_pce_quad_points(use_pce_quad_points_)
66template <
typename ordinal_type,
typename value_type,
typename func_type>
72template <
typename ordinal_type,
typename value_type,
typename func_type>
76 Teuchos::Array<value_type>& quad_points,
77 Teuchos::Array<value_type>& quad_weights,
78 Teuchos::Array< Teuchos::Array<value_type> >& quad_values)
const
80#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
81 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::StieltjesBasis -- compute Gauss points");
85 if (use_pce_quad_points) {
86 quad_points = func_vals;
87 quad_weights = pce_weights;
88 quad_values = phi_vals;
93 ordinal_type num_points =
94 static_cast<ordinal_type
>(std::ceil((quad_order+1)/2.0));
97 if (quad_order > 2*this->p)
98 quad_order = 2*this->p;
105 if (quad_weights.size() < num_points) {
106 ordinal_type old_size = quad_weights.size();
107 quad_weights.resize(num_points);
108 quad_points.resize(num_points);
109 quad_values.resize(num_points);
110 for (ordinal_type i=old_size; i<num_points; i++) {
111 quad_weights[i] = value_type(0);
112 quad_points[i] = quad_points[0];
113 quad_values[i].resize(this->p+1);
114 this->evaluateBases(quad_points[i], quad_values[i]);
119template <
typename ordinal_type,
typename value_type,
typename func_type>
123 Teuchos::Array<value_type>& alpha,
124 Teuchos::Array<value_type>& beta,
125 Teuchos::Array<value_type>& delta,
126 Teuchos::Array<value_type>& gamma)
const
128 ordinal_type nqp = phi_vals.size();
129 Teuchos::Array<value_type> nrm(n);
130 Teuchos::Array< Teuchos::Array<value_type> > vals(nqp);
131 for (ordinal_type i=0; i<nqp; i++)
133 stieltjes(0, n, pce_weights, func_vals, alpha, beta, nrm, vals);
134 for (ordinal_type i=0; i<n; i++) {
135 delta[i] = value_type(1.0);
136 gamma[i] = value_type(1.0);
146template <
typename ordinal_type,
typename value_type,
typename func_type>
152 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
153 quad->getQuadPoints();
154 ordinal_type nqp = pce_weights.size();
155 func_vals.resize(nqp);
156 phi_vals.resize(nqp);
157 for (ordinal_type i=0; i<nqp; i++) {
158 func_vals[i] = (*func)(quad_points[i]);
159 phi_vals[i].resize(this->p+1);
165template <
typename ordinal_type,
typename value_type,
typename func_type>
169 ordinal_type nfinish,
170 const Teuchos::Array<value_type>& weights,
171 const Teuchos::Array<value_type>& points,
172 Teuchos::Array<value_type>& a,
173 Teuchos::Array<value_type>& b,
174 Teuchos::Array<value_type>& nrm,
175 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals)
const
177#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
178 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos::StieltjesBasis -- Discretized Stieltjes Procedure");
181 value_type val1, val2;
182 ordinal_type start = nstart;
184 integrateBasisSquared(0, a, b, weights, points, phi_vals, val1, val2);
187 b[0] = value_type(1);
190 for (ordinal_type i=start; i<nfinish; i++) {
191 integrateBasisSquared(i, a, b, weights, points, phi_vals, val1, val2);
194 TEUCHOS_TEST_FOR_EXCEPTION(val1 < 0.0, std::logic_error,
195 "Stokhos::StieltjesBasis::stieltjes(): "
196 <<
" Polynomial " << i <<
" out of " << nfinish
197 <<
" has norm " << val1
198 <<
"! Try increasing number of quadrature points");
201 b[i] = nrm[i]/nrm[i-1];
205template <
typename ordinal_type,
typename value_type,
typename func_type>
209 const Teuchos::Array<value_type>& a,
210 const Teuchos::Array<value_type>& b,
211 const Teuchos::Array<value_type>& weights,
212 const Teuchos::Array<value_type>& points,
213 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals,
214 value_type& val1, value_type& val2)
const
216 evaluateRecurrence(k, a, b, points, phi_vals);
217 ordinal_type nqp = weights.size();
218 val1 = value_type(0);
219 val2 = value_type(0);
220 for (ordinal_type i=0; i<nqp; i++) {
221 val1 += weights[i]*phi_vals[i][k]*phi_vals[i][k];
222 val2 += weights[i]*phi_vals[i][k]*phi_vals[i][k]*points[i];
226template <
typename ordinal_type,
typename value_type,
typename func_type>
230 const Teuchos::Array<value_type>& a,
231 const Teuchos::Array<value_type>& b,
232 const Teuchos::Array<value_type>& points,
233 Teuchos::Array< Teuchos::Array<value_type> >& values)
const
235 ordinal_type np = points.size();
237 for (ordinal_type i=0; i<np; i++)
240 for (ordinal_type i=0; i<np; i++)
241 values[i][k] = points[i] - a[k-1];
243 for (ordinal_type i=0; i<np; i++)
245 (points[i] - a[k-1])*values[i][k-1] - b[k-1]*values[i][k-2];
248template <
typename ordinal_type,
typename value_type,
typename func_type>
249Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
256template <
typename ordinal_type,
typename value_type,
typename func_type>
262 pce_weights(quad->getQuadWeights()),
263 basis_values(quad->getBasisAtQuadPoints()),
266 use_pce_quad_points(basis.use_pce_quad_points)
Abstract base class for quadrature methods.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship:
virtual void setup()
Setup basis after computing recurrence coefficients.
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.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a functional map...
~StieltjesBasis()
Destructor.
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.
StieltjesBasis(ordinal_type p, const Teuchos::RCP< const func_type > &func, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, bool use_pce_quad_points, bool normalize=false)
Constructor.
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.
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.
virtual void setup()
Setup basis after computing recurrence coefficients.
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.
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 .
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.