Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_StieltjesBasisImp.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Stokhos Package
5// Copyright (2009) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38//
39// ***********************************************************************
40// @HEADER
41
42#include "Teuchos_Assert.hpp"
43#include "Teuchos_TimeMonitor.hpp"
44
45template <typename ordinal_type, typename value_type, typename func_type>
48 ordinal_type p,
49 const Teuchos::RCP<const func_type>& func_,
50 const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad_,
51 bool use_pce_quad_points_,
52 bool normalize) :
53 RecurrenceBasis<ordinal_type, value_type>("Stieltjes PCE", p, normalize),
54 func(func_),
55 quad(quad_),
56 pce_weights(quad->getQuadWeights()),
57 basis_values(quad->getBasisAtQuadPoints()),
58 func_vals(),
59 phi_vals(),
60 use_pce_quad_points(use_pce_quad_points_)
61{
62 // Setup rest of recurrence basis
63 this->setup();
64}
65
66template <typename ordinal_type, typename value_type, typename func_type>
69{
70}
71
72template <typename ordinal_type, typename value_type, typename func_type>
73void
75getQuadPoints(ordinal_type quad_order,
76 Teuchos::Array<value_type>& quad_points,
77 Teuchos::Array<value_type>& quad_weights,
78 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
79{
80#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
81 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::StieltjesBasis -- compute Gauss points");
82#endif
83
84 // Use underlying pce's quad points, weights, values
85 if (use_pce_quad_points) {
86 quad_points = func_vals;
87 quad_weights = pce_weights;
88 quad_values = phi_vals;
89 return;
90 }
91
92 // Call base class
93 ordinal_type num_points =
94 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
95
96 // We can't reliably generate quadrature points of order > 2*p
97 if (quad_order > 2*this->p)
98 quad_order = 2*this->p;
100 quad_points,
101 quad_weights,
102 quad_values);
103
104 // Fill in the rest of the points with zero weight
105 if (quad_weights.size() < num_points) {
106 ordinal_type old_size = quad_weights.size();
107 quad_weights.resize(num_points);
108 quad_points.resize(num_points);
109 quad_values.resize(num_points);
110 for (ordinal_type i=old_size; i<num_points; i++) {
111 quad_weights[i] = value_type(0);
112 quad_points[i] = quad_points[0];
113 quad_values[i].resize(this->p+1);
114 this->evaluateBases(quad_points[i], quad_values[i]);
115 }
116 }
117}
118
119template <typename ordinal_type, typename value_type, typename func_type>
120bool
122computeRecurrenceCoefficients(ordinal_type n,
123 Teuchos::Array<value_type>& alpha,
124 Teuchos::Array<value_type>& beta,
125 Teuchos::Array<value_type>& delta,
126 Teuchos::Array<value_type>& gamma) const
127{
128 ordinal_type nqp = phi_vals.size();
129 Teuchos::Array<value_type> nrm(n);
130 Teuchos::Array< Teuchos::Array<value_type> > vals(nqp);
131 for (ordinal_type i=0; i<nqp; i++)
132 vals[i].resize(n);
133 stieltjes(0, n, pce_weights, func_vals, alpha, beta, nrm, vals);
134 for (ordinal_type i=0; i<n; i++) {
135 delta[i] = value_type(1.0);
136 gamma[i] = value_type(1.0);
137 }
138
139 // Save basis functions at quad point values
140 if (n == this->p+1)
141 phi_vals = vals;
142
143 return false;
144}
145
146template <typename ordinal_type, typename value_type, typename func_type>
147void
149setup()
150{
151 // Evaluate func at quad points
152 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
153 quad->getQuadPoints();
154 ordinal_type nqp = pce_weights.size();
155 func_vals.resize(nqp);
156 phi_vals.resize(nqp);
157 for (ordinal_type i=0; i<nqp; i++) {
158 func_vals[i] = (*func)(quad_points[i]);
159 phi_vals[i].resize(this->p+1);
160 }
161
163}
164
165template <typename ordinal_type, typename value_type, typename func_type>
166void
168stieltjes(ordinal_type nstart,
169 ordinal_type nfinish,
170 const Teuchos::Array<value_type>& weights,
171 const Teuchos::Array<value_type>& points,
172 Teuchos::Array<value_type>& a,
173 Teuchos::Array<value_type>& b,
174 Teuchos::Array<value_type>& nrm,
175 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals) const
176{
177#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
178 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::StieltjesBasis -- Discretized Stieltjes Procedure");
179#endif
180
181 value_type val1, val2;
182 ordinal_type start = nstart;
183 if (nstart == 0) {
184 integrateBasisSquared(0, a, b, weights, points, phi_vals, val1, val2);
185 nrm[0] = val1;
186 a[0] = val2/val1;
187 b[0] = value_type(1);
188 start = 1;
189 }
190 for (ordinal_type i=start; i<nfinish; i++) {
191 integrateBasisSquared(i, a, b, weights, points, phi_vals, val1, val2);
192 // std::cout << "i = " << i << " val1 = " << val1 << " val2 = " << val2
193 // << std::endl;
194 TEUCHOS_TEST_FOR_EXCEPTION(val1 < 0.0, std::logic_error,
195 "Stokhos::StieltjesBasis::stieltjes(): "
196 << " Polynomial " << i << " out of " << nfinish
197 << " has norm " << val1
198 << "! Try increasing number of quadrature points");
199 nrm[i] = val1;
200 a[i] = val2/val1;
201 b[i] = nrm[i]/nrm[i-1];
202 }
203}
204
205template <typename ordinal_type, typename value_type, typename func_type>
206void
208integrateBasisSquared(ordinal_type k,
209 const Teuchos::Array<value_type>& a,
210 const Teuchos::Array<value_type>& b,
211 const Teuchos::Array<value_type>& weights,
212 const Teuchos::Array<value_type>& points,
213 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals,
214 value_type& val1, value_type& val2) const
215{
216 evaluateRecurrence(k, a, b, points, phi_vals);
217 ordinal_type nqp = weights.size();
218 val1 = value_type(0);
219 val2 = value_type(0);
220 for (ordinal_type i=0; i<nqp; i++) {
221 val1 += weights[i]*phi_vals[i][k]*phi_vals[i][k];
222 val2 += weights[i]*phi_vals[i][k]*phi_vals[i][k]*points[i];
223 }
224}
225
226template <typename ordinal_type, typename value_type, typename func_type>
227void
229evaluateRecurrence(ordinal_type k,
230 const Teuchos::Array<value_type>& a,
231 const Teuchos::Array<value_type>& b,
232 const Teuchos::Array<value_type>& points,
233 Teuchos::Array< Teuchos::Array<value_type> >& values) const
234{
235 ordinal_type np = points.size();
236 if (k == 0)
237 for (ordinal_type i=0; i<np; i++)
238 values[i][k] = 1.0;
239 else if (k == 1)
240 for (ordinal_type i=0; i<np; i++)
241 values[i][k] = points[i] - a[k-1];
242 else
243 for (ordinal_type i=0; i<np; i++)
244 values[i][k] =
245 (points[i] - a[k-1])*values[i][k-1] - b[k-1]*values[i][k-2];
246}
247
248template <typename ordinal_type, typename value_type, typename func_type>
249Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
251cloneWithOrder(ordinal_type p) const
252{
253 return Teuchos::rcp(new Stokhos::StieltjesBasis<ordinal_type,value_type,func_type>(p,*this));
254}
255
256template <typename ordinal_type, typename value_type, typename func_type>
258StieltjesBasis(ordinal_type p, const StieltjesBasis& basis) :
259 RecurrenceBasis<ordinal_type, value_type>(p, basis),
260 func(basis.func),
261 quad(basis.quad),
262 pce_weights(quad->getQuadWeights()),
263 basis_values(quad->getBasisAtQuadPoints()),
264 func_vals(),
265 phi_vals(),
266 use_pce_quad_points(basis.use_pce_quad_points)
267{
268 this->setup();
269}
Abstract base class for quadrature methods.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship:
virtual void setup()
Setup basis after computing recurrence coefficients.
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.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a functional map...
virtual bool computeRecurrenceCoefficients(ordinal_type n, Teuchos::Array< value_type > &alpha, Teuchos::Array< value_type > &beta, Teuchos::Array< value_type > &delta, Teuchos::Array< value_type > &gamma) const
Compute recurrence coefficients.
StieltjesBasis(ordinal_type p, const Teuchos::RCP< const func_type > &func, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, bool use_pce_quad_points, bool normalize=false)
Constructor.
void stieltjes(ordinal_type nstart, ordinal_type nfinish, const Teuchos::Array< value_type > &weights, const Teuchos::Array< value_type > &points, Teuchos::Array< value_type > &a, Teuchos::Array< value_type > &b, Teuchos::Array< value_type > &nrm, Teuchos::Array< Teuchos::Array< value_type > > &phi_vals) const
Compute 3-term recurrence using Stieljtes procedure.
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
Get Gauss quadrature points, weights, and values of basis at points.
virtual void setup()
Setup basis after computing recurrence coefficients.
virtual Teuchos::RCP< OneDOrthogPolyBasis< ordinal_type, value_type > > cloneWithOrder(ordinal_type p) const
Clone this object with the option of building a higher order basis.
void integrateBasisSquared(ordinal_type k, const Teuchos::Array< value_type > &a, const Teuchos::Array< value_type > &b, const Teuchos::Array< value_type > &weights, const Teuchos::Array< value_type > &points, Teuchos::Array< Teuchos::Array< value_type > > &phi_vals, value_type &val1, value_type &val2) const
Compute and .
void evaluateRecurrence(ordinal_type k, const Teuchos::Array< value_type > &a, const Teuchos::Array< value_type > &b, const Teuchos::Array< value_type > &points, Teuchos::Array< Teuchos::Array< value_type > > &values) const
Evaluate polynomials via 3-term recurrence.