44#include "Teuchos_LAPACK.hpp"
45#include "Teuchos_SerialDenseMatrix.hpp"
47template <
typename ordinal_type,
typename value_type>
49RecurrenceBasis(
const std::string& name_, ordinal_type p_,
bool normalize_,
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),
61 alpha(p+1, value_type(0.0)),
62 beta(p+1, value_type(0.0)),
63 delta(p+1, value_type(0.0)),
64 gamma(p+1, value_type(0.0)),
65 norms(p+1, value_type(0.0))
69template <
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),
78 alpha(p+1, value_type(0.0)),
79 beta(p+1, value_type(0.0)),
80 delta(p+1, value_type(0.0)),
81 gamma(p+1, value_type(0.0)),
82 norms(p+1, value_type(0.0))
86template <
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]);
101 for (ordinal_type k=1; k<=p; k++) {
102 norms[k] = (beta[k]/gamma[k])*(delta[k-1]/delta[k])*norms[k-1];
106template <
typename ordinal_type,
typename value_type>
112template <
typename ordinal_type,
typename value_type>
120template <
typename ordinal_type,
typename value_type>
128template <
typename ordinal_type,
typename value_type>
129const Teuchos::Array<value_type>&
136template <
typename ordinal_type,
typename value_type>
144template <
typename ordinal_type,
typename value_type>
145Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> >
150 ordinal_type sz = size();
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);
157 for (ordinal_type i=0; i<sz; i++) {
158 for (ordinal_type
j=0;
j<sz;
j++) {
159 for (ordinal_type k=0; k<sz; k++) {
160 value_type triple_product = 0;
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;
174template <
typename ordinal_type,
typename value_type>
175Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
180 value_type sparse_tol = 1.0e-15;
181 ordinal_type sz = size();
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);
188 for (ordinal_type i=0; i<sz; i++) {
189 for (ordinal_type
j=0;
j<sz;
j++) {
190 for (ordinal_type k=0; k<theOrder; k++) {
191 value_type triple_product = 0;
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();
207template <
typename ordinal_type,
typename value_type>
208Teuchos::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);
216 ordinal_type nqp = weights.size();
218 ordinal_type sz = size();
219 for (ordinal_type i=0; i<nqp; i++) {
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));
225 for (ordinal_type i=0; i<sz; i++) {
226 for (ordinal_type
j=0;
j<sz;
j++) {
227 value_type b = value_type(0.0);
228 for (
int qp=0; qp<nqp; qp++)
229 b += weights[qp]*derivs[qp][i]*values[qp][
j];
237template <
typename ordinal_type,
typename value_type>
240evaluateBases(
const value_type& x, Teuchos::Array<value_type>& basis_pts)
const
247 basis_pts[0] = value_type(1)/gamma[0];
249 basis_pts[1] = (delta[0]*x-alpha[0])*basis_pts[0]/gamma[1];
250 for (ordinal_type i=2; i<=p; i++)
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];
255template <
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]);
266 for (ordinal_type i=2; i<=p; i++)
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];
271template <
typename ordinal_type,
typename value_type>
274evaluate(
const value_type& x, ordinal_type k)
const
282 value_type v0 = value_type(1.0)/gamma[0];
286 value_type v1 = (x*delta[0]-alpha[0])*v0/gamma[1];
290 value_type v2 = value_type(0.0);
291 for (ordinal_type i=2; i<=k; i++) {
292 v2 = ((delta[i-1]*x-alpha[i-1])*v1 - beta[i-1]*v0)/gamma[i];
300template <
typename ordinal_type,
typename value_type>
303print(std::ostream& os)
const
305 os << name <<
" basis of order " << p <<
"." << std::endl;
307 os <<
"Alpha recurrence coefficients:\n\t";
308 for (ordinal_type i=0; i<=p; i++)
309 os << alpha[i] <<
" ";
312 os <<
"Beta recurrence coefficients:\n\t";
313 for (ordinal_type i=0; i<=p; i++)
314 os << beta[i] <<
" ";
317 os <<
"Delta recurrence coefficients:\n\t";
318 for (ordinal_type i=0; i<=p; i++)
319 os << delta[i] <<
" ";
322 os <<
"Gamma recurrence coefficients:\n\t";
323 for (ordinal_type i=0; i<=p; i++)
324 os << gamma[i] <<
" ";
327 os <<
"Basis polynomial norms (squared):\n\t";
328 for (ordinal_type i=0; i<=p; i++)
329 os << norms[i] <<
" ";
333template <
typename ordinal_type,
typename value_type>
341template <
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
354 ordinal_type num_points =
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);
368 for(ordinal_type n = 0; n<num_points; n++){
375 normalizeRecurrenceCoefficients(a, b, c, d);
378 quad_points.resize(num_points);
379 quad_weights.resize(num_points);
381 if (num_points == 1) {
382 quad_points[0] = a[0];
383 quad_weights[0] = beta[0];
394 Teuchos::SerialDenseMatrix<ordinal_type,value_type> eig_vectors(num_points,
396 Teuchos::Array<value_type> workspace(4*num_points);
397 ordinal_type info_flag;
401 value_type max_a = 0.0;
402 value_type max_b = 0.0;
403 for (ordinal_type n = 0; n<num_points; n++) {
404 if (std::abs(a[n]) > max_a)
407 for (ordinal_type n = 1; n<num_points; n++) {
408 if (std::abs(b[n]) > max_b)
411 value_type shift = 1.0 + max_a + 2.0*max_b;
414 for (ordinal_type n = 0; n<num_points; n++)
418 my_lapack.PTEQR(
'I', num_points, &a[0], &b[1], eig_vectors.values(),
419 num_points, &workspace[0], &info_flag);
420 TEUCHOS_TEST_FOR_EXCEPTION(info_flag != 0, std::logic_error,
421 "PTEQR returned info = " << info_flag);
425 for (ordinal_type i=0; i<num_points; i++) {
426 quad_points[i] = a[num_points-1-i]-shift;
427 if (std::abs(quad_points[i]) < quad_zero_tol)
428 quad_points[i] = 0.0;
429 quad_weights[i] = beta[0]*eig_vectors[num_points-1-i][0]*eig_vectors[num_points-1-i][0];
434 quad_values.resize(num_points);
435 for (ordinal_type i=0; i<num_points; i++) {
436 quad_values[i].resize(p+1);
437 evaluateBases(quad_points[i], quad_values[i]);
441template <
typename ordinal_type,
typename value_type>
446 return ordinal_type(2)*n-ordinal_type(1);
449template <
typename ordinal_type,
typename value_type>
458 return ordinal_type(2)*n;
461template <
typename ordinal_type,
typename value_type>
467 if (n % ordinal_type(2) == ordinal_type(1))
468 return n + ordinal_type(1);
477template <
typename ordinal_type,
typename value_type>
481 Teuchos::Array<value_type>& b,
482 Teuchos::Array<value_type>& c,
483 Teuchos::Array<value_type>&
g)
const
491template <
typename ordinal_type,
typename value_type>
495 Teuchos::Array<value_type>& a,
496 Teuchos::Array<value_type>& b,
497 Teuchos::Array<value_type>& c,
498 Teuchos::Array<value_type>&
g)
const
500 ordinal_type n = a.size();
502 b[0] = std::sqrt(b[0]);
504 for (ordinal_type k=1; k<n; k++) {
506 b[k] = std::sqrt((b[k]*
g[k])/(c[k]*c[k-1]));
509 for (ordinal_type k=0; k<n; k++)
510 c[k] = value_type(1);
Data structure storing a dense 3-tensor C(i,j,k).
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship:
virtual ordinal_type order() const
Return order of basis (largest monomial degree ).
virtual ordinal_type pointGrowth(ordinal_type n) const
Evaluate point growth rule for Smolyak-type bases.
virtual Teuchos::RCP< Stokhos::Dense3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual void print(std::ostream &os) const
Print basis to stream os.
virtual value_type evaluate(const value_type &point, ordinal_type order) const
Evaluate basis polynomial given by order order at given point point.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
virtual ordinal_type quadDegreeOfExactness(ordinal_type n) const
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeSparseTripleProductTensor(ordinal_type order) const
Compute triple product tensor.
virtual Teuchos::RCP< Teuchos::SerialDenseMatrix< ordinal_type, value_type > > computeDerivDoubleProductTensor() const
Compute derivative double product tensor.
virtual ordinal_type size() const
Return total size of basis (given by order() + 1).
virtual void evaluateBases(const value_type &point, Teuchos::Array< value_type > &basis_pts) const
Evaluate each basis polynomial at given point point.
virtual ordinal_type coefficientGrowth(ordinal_type n) const
Evaluate coefficient growth rule for Smolyak-type bases.
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 void setup()
Setup basis after computing recurrence coefficients.
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.
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 ~RecurrenceBasis()
Destructor.
RecurrenceBasis(const std::string &name, ordinal_type p, bool normalize, GrowthPolicy growth=SLOW_GROWTH)
Constructor to be called by derived classes.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Specialization for Sacado::UQ::PCE< Storage<...> >
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
GrowthPolicy
Enumerated type for determining Smolyak growth policies.