Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_RecurrenceBasisImp.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Teuchos_LAPACK.hpp"
12 
13 template <typename ordinal_type, typename value_type>
15 RecurrenceBasis(const std::string& name_, ordinal_type p_, bool normalize_,
16  Stokhos::GrowthPolicy growth_) :
17  name(name_),
18  p(p_),
19  normalize(normalize_),
20  growth(growth_),
21  quad_zero_tol(1.0e-14),
22 #ifdef HAVE_STOKHOS_DAKOTA
23  sparse_grid_growth_rule(webbur::level_to_order_linear_nn),
24 #else
25  sparse_grid_growth_rule(NULL),
26 #endif
27  alpha(p+1, value_type(0.0)),
28  beta(p+1, value_type(0.0)),
29  delta(p+1, value_type(0.0)),
30  gamma(p+1, value_type(0.0)),
31  norms(p+1, value_type(0.0))
32 {
33 }
34 
35 template <typename ordinal_type, typename value_type>
38  name(basis.name),
39  p(p_),
40  normalize(basis.normalize),
41  growth(basis.growth),
42  quad_zero_tol(basis.quad_zero_tol),
43  sparse_grid_growth_rule(basis.sparse_grid_growth_rule),
44  alpha(p+1, value_type(0.0)),
45  beta(p+1, value_type(0.0)),
46  delta(p+1, value_type(0.0)),
47  gamma(p+1, value_type(0.0)),
48  norms(p+1, value_type(0.0))
49 {
50 }
51 
52 template <typename ordinal_type, typename value_type>
53 void
56 {
57  bool is_normalized =
58  computeRecurrenceCoefficients(p+1, alpha, beta, delta, gamma);
59  if (normalize && !is_normalized) {
60  normalizeRecurrenceCoefficients(alpha, beta, delta, gamma);
61  }
62 
63 
64  // Compute norms -- when polynomials are normalized, this should work
65  // out to one (norms[0] == 1, delta[k] == 1, beta[k] == gamma[k]
66  norms[0] = beta[0]/(gamma[0]*gamma[0]);
67  for (ordinal_type k=1; k<=p; k++) {
68  norms[k] = (beta[k]/gamma[k])*(delta[k-1]/delta[k])*norms[k-1];
69  }
70 }
71 
72 template <typename ordinal_type, typename value_type>
75 {
76 }
77 
78 template <typename ordinal_type, typename value_type>
81 order() const
82 {
83  return p;
84 }
85 
86 template <typename ordinal_type, typename value_type>
89 size() const
90 {
91  return p+1;
92 }
93 
94 template <typename ordinal_type, typename value_type>
97 norm_squared() const
98 {
99  return norms;
100 }
101 
102 template <typename ordinal_type, typename value_type>
103 const value_type&
106 {
107  return norms[i];
108 }
109 
110 template <typename ordinal_type, typename value_type>
114 {
115  // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
116  ordinal_type sz = size();
119  Teuchos::Array<value_type> points, weights;
121  getQuadPoints(3*p, points, weights, values);
122 
123  for (ordinal_type i=0; i<sz; i++) {
124  for (ordinal_type j=0; j<sz; j++) {
125  for (ordinal_type k=0; k<sz; k++) {
126  value_type triple_product = 0;
127  for (ordinal_type l=0; l<static_cast<ordinal_type>(points.size());
128  l++){
129  triple_product +=
130  weights[l]*(values[l][i])*(values[l][j])*(values[l][k]);
131  }
132  (*Cijk)(i,j,k) = triple_product;
133  }
134  }
135  }
136 
137  return Cijk;
138 }
139 
140 template <typename ordinal_type, typename value_type>
144 {
145  // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
146  value_type sparse_tol = 1.0e-15;
147  ordinal_type sz = size();
150  Teuchos::Array<value_type> points, weights;
152  getQuadPoints(3*p, points, weights, values);
153 
154  for (ordinal_type i=0; i<sz; i++) {
155  for (ordinal_type j=0; j<sz; j++) {
156  for (ordinal_type k=0; k<theOrder; k++) {
157  value_type triple_product = 0;
158  for (ordinal_type l=0; l<static_cast<ordinal_type>(points.size());
159  l++){
160  triple_product +=
161  weights[l]*(values[l][i])*(values[l][j])*(values[l][k]);
162  }
163  if (std::abs(triple_product/norms[i]) > sparse_tol)
164  Cijk->add_term(i,j,k,triple_product);
165  }
166  }
167  }
168  Cijk->fillComplete();
169 
170  return Cijk;
171 }
172 
173 template <typename ordinal_type, typename value_type>
177 {
178  // Compute Bij = < \Psi_i' \Psi_j >
179  Teuchos::Array<value_type> points, weights;
181  getQuadPoints(2*p, points, weights, values);
182  ordinal_type nqp = weights.size();
183  derivs.resize(nqp);
184  ordinal_type sz = size();
185  for (ordinal_type i=0; i<nqp; i++) {
186  derivs[i].resize(sz);
187  evaluateBasesAndDerivatives(points[i], values[i], derivs[i]);
188  }
191  for (ordinal_type i=0; i<sz; i++) {
192  for (ordinal_type j=0; j<sz; j++) {
193  value_type b = value_type(0.0);
194  for (int qp=0; qp<nqp; qp++)
195  b += weights[qp]*derivs[qp][i]*values[qp][j];
196  (*Bij)(i,j) = b;
197  }
198  }
199 
200  return Bij;
201 }
202 
203 template <typename ordinal_type, typename value_type>
204 void
207 {
208  // Evaluate basis polynomials P(x) using 3 term recurrence
209  // P_0(x) = 1/gamma[0], P_-1 = 0
210  // P_i(x) = [(delta[i-1]*x-alpha[i-1])*P_{i-1}(x) -
211  // beta[i-1]*P_{i-2}(x)]/gamma[i],
212  // i=2,3,...
213  basis_pts[0] = value_type(1)/gamma[0];
214  if (p >= 1)
215  basis_pts[1] = (delta[0]*x-alpha[0])*basis_pts[0]/gamma[1];
216  for (ordinal_type i=2; i<=p; i++)
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];
219 }
220 
221 template <typename ordinal_type, typename value_type>
222 void
226  Teuchos::Array<value_type>& derivs) const
227 {
228  evaluateBases(x, vals);
229  derivs[0] = 0.0;
230  if (p >= 1)
231  derivs[1] = delta[0]/(gamma[0]*gamma[1]);
232  for (ordinal_type i=2; i<=p; i++)
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];
235 }
236 
237 template <typename ordinal_type, typename value_type>
240 evaluate(const value_type& x, ordinal_type k) const
241 {
242  // Evaluate basis polynomials P(x) using 3 term recurrence
243  // P_0(x) = 1/gamma[0], P_-1 = 0
244  // P_i(x) = [(delta[i-1]*x-alpha[i-1])*P_{i-1}(x) -
245  // beta[i-1]*P_{i-2}(x)]/gamma[i],
246  // i=2,3,...
247 
248  value_type v0 = value_type(1.0)/gamma[0];
249  if (k == 0)
250  return v0;
251 
252  value_type v1 = (x*delta[0]-alpha[0])*v0/gamma[1];
253  if (k == 1)
254  return v1;
255 
256  value_type v2 = value_type(0.0);
257  for (ordinal_type i=2; i<=k; i++) {
258  v2 = ((delta[i-1]*x-alpha[i-1])*v1 - beta[i-1]*v0)/gamma[i];
259  v0 = v1;
260  v1 = v2;
261  }
262 
263  return v2;
264 }
265 
266 template <typename ordinal_type, typename value_type>
267 void
269 print(std::ostream& os) const
270 {
271  os << name << " basis of order " << p << "." << std::endl;
272 
273  os << "Alpha recurrence coefficients:\n\t";
274  for (ordinal_type i=0; i<=p; i++)
275  os << alpha[i] << " ";
276  os << std::endl;
277 
278  os << "Beta recurrence coefficients:\n\t";
279  for (ordinal_type i=0; i<=p; i++)
280  os << beta[i] << " ";
281  os << std::endl;
282 
283  os << "Delta recurrence coefficients:\n\t";
284  for (ordinal_type i=0; i<=p; i++)
285  os << delta[i] << " ";
286  os << std::endl;
287 
288  os << "Gamma recurrence coefficients:\n\t";
289  for (ordinal_type i=0; i<=p; i++)
290  os << gamma[i] << " ";
291  os << std::endl;
292 
293  os << "Basis polynomial norms (squared):\n\t";
294  for (ordinal_type i=0; i<=p; i++)
295  os << norms[i] << " ";
296  os << std::endl;
297 }
298 
299 template <typename ordinal_type, typename value_type>
300 const std::string&
302 getName() const
303 {
304  return name;
305 }
306 
307 template <typename ordinal_type, typename value_type>
308 void
311  Teuchos::Array<value_type>& quad_points,
312  Teuchos::Array<value_type>& quad_weights,
313  Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
314 {
315 
316  //This is a transposition into C++ of Gautschi's code for taking the first
317  // N recurrance coefficient and generating a N point quadrature rule.
318  // The MATLAB version is available at
319  // http://www.cs.purdue.edu/archives/2002/wxg/codes/gauss.m
320  ordinal_type num_points =
321  static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
322  Teuchos::Array<value_type> a(num_points,0);
323  Teuchos::Array<value_type> b(num_points,0);
324  Teuchos::Array<value_type> c(num_points,0);
325  Teuchos::Array<value_type> d(num_points,0);
326 
327  // If we don't have enough recurrance coefficients, get some more.
328  if(num_points > p+1){
329  bool is_normalized = computeRecurrenceCoefficients(num_points, a, b, c, d);
330  if (!is_normalized)
331  normalizeRecurrenceCoefficients(a, b, c, d);
332  }
333  else { //else just take the ones we already have.
334  for(ordinal_type n = 0; n<num_points; n++){
335  a[n] = alpha[n];
336  b[n] = beta[n];
337  c[n] = delta[n];
338  d[n] = gamma[n];
339  }
340  if (!normalize)
341  normalizeRecurrenceCoefficients(a, b, c, d);
342  }
343 
344  quad_points.resize(num_points);
345  quad_weights.resize(num_points);
346 
347  if (num_points == 1) {
348  quad_points[0] = a[0];
349  quad_weights[0] = beta[0];
350  }
351  else {
352 
353  // With normalized coefficients, A is symmetric and tri-diagonal, with
354  // diagonal = a, and off-diagonal = b, starting at b[1]. We use
355  // LAPACK'S PTEQR function to compute the eigenvalues/vectors of A
356  // to compute the quadrature points/weights respectively. PTEQR requires
357  // an SPD matrix, so we need to shift the matrix to make it diagonally
358  // dominant (We could use STEQR which works for indefinite matrices, but
359  // this function has proven to be problematic on some platforms).
361  num_points);
362  Teuchos::Array<value_type> workspace(4*num_points);
363  ordinal_type info_flag;
365 
366  // Compute a shift to make matrix positive-definite
367  value_type max_a = 0.0;
368  value_type max_b = 0.0;
369  for (ordinal_type n = 0; n<num_points; n++) {
370  if (std::abs(a[n]) > max_a)
371  max_a = a[n];
372  }
373  for (ordinal_type n = 1; n<num_points; n++) {
374  if (std::abs(b[n]) > max_b)
375  max_b = b[n];
376  }
377  value_type shift = 1.0 + max_a + 2.0*max_b;
378 
379  // Shift diagonal
380  for (ordinal_type n = 0; n<num_points; n++)
381  a[n] += shift;
382 
383  // compute the eigenvalues (stored in a) and right eigenvectors.
384  my_lapack.PTEQR('I', num_points, &a[0], &b[1], eig_vectors.values(),
385  num_points, &workspace[0], &info_flag);
386  TEUCHOS_TEST_FOR_EXCEPTION(info_flag != 0, std::logic_error,
387  "PTEQR returned info = " << info_flag);
388 
389  // (shifted) eigenvalues are sorted in descending order by PTEQR.
390  // Reorder them to ascending as in STEQR
391  for (ordinal_type i=0; i<num_points; i++) {
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];
396  }
397  }
398 
399  // Evalute basis at gauss points
400  quad_values.resize(num_points);
401  for (ordinal_type i=0; i<num_points; i++) {
402  quad_values[i].resize(p+1);
403  evaluateBases(quad_points[i], quad_values[i]);
404  }
405 }
406 
407 template <typename ordinal_type, typename value_type>
411 {
412  return ordinal_type(2)*n-ordinal_type(1);
413 }
414 
415 template <typename ordinal_type, typename value_type>
419 {
420  if (growth == SLOW_GROWTH)
421  return n;
422 
423  // else moderate
424  return ordinal_type(2)*n;
425 }
426 
427 template <typename ordinal_type, typename value_type>
431 {
432  if (growth == SLOW_GROWTH) {
433  if (n % ordinal_type(2) == ordinal_type(1))
434  return n + ordinal_type(1);
435  else
436  return n;
437  }
438 
439  // else moderate
440  return n;
441 }
442 
443 template <typename ordinal_type, typename value_type>
444 void
450 {
451  a = alpha;
452  b = beta;
453  c = delta;
454  g = gamma;
455 }
456 
457 template <typename ordinal_type, typename value_type>
458 void
465 {
466  ordinal_type n = a.size();
467  a[0] = a[0] / c[0];
468  b[0] = std::sqrt(b[0]);
469  g[0] = b[0];
470  for (ordinal_type k=1; k<n; k++) {
471  a[k] = a[k] / c[k];
472  b[k] = std::sqrt((b[k]*g[k])/(c[k]*c[k-1]));
473  g[k] = b[k];
474  }
475  for (ordinal_type k=0; k<n; k++)
476  c[k] = value_type(1);
477 }
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.
size_type size() const
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.
int n
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)