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>
303 print(std::ostream& os)
const
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