Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_RecurrenceBasisImp.hpp
Go to the documentation of this file.
1// $Id$
2// $Source$
3// @HEADER
4// ***********************************************************************
5//
6// Stokhos Package
7// Copyright (2009) Sandia Corporation
8//
9// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
10// license for use of this work by or on behalf of the U.S. Government.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
40//
41// ***********************************************************************
42// @HEADER
43
44#include "Teuchos_LAPACK.hpp"
45#include "Teuchos_SerialDenseMatrix.hpp"
46
47template <typename ordinal_type, typename value_type>
49RecurrenceBasis(const std::string& name_, ordinal_type p_, bool normalize_,
50 Stokhos::GrowthPolicy growth_) :
51 name(name_),
52 p(p_),
53 normalize(normalize_),
54 growth(growth_),
55 quad_zero_tol(1.0e-14),
56#ifdef HAVE_STOKHOS_DAKOTA
57 sparse_grid_growth_rule(webbur::level_to_order_linear_nn),
58#else
59 sparse_grid_growth_rule(NULL),
60#endif
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))
66{
67}
68
69template <typename ordinal_type, typename value_type>
71RecurrenceBasis(ordinal_type p_, const RecurrenceBasis& basis) :
72 name(basis.name),
73 p(p_),
74 normalize(basis.normalize),
75 growth(basis.growth),
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))
83{
84}
85
86template <typename ordinal_type, typename value_type>
87void
89setup()
91 bool is_normalized =
92 computeRecurrenceCoefficients(p+1, alpha, beta, delta, gamma);
93 if (normalize && !is_normalized) {
94 normalizeRecurrenceCoefficients(alpha, beta, delta, gamma);
95 }
97
98 // Compute norms -- when polynomials are normalized, this should work
99 // out to one (norms[0] == 1, delta[k] == 1, beta[k] == gamma[k]
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];
103 }
104}
105
106template <typename ordinal_type, typename value_type>
110}
111
112template <typename ordinal_type, typename value_type>
113ordinal_type
119
120template <typename ordinal_type, typename value_type>
121ordinal_type
123size() const
124{
125 return p+1;
126}
127
128template <typename ordinal_type, typename value_type>
129const Teuchos::Array<value_type>&
135
136template <typename ordinal_type, typename value_type>
137const value_type&
139norm_squared(ordinal_type i) const
140{
141 return norms[i];
142}
143
144template <typename ordinal_type, typename value_type>
145Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> >
148{
149 // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
150 ordinal_type sz = size();
151 Teuchos::RCP< Stokhos::Dense3Tensor<ordinal_type, value_type> > Cijk =
152 Teuchos::rcp(new Dense3Tensor<ordinal_type, value_type>(sz));
153 Teuchos::Array<value_type> points, weights;
154 Teuchos::Array< Teuchos::Array<value_type> > values;
155 getQuadPoints(3*p, points, weights, values);
156
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());
162 l++){
163 triple_product +=
164 weights[l]*(values[l][i])*(values[l][j])*(values[l][k]);
165 }
166 (*Cijk)(i,j,k) = triple_product;
167 }
169 }
170
171 return Cijk;
172}
173
174template <typename ordinal_type, typename value_type>
175Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
177computeSparseTripleProductTensor(ordinal_type theOrder) const
178{
179 // Compute Cijk = < \Psi_i \Psi_j \Psi_k >
180 value_type sparse_tol = 1.0e-15;
181 ordinal_type sz = size();
182 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
183 Teuchos::rcp(new Sparse3Tensor<ordinal_type, value_type>());
184 Teuchos::Array<value_type> points, weights;
185 Teuchos::Array< Teuchos::Array<value_type> > values;
186 getQuadPoints(3*p, points, weights, values);
187
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());
193 l++){
194 triple_product +=
195 weights[l]*(values[l][i])*(values[l][j])*(values[l][k]);
196 }
197 if (std::abs(triple_product/norms[i]) > sparse_tol)
198 Cijk->add_term(i,j,k,triple_product);
199 }
201 }
202 Cijk->fillComplete();
204 return Cijk;
205}
207template <typename ordinal_type, typename value_type>
208Teuchos::RCP< Teuchos::SerialDenseMatrix<ordinal_type, value_type> >
211{
212 // Compute Bij = < \Psi_i' \Psi_j >
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();
217 derivs.resize(nqp);
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]);
222 }
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];
230 (*Bij)(i,j) = b;
232 }
233
234 return Bij;
235}
236
237template <typename ordinal_type, typename value_type>
238void
240evaluateBases(const value_type& x, Teuchos::Array<value_type>& basis_pts) 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 basis_pts[0] = value_type(1)/gamma[0];
248 if (p >= 1)
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];
253}
254
255template <typename ordinal_type, typename value_type>
256void
258evaluateBasesAndDerivatives(const value_type& x,
259 Teuchos::Array<value_type>& vals,
260 Teuchos::Array<value_type>& derivs) const
261{
262 evaluateBases(x, vals);
263 derivs[0] = 0.0;
264 if (p >= 1)
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];
269}
270
271template <typename ordinal_type, typename value_type>
272value_type
274evaluate(const value_type& x, ordinal_type k) const
275{
276 // Evaluate basis polynomials P(x) using 3 term recurrence
277 // P_0(x) = 1/gamma[0], P_-1 = 0
278 // P_i(x) = [(delta[i-1]*x-alpha[i-1])*P_{i-1}(x) -
279 // beta[i-1]*P_{i-2}(x)]/gamma[i],
280 // i=2,3,...
281
282 value_type v0 = value_type(1.0)/gamma[0];
283 if (k == 0)
284 return v0;
285
286 value_type v1 = (x*delta[0]-alpha[0])*v0/gamma[1];
287 if (k == 1)
288 return v1;
289
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];
293 v0 = v1;
294 v1 = v2;
295 }
296
297 return v2;
298}
299
300template <typename ordinal_type, typename value_type>
301void
303print(std::ostream& os) const
304{
305 os << name << " basis of order " << p << "." << std::endl;
306
307 os << "Alpha recurrence coefficients:\n\t";
308 for (ordinal_type i=0; i<=p; i++)
309 os << alpha[i] << " ";
310 os << std::endl;
311
312 os << "Beta recurrence coefficients:\n\t";
313 for (ordinal_type i=0; i<=p; i++)
314 os << beta[i] << " ";
315 os << std::endl;
316
317 os << "Delta recurrence coefficients:\n\t";
318 for (ordinal_type i=0; i<=p; i++)
319 os << delta[i] << " ";
320 os << std::endl;
321
322 os << "Gamma recurrence coefficients:\n\t";
323 for (ordinal_type i=0; i<=p; i++)
324 os << gamma[i] << " ";
325 os << std::endl;
326
327 os << "Basis polynomial norms (squared):\n\t";
328 for (ordinal_type i=0; i<=p; i++)
329 os << norms[i] << " ";
330 os << std::endl;
331}
332
333template <typename ordinal_type, typename value_type>
334const std::string&
340
341template <typename ordinal_type, typename value_type>
342void
344getQuadPoints(ordinal_type quad_order,
345 Teuchos::Array<value_type>& quad_points,
346 Teuchos::Array<value_type>& quad_weights,
347 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
348{
349
350 //This is a transposition into C++ of Gautschi's code for taking the first
351 // N recurrance coefficient and generating a N point quadrature rule.
352 // The MATLAB version is available at
353 // http://www.cs.purdue.edu/archives/2002/wxg/codes/gauss.m
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);
360
361 // If we don't have enough recurrance coefficients, get some more.
362 if(num_points > p+1){
363 bool is_normalized = computeRecurrenceCoefficients(num_points, a, b, c, d);
364 if (!is_normalized)
365 normalizeRecurrenceCoefficients(a, b, c, d);
366 }
367 else { //else just take the ones we already have.
368 for(ordinal_type n = 0; n<num_points; n++){
369 a[n] = alpha[n];
370 b[n] = beta[n];
371 c[n] = delta[n];
372 d[n] = gamma[n];
373 }
374 if (!normalize)
375 normalizeRecurrenceCoefficients(a, b, c, d);
376 }
377
378 quad_points.resize(num_points);
379 quad_weights.resize(num_points);
380
381 if (num_points == 1) {
382 quad_points[0] = a[0];
383 quad_weights[0] = beta[0];
384 }
385 else {
386
387 // With normalized coefficients, A is symmetric and tri-diagonal, with
388 // diagonal = a, and off-diagonal = b, starting at b[1]. We use
389 // LAPACK'S PTEQR function to compute the eigenvalues/vectors of A
390 // to compute the quadrature points/weights respectively. PTEQR requires
391 // an SPD matrix, so we need to shift the matrix to make it diagonally
392 // dominant (We could use STEQR which works for indefinite matrices, but
393 // this function has proven to be problematic on some platforms).
394 Teuchos::SerialDenseMatrix<ordinal_type,value_type> eig_vectors(num_points,
395 num_points);
396 Teuchos::Array<value_type> workspace(4*num_points);
397 ordinal_type info_flag;
399
400 // Compute a shift to make matrix positive-definite
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)
405 max_a = a[n];
406 }
407 for (ordinal_type n = 1; n<num_points; n++) {
408 if (std::abs(b[n]) > max_b)
409 max_b = b[n];
410 }
411 value_type shift = 1.0 + max_a + 2.0*max_b;
412
413 // Shift diagonal
414 for (ordinal_type n = 0; n<num_points; n++)
415 a[n] += shift;
416
417 // compute the eigenvalues (stored in a) and right eigenvectors.
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);
422
423 // (shifted) eigenvalues are sorted in descending order by PTEQR.
424 // Reorder them to ascending as in STEQR
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];
430 }
431 }
432
433 // Evalute basis at gauss points
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]);
438 }
439}
440
441template <typename ordinal_type, typename value_type>
442ordinal_type
444quadDegreeOfExactness(ordinal_type n) const
445{
446 return ordinal_type(2)*n-ordinal_type(1);
447}
448
449template <typename ordinal_type, typename value_type>
450ordinal_type
452coefficientGrowth(ordinal_type n) const
453{
454 if (growth == SLOW_GROWTH)
455 return n;
456
457 // else moderate
458 return ordinal_type(2)*n;
459}
460
461template <typename ordinal_type, typename value_type>
462ordinal_type
464pointGrowth(ordinal_type n) const
465{
466 if (growth == SLOW_GROWTH) {
467 if (n % ordinal_type(2) == ordinal_type(1))
468 return n + ordinal_type(1);
469 else
470 return n;
471 }
472
473 // else moderate
474 return n;
475}
476
477template <typename ordinal_type, typename value_type>
478void
480getRecurrenceCoefficients(Teuchos::Array<value_type>& a,
481 Teuchos::Array<value_type>& b,
482 Teuchos::Array<value_type>& c,
483 Teuchos::Array<value_type>& g) const
484{
485 a = alpha;
486 b = beta;
487 c = delta;
488 g = gamma;
489}
490
491template <typename ordinal_type, typename value_type>
492void
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
499{
500 ordinal_type n = a.size();
501 a[0] = a[0] / c[0];
502 b[0] = std::sqrt(b[0]);
503 g[0] = b[0];
504 for (ordinal_type k=1; k<n; k++) {
505 a[k] = a[k] / c[k];
506 b[k] = std::sqrt((b[k]*g[k])/(c[k]*c[k-1]));
507 g[k] = b[k];
508 }
509 for (ordinal_type k=0; k<n; k++)
510 c[k] = value_type(1);
511}
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.
RecurrenceBasis(const std::string &name, ordinal_type p, bool normalize, GrowthPolicy growth=SLOW_GROWTH)
Constructor to be called by derived classes.
Specialization for Sacado::UQ::PCE< Storage<...> >
ScalarType g(const Teuchos::Array< ScalarType > &x, const ScalarType &y)
GrowthPolicy
Enumerated type for determining Smolyak growth policies.