13 template <
typename ordinal_type,
typename value_type>
19 normalize(normalize_),
21 quad_zero_tol(1.0e-14),
22 #ifdef HAVE_STOKHOS_DAKOTA
23 sparse_grid_growth_rule(webbur::level_to_order_linear_nn),
25 sparse_grid_growth_rule(NULL),
35 template <
typename ordinal_type,
typename value_type>
40 normalize(basis.normalize),
42 quad_zero_tol(basis.quad_zero_tol),
43 sparse_grid_growth_rule(basis.sparse_grid_growth_rule),
52 template <
typename ordinal_type,
typename value_type>
58 computeRecurrenceCoefficients(p+1, alpha, beta, delta, gamma);
59 if (normalize && !is_normalized) {
60 normalizeRecurrenceCoefficients(alpha, beta, delta, gamma);
66 norms[0] = beta[0]/(gamma[0]*gamma[0]);
68 norms[k] = (beta[k]/gamma[k])*(delta[k-1]/delta[k])*norms[k-1];
72 template <
typename ordinal_type,
typename value_type>
78 template <
typename ordinal_type,
typename value_type>
86 template <
typename ordinal_type,
typename value_type>
94 template <
typename ordinal_type,
typename value_type>
102 template <
typename ordinal_type,
typename value_type>
110 template <
typename ordinal_type,
typename value_type>
121 getQuadPoints(3*p, points, weights, values);
130 weights[l]*(values[l][i])*(values[l][
j])*(values[l][k]);
132 (*Cijk)(i,
j,k) = triple_product;
140 template <
typename ordinal_type,
typename value_type>
152 getQuadPoints(3*p, points, weights, values);
161 weights[l]*(values[l][i])*(values[l][
j])*(values[l][k]);
163 if (
std::abs(triple_product/norms[i]) > sparse_tol)
164 Cijk->add_term(i,
j,k,triple_product);
168 Cijk->fillComplete();
173 template <
typename ordinal_type,
typename value_type>
181 getQuadPoints(2*p, points, weights, values);
187 evaluateBasesAndDerivatives(points[i], values[i], derivs[i]);
194 for (
int qp=0; qp<nqp; qp++)
195 b += weights[qp]*derivs[qp][i]*values[qp][
j];
203 template <
typename ordinal_type,
typename value_type>
215 basis_pts[1] = (delta[0]*x-alpha[0])*basis_pts[0]/gamma[1];
217 basis_pts[i] = ((delta[i-1]*x-alpha[i-1])*basis_pts[i-1] -
218 beta[i-1]*basis_pts[i-2])/gamma[i];
221 template <
typename ordinal_type,
typename value_type>
228 evaluateBases(x, vals);
231 derivs[1] = delta[0]/(gamma[0]*gamma[1]);
233 derivs[i] = (delta[i-1]*vals[i-1] + (delta[i-1]*x-alpha[i-1])*derivs[i-1] -
234 beta[i-1]*derivs[i-2])/gamma[i];
237 template <
typename ordinal_type,
typename value_type>
252 value_type v1 = (x*delta[0]-alpha[0])*v0/gamma[1];
258 v2 = ((delta[i-1]*x-alpha[i-1])*v1 - beta[i-1]*v0)/gamma[i];
266 template <
typename ordinal_type,
typename value_type>
271 os << name <<
" basis of order " << p <<
"." << std::endl;
273 os <<
"Alpha recurrence coefficients:\n\t";
275 os << alpha[i] <<
" ";
278 os <<
"Beta recurrence coefficients:\n\t";
280 os << beta[i] <<
" ";
283 os <<
"Delta recurrence coefficients:\n\t";
285 os << delta[i] <<
" ";
288 os <<
"Gamma recurrence coefficients:\n\t";
290 os << gamma[i] <<
" ";
293 os <<
"Basis polynomial norms (squared):\n\t";
295 os << norms[i] <<
" ";
299 template <
typename ordinal_type,
typename value_type>
307 template <
typename ordinal_type,
typename value_type>
328 if(num_points > p+1){
329 bool is_normalized = computeRecurrenceCoefficients(num_points, a, b, c, d);
331 normalizeRecurrenceCoefficients(a, b, c, d);
341 normalizeRecurrenceCoefficients(a, b, c, d);
344 quad_points.
resize(num_points);
345 quad_weights.
resize(num_points);
347 if (num_points == 1) {
348 quad_points[0] = a[0];
349 quad_weights[0] = beta[0];
384 my_lapack.
PTEQR(
'I', num_points, &a[0], &b[1], eig_vectors.
values(),
385 num_points, &workspace[0], &info_flag);
387 "PTEQR returned info = " << info_flag);
392 quad_points[i] = a[num_points-1-i]-shift;
393 if (
std::abs(quad_points[i]) < quad_zero_tol)
394 quad_points[i] = 0.0;
395 quad_weights[i] = beta[0]*eig_vectors[num_points-1-i][0]*eig_vectors[num_points-1-i][0];
400 quad_values.resize(num_points);
402 quad_values[i].resize(p+1);
403 evaluateBases(quad_points[i], quad_values[i]);
407 template <
typename ordinal_type,
typename value_type>
415 template <
typename ordinal_type,
typename value_type>
427 template <
typename ordinal_type,
typename value_type>
443 template <
typename ordinal_type,
typename value_type>
457 template <
typename ordinal_type,
typename value_type>
472 b[k] =
std::sqrt((b[k]*g[k])/(c[k]*c[k-1]));
ScalarType * values() const
KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
virtual void evaluateBases(const value_type &point, Teuchos::Array< value_type > &basis_pts) const
Evaluate each basis polynomial at given point point.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeSparseTripleProductTensor(ordinal_type order) const
Compute triple product tensor.
virtual void getRecurrenceCoefficients(Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &delta, Teuchos::Array< value_type > &gamma) const
Return recurrence coefficients defined by above formula.
virtual const std::string & getName() const
Return string name of basis.
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 value_type evaluate(const value_type &point, ordinal_type order) const
Evaluate basis polynomial given by order order at given point point.
virtual ~RecurrenceBasis()
Destructor.
virtual ordinal_type pointGrowth(ordinal_type n) const
Evaluate point growth rule for Smolyak-type bases.
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...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
virtual Teuchos::RCP< Stokhos::Dense3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
void PTEQR(const char &COMPZ, const OrdinalType &n, MagnitudeType *D, MagnitudeType *E, ScalarType *Z, const OrdinalType &ldz, MagnitudeType *WORK, OrdinalType *info) const
void normalizeRecurrenceCoefficients(Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &delta, Teuchos::Array< value_type > &gamma) const
Normalize coefficients.
virtual void evaluateBasesAndDerivatives(const value_type &point, Teuchos::Array< value_type > &vals, Teuchos::Array< value_type > &derivs) const
Evaluate basis polynomials and their derivatives at given point point.
virtual ordinal_type quadDegreeOfExactness(ordinal_type n) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
KOKKOS_INLINE_FUNCTION PCE< Storage > ceil(const PCE< Storage > &a)
virtual Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > computeDerivDoubleProductTensor() const
Compute derivative double product tensor.
virtual ordinal_type order() const
Return order of basis (largest monomial degree ).
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
void resize(size_type new_size, const value_type &x=value_type())
virtual void print(std::ostream &os) const
Print basis to stream os.
virtual ordinal_type coefficientGrowth(ordinal_type n) const
Evaluate coefficient growth rule for Smolyak-type bases.
virtual void setup()
Setup basis after computing recurrence coefficients.
virtual ordinal_type size() const
Return total size of basis (given by order() + 1).
Data structure storing a dense 3-tensor C(i,j,k).
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
RecurrenceBasis(const std::string &name, ordinal_type p, bool normalize, GrowthPolicy growth=SLOW_GROWTH)
Constructor to be called by derived classes.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)