44#include "Teuchos_Assert.hpp"
45#include "Teuchos_BLAS.hpp"
46#include "Teuchos_TimeMonitor.hpp"
48template <
typename ordinal_type,
typename value_type>
54 bool use_pce_quad_points_,
56 bool project_integrals_,
58 RecurrenceBasis<ordinal_type, value_type>(
"Stieltjes PCE", p, normalize),
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_),
76template <
typename ordinal_type,
typename value_type>
82template <
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;
103 ordinal_type num_points =
104 static_cast<ordinal_type
>(std::ceil((quad_order+1)/2.0));
108 if (quad_order > 2*this->p)
109 quad_order = 2*this->p;
116 if (quad_weights.size() < num_points) {
117 ordinal_type old_size = quad_weights.size();
118 quad_weights.resize(num_points);
119 quad_points.resize(num_points);
120 quad_values.resize(num_points);
121 for (ordinal_type i=old_size; i<num_points; i++) {
122 quad_weights[i] = value_type(0);
123 quad_points[i] = quad_points[0];
124 quad_values[i].resize(this->p+1);
125 this->evaluateBases(quad_points[i], quad_values[i]);
130template <
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
139 ordinal_type nqp = phi_vals.size();
140 Teuchos::Array<value_type> nrm(n);
141 Teuchos::Array< Teuchos::Array<value_type> > vals(nqp);
142 for (ordinal_type i=0; i<nqp; i++)
144 stieltjes(0, n, pce_weights, pce_vals, alpha, beta, nrm, vals);
145 for (ordinal_type i=0; i<n; i++) {
146 delta[i] = value_type(1.0);
147 gamma[i] = value_type(1.0);
157template <
typename ordinal_type,
typename value_type>
163 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
164 quad->getQuadPoints();
165 ordinal_type nqp = pce_weights.size();
166 pce_vals.resize(nqp);
167 phi_vals.resize(nqp);
168 for (ordinal_type i=0; i<nqp; i++) {
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());
178 ordinal_type sz = pce->size();
179 fromStieltjesMat.putScalar(0.0);
180 for (ordinal_type i=0; i<=this->p; i++) {
181 for (ordinal_type
j=0;
j<sz;
j++) {
182 for (ordinal_type k=0; k<nqp; k++)
183 fromStieltjesMat(i,
j) +=
184 pce_weights[k]*phi_vals[k][i]*basis_values[k][
j];
185 fromStieltjesMat(i,
j) /= basis->norm_squared(
j);
190template <
typename ordinal_type,
typename value_type>
194 ordinal_type nfinish,
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");
206 value_type val1, val2;
207 ordinal_type start = nstart;
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);
215 b[0] = value_type(1);
218 for (ordinal_type i=start; i<nfinish; i++) {
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];
239template <
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,
248 value_type& val1, value_type& val2)
const
250 evaluateRecurrence(k, a, b, points, phi_vals);
251 ordinal_type nqp = weights.size();
252 val1 = value_type(0);
253 val2 = value_type(0);
254 for (ordinal_type i=0; i<nqp; i++) {
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];
260template <
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
269 ordinal_type np = points.size();
271 for (ordinal_type i=0; i<np; i++)
274 for (ordinal_type i=0; i<np; i++)
275 values[i][k] = points[i] - a[k-1];
277 for (ordinal_type i=0; i<np; i++)
279 (points[i] - a[k-1])*values[i][k-1] - b[k-1]*values[i][k-2];
282template <
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,
292 value_type& val1, value_type& val2)
const
294 ordinal_type nqp = weights.size();
295 ordinal_type npc = basis->size();
296 const Teuchos::Array<value_type>& norms = basis->norm_squared();
299 evaluateRecurrence(k, a, b, points, phi_vals);
300 for (ordinal_type
j=0;
j<npc;
j++) {
301 value_type c = value_type(0);
302 for (ordinal_type i=0; i<nqp; i++)
303 c += weights[i]*phi_vals[i][k]*basis_values[i][
j];
305 phi_pce_coeffs[
j] = c;
309 val1 = value_type(0);
310 for (ordinal_type
j=0;
j<npc;
j++)
311 val1 += phi_pce_coeffs[
j]*phi_pce_coeffs[
j]*norms[
j];
314 val2 = value_type(0);
316 k_it != Cijk->k_end(); ++k_it) {
317 ordinal_type l = index(k_it);
319 j_it != Cijk->j_end(k_it); ++j_it) {
320 ordinal_type
j = index(j_it);
322 i_it != Cijk->i_end(j_it); ++i_it) {
323 ordinal_type i = index(i_it);
324 value_type c = value(i_it);
325 val2 += phi_pce_coeffs[i]*phi_pce_coeffs[
j]*(*pce)[l]*c;
331template <
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);
342template <
typename ordinal_type,
typename value_type>
343Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
350template <
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),
Class to store coefficients of a projection onto an orthogonal polynomial basis.
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.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
j_sparse_array::const_iterator kji_iterator
Iterator for looping over i entries given k and j.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
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.
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.
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 .
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 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.
virtual void setup()
Setup basis after computing recurrence coefficients.
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.
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.
~StieltjesPCEBasis()
Destructor.
void transformCoeffsFromStieltjes(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
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.
Bi-directional iterator for traversing a sparse array.