44 #include "Teuchos_LAPACK.hpp" 45 #include "Teuchos_SerialDenseMatrix.hpp" 47 template <
typename ordinal_type,
typename value_type>
53 normalize(normalize_),
55 quad_zero_tol(1.0e-14),
56 #ifdef HAVE_STOKHOS_DAKOTA
57 sparse_grid_growth_rule(webbur::level_to_order_linear_nn),
59 sparse_grid_growth_rule(NULL),
69 template <
typename ordinal_type,
typename value_type>
74 normalize(basis.normalize),
76 quad_zero_tol(basis.quad_zero_tol),
77 sparse_grid_growth_rule(basis.sparse_grid_growth_rule),
86 template <
typename ordinal_type,
typename value_type>
92 computeRecurrenceCoefficients(p+1, alpha, beta, delta, gamma);
93 if (normalize && !is_normalized) {
94 normalizeRecurrenceCoefficients(alpha, beta, delta, gamma);
100 norms[0] = beta[0]/(gamma[0]*gamma[0]);
102 norms[k] = (beta[k]/gamma[k])*(delta[k-1]/delta[k])*norms[k-1];
106 template <
typename ordinal_type,
typename value_type>
112 template <
typename ordinal_type,
typename value_type>
120 template <
typename ordinal_type,
typename value_type>
128 template <
typename ordinal_type,
typename value_type>
129 const Teuchos::Array<value_type>&
136 template <
typename ordinal_type,
typename value_type>
144 template <
typename ordinal_type,
typename value_type>
145 Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> >
151 Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> > Cijk =
153 Teuchos::Array<value_type> points, weights;
154 Teuchos::Array< Teuchos::Array<value_type> > values;
155 getQuadPoints(3*p, points, weights, values);
161 for (
ordinal_type l=0; l<static_cast<ordinal_type>(points.size());
164 weights[l]*(values[l][i])*(values[l][
j])*(values[l][k]);
166 (*Cijk)(i,
j,k) = triple_product;
174 template <
typename ordinal_type,
typename value_type>
175 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
182 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
184 Teuchos::Array<value_type> points, weights;
185 Teuchos::Array< Teuchos::Array<value_type> > values;
186 getQuadPoints(3*p, points, weights, values);
192 for (
ordinal_type l=0; l<static_cast<ordinal_type>(points.size());
195 weights[l]*(values[l][i])*(values[l][
j])*(values[l][k]);
197 if (
std::abs(triple_product/norms[i]) > sparse_tol)
198 Cijk->add_term(i,
j,k,triple_product);
202 Cijk->fillComplete();
207 template <
typename ordinal_type,
typename value_type>
208 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> >
213 Teuchos::Array<value_type> points, weights;
214 Teuchos::Array< Teuchos::Array<value_type> > values, derivs;
215 getQuadPoints(2*p, points, weights, values);
220 derivs[i].resize(sz);
221 evaluateBasesAndDerivatives(points[i], values[i], derivs[i]);
223 Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> > Bij =
224 Teuchos::rcp(
new Teuchos::SerialDenseMatrix<ordinal_type, value_type>(sz,sz));
228 for (
int qp=0; qp<nqp; qp++)
229 b += weights[qp]*derivs[qp][i]*values[qp][
j];
237 template <
typename ordinal_type,
typename value_type>
249 basis_pts[1] = (delta[0]*
x-alpha[0])*basis_pts[0]/gamma[1];
251 basis_pts[i] = ((delta[i-1]*
x-alpha[i-1])*basis_pts[i-1] -
252 beta[i-1]*basis_pts[i-2])/gamma[i];
255 template <
typename ordinal_type,
typename value_type>
259 Teuchos::Array<value_type>& vals,
260 Teuchos::Array<value_type>& derivs)
const 262 evaluateBases(
x, vals);
265 derivs[1] = delta[0]/(gamma[0]*gamma[1]);
267 derivs[i] = (delta[i-1]*vals[i-1] + (delta[i-1]*
x-alpha[i-1])*derivs[i-1] -
268 beta[i-1]*derivs[i-2])/gamma[i];
271 template <
typename ordinal_type,
typename value_type>
292 v2 = ((delta[i-1]*
x-alpha[i-1])*v1 - beta[i-1]*v0)/gamma[i];
300 template <
typename ordinal_type,
typename value_type>
305 os << name <<
" basis of order " << p <<
"." << std::endl;
307 os <<
"Alpha recurrence coefficients:\n\t";
309 os << alpha[i] <<
" ";
312 os <<
"Beta recurrence coefficients:\n\t";
314 os << beta[i] <<
" ";
317 os <<
"Delta recurrence coefficients:\n\t";
319 os << delta[i] <<
" ";
322 os <<
"Gamma recurrence coefficients:\n\t";
324 os << gamma[i] <<
" ";
327 os <<
"Basis polynomial norms (squared):\n\t";
329 os << norms[i] <<
" ";
333 template <
typename ordinal_type,
typename value_type>
341 template <
typename ordinal_type,
typename value_type>
345 Teuchos::Array<value_type>& quad_points,
346 Teuchos::Array<value_type>& quad_weights,
347 Teuchos::Array< Teuchos::Array<value_type> >& quad_values)
const 355 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
356 Teuchos::Array<value_type> a(num_points,0);
357 Teuchos::Array<value_type> b(num_points,0);
358 Teuchos::Array<value_type> c(num_points,0);
359 Teuchos::Array<value_type> d(num_points,0);
362 if(num_points > p+1){
363 bool is_normalized = computeRecurrenceCoefficients(num_points, a, b, c, d);
365 normalizeRecurrenceCoefficients(a, b, c, d);
375 normalizeRecurrenceCoefficients(a, b, c, d);
380 Teuchos::SerialDenseMatrix<ordinal_type,value_type> eig_vectors(num_points,
382 Teuchos::Array<value_type> workspace(2*num_points);
388 my_lapack.STEQR(
'I', num_points, &a[0], &b[0], eig_vectors.values(),
389 num_points, &workspace[0], &info_flag);
391 my_lapack.STEQR(
'I', num_points, &a[0], &b[1], eig_vectors.values(),
392 num_points, &workspace[0], &info_flag);
395 quad_points.resize(num_points);
396 quad_weights.resize(num_points);
398 quad_points[i] = a[i];
399 if (
std::abs(quad_points[i]) < quad_zero_tol)
400 quad_points[i] = 0.0;
401 quad_weights[i] = beta[0]*eig_vectors[i][0]*eig_vectors[i][0];
405 quad_values.resize(num_points);
407 quad_values[i].resize(p+1);
408 evaluateBases(quad_points[i], quad_values[i]);
412 template <
typename ordinal_type,
typename value_type>
420 template <
typename ordinal_type,
typename value_type>
432 template <
typename ordinal_type,
typename value_type>
448 template <
typename ordinal_type,
typename value_type>
452 Teuchos::Array<value_type>& b,
453 Teuchos::Array<value_type>& c,
454 Teuchos::Array<value_type>&
g)
const 462 template <
typename ordinal_type,
typename value_type>
466 Teuchos::Array<value_type>& a,
467 Teuchos::Array<value_type>& b,
468 Teuchos::Array<value_type>& c,
469 Teuchos::Array<value_type>&
g)
const KOKKOS_INLINE_FUNCTION PCE< Storage > sqrt(const PCE< Storage > &a)
virtual value_type evaluate(const value_type &point, ordinal_type order) const
Evaluate basis polynomial given by order order at given point point.
virtual ordinal_type pointGrowth(ordinal_type n) const
Evaluate point growth rule for Smolyak-type bases.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship: for ...
virtual ordinal_type coefficientGrowth(ordinal_type n) const
Evaluate coefficient growth rule for Smolyak-type bases.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
virtual ~RecurrenceBasis()
Destructor.
virtual ordinal_type size() const
Return total size of basis (given by order() + 1).
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
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...
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
virtual Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > computeDerivDoubleProductTensor() const
Compute derivative double product tensor.
virtual void print(std::ostream &os) const
Print basis to stream os.
virtual const std::string & getName() const
Return string name of basis.
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
virtual Teuchos::RCP< Stokhos::Dense3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual void setup()
Setup basis after computing recurrence coefficients.
virtual ordinal_type order() const
Return order of basis (largest monomial degree ).
Data structure storing a dense 3-tensor C(i,j,k).
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 Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeSparseTripleProductTensor(ordinal_type order) const
Compute triple product tensor.
RecurrenceBasis(const std::string &name, ordinal_type p, bool normalize, GrowthPolicy growth=SLOW_GROWTH)
Constructor to be called by derived classes.
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
virtual void evaluateBases(const value_type &point, Teuchos::Array< value_type > &basis_pts) const
Evaluate each basis polynomial at given point point.
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)