Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_StieltjesPCEBasisImp.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_Assert.hpp"
45#include "Teuchos_BLAS.hpp"
46#include "Teuchos_TimeMonitor.hpp"
47
48template <typename ordinal_type, typename value_type>
51 ordinal_type p,
52 const Teuchos::RCP<const Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce_,
53 const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad_,
54 bool use_pce_quad_points_,
55 bool normalize,
56 bool project_integrals_,
57 const Teuchos::RCP<const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk_) :
58 RecurrenceBasis<ordinal_type, value_type>("Stieltjes PCE", p, normalize),
59 pce(pce_),
60 quad(quad_),
61 pce_weights(quad->getQuadWeights()),
62 basis_values(quad->getBasisAtQuadPoints()),
63 pce_vals(),
64 phi_vals(),
65 use_pce_quad_points(use_pce_quad_points_),
66 fromStieltjesMat(p+1,pce->size()),
67 project_integrals(project_integrals_),
68 basis(pce->basis()),
69 Cijk(Cijk_),
70 phi_pce_coeffs()
71{
72 // Setup rest of recurrence basis
73 this->setup();
74}
75
76template <typename ordinal_type, typename value_type>
79{
80}
81
82template <typename ordinal_type, typename value_type>
83void
85getQuadPoints(ordinal_type quad_order,
86 Teuchos::Array<value_type>& quad_points,
87 Teuchos::Array<value_type>& quad_weights,
88 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
89{
90#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
91 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::StieltjesPCEBasis -- compute Gauss points");
92#endif
93
94 // Use underlying pce's quad points, weights, values
95 if (use_pce_quad_points) {
96 quad_points = pce_vals;
97 quad_weights = pce_weights;
98 quad_values = phi_vals;
99 return;
100 }
101
102 // Call base class
103 ordinal_type num_points =
104 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
105
106 // We can't reliably generate quadrature points of order > 2*p
107 //if (!project_integrals && quad_order > 2*this->p)
108 if (quad_order > 2*this->p)
109 quad_order = 2*this->p;
111 quad_points,
112 quad_weights,
113 quad_values);
114
115 // Fill in the rest of the points with zero weight
116 if (quad_weights.size() < num_points) {
117 ordinal_type old_size = quad_weights.size();
118 quad_weights.resize(num_points);
119 quad_points.resize(num_points);
120 quad_values.resize(num_points);
121 for (ordinal_type i=old_size; i<num_points; i++) {
122 quad_weights[i] = value_type(0);
123 quad_points[i] = quad_points[0];
124 quad_values[i].resize(this->p+1);
125 this->evaluateBases(quad_points[i], quad_values[i]);
126 }
127 }
128}
129
130template <typename ordinal_type, typename value_type>
131bool
133computeRecurrenceCoefficients(ordinal_type n,
134 Teuchos::Array<value_type>& alpha,
135 Teuchos::Array<value_type>& beta,
136 Teuchos::Array<value_type>& delta,
137 Teuchos::Array<value_type>& gamma) const
138{
139 ordinal_type nqp = phi_vals.size();
140 Teuchos::Array<value_type> nrm(n);
141 Teuchos::Array< Teuchos::Array<value_type> > vals(nqp);
142 for (ordinal_type i=0; i<nqp; i++)
143 vals[i].resize(n);
144 stieltjes(0, n, pce_weights, pce_vals, alpha, beta, nrm, vals);
145 for (ordinal_type i=0; i<n; i++) {
146 delta[i] = value_type(1.0);
147 gamma[i] = value_type(1.0);
148 }
149
150 // Save basis functions at quad point values
151 if (n == this->p+1)
152 phi_vals = vals;
153
154 return false;
155}
156
157template <typename ordinal_type, typename value_type>
158void
160setup()
161{
162 // Evaluate PCE at quad points
163 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
164 quad->getQuadPoints();
165 ordinal_type nqp = pce_weights.size();
166 pce_vals.resize(nqp);
167 phi_vals.resize(nqp);
168 for (ordinal_type i=0; i<nqp; i++) {
169 pce_vals[i] = pce->evaluate(quad_points[i], basis_values[i]);
170 phi_vals[i].resize(this->p+1);
171 }
172
173 if (project_integrals)
174 phi_pce_coeffs.resize(basis->size());
175
177
178 ordinal_type sz = pce->size();
179 fromStieltjesMat.putScalar(0.0);
180 for (ordinal_type i=0; i<=this->p; i++) {
181 for (ordinal_type j=0; j<sz; j++) {
182 for (ordinal_type k=0; k<nqp; k++)
183 fromStieltjesMat(i,j) +=
184 pce_weights[k]*phi_vals[k][i]*basis_values[k][j];
185 fromStieltjesMat(i,j) /= basis->norm_squared(j);
186 }
187 }
188}
189
190template <typename ordinal_type, typename value_type>
191void
193stieltjes(ordinal_type nstart,
194 ordinal_type nfinish,
195 const Teuchos::Array<value_type>& weights,
196 const Teuchos::Array<value_type>& points,
197 Teuchos::Array<value_type>& a,
198 Teuchos::Array<value_type>& b,
199 Teuchos::Array<value_type>& nrm,
200 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals) const
201{
202#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
203 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::StieltjesPCEBasis -- Discretized Stieltjes Procedure");
204#endif
205
206 value_type val1, val2;
207 ordinal_type start = nstart;
208 if (nstart == 0) {
209 if (project_integrals)
210 integrateBasisSquaredProj(0, a, b, weights, points, phi_vals, val1, val2);
211 else
212 integrateBasisSquared(0, a, b, weights, points, phi_vals, val1, val2);
213 nrm[0] = val1;
214 a[0] = val2/val1;
215 b[0] = value_type(1);
216 start = 1;
217 }
218 for (ordinal_type i=start; i<nfinish; i++) {
219 if (project_integrals)
220 integrateBasisSquaredProj(i, a, b, weights, points, phi_vals, val1, val2);
221 else
222 integrateBasisSquared(i, a, b, weights, points, phi_vals, val1, val2);
223 // std::cout << "i = " << i << " val1 = " << val1 << " val2 = " << val2
224 // << std::endl;
225 TEUCHOS_TEST_FOR_EXCEPTION(val1 < 0.0, std::logic_error,
226 "Stokhos::StieltjesPCEBasis::stieltjes(): "
227 << " Polynomial " << i << " out of " << nfinish
228 << " has norm " << val1
229 << "! Try increasing number of quadrature points");
230 nrm[i] = val1;
231 a[i] = val2/val1;
232 b[i] = nrm[i]/nrm[i-1];
233
234 // std::cout << "i = " << i << " alpha = " << a[i] << " beta = " << b[i]
235 // << " nrm = " << nrm[i] << std::endl;
236 }
237}
238
239template <typename ordinal_type, typename value_type>
240void
242integrateBasisSquared(ordinal_type k,
243 const Teuchos::Array<value_type>& a,
244 const Teuchos::Array<value_type>& b,
245 const Teuchos::Array<value_type>& weights,
246 const Teuchos::Array<value_type>& points,
247 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals,
248 value_type& val1, value_type& val2) const
249{
250 evaluateRecurrence(k, a, b, points, phi_vals);
251 ordinal_type nqp = weights.size();
252 val1 = value_type(0);
253 val2 = value_type(0);
254 for (ordinal_type i=0; i<nqp; i++) {
255 val1 += weights[i]*phi_vals[i][k]*phi_vals[i][k];
256 val2 += weights[i]*phi_vals[i][k]*phi_vals[i][k]*points[i];
257 }
258}
259
260template <typename ordinal_type, typename value_type>
261void
263evaluateRecurrence(ordinal_type k,
264 const Teuchos::Array<value_type>& a,
265 const Teuchos::Array<value_type>& b,
266 const Teuchos::Array<value_type>& points,
267 Teuchos::Array< Teuchos::Array<value_type> >& values) const
268{
269 ordinal_type np = points.size();
270 if (k == 0)
271 for (ordinal_type i=0; i<np; i++)
272 values[i][k] = 1.0;
273 else if (k == 1)
274 for (ordinal_type i=0; i<np; i++)
275 values[i][k] = points[i] - a[k-1];
276 else
277 for (ordinal_type i=0; i<np; i++)
278 values[i][k] =
279 (points[i] - a[k-1])*values[i][k-1] - b[k-1]*values[i][k-2];
280}
281
282template <typename ordinal_type, typename value_type>
283void
286 ordinal_type k,
287 const Teuchos::Array<value_type>& a,
288 const Teuchos::Array<value_type>& b,
289 const Teuchos::Array<value_type>& weights,
290 const Teuchos::Array<value_type>& points,
291 Teuchos::Array< Teuchos::Array<value_type> >& phi_vals,
292 value_type& val1, value_type& val2) const
293{
294 ordinal_type nqp = weights.size();
295 ordinal_type npc = basis->size();
296 const Teuchos::Array<value_type>& norms = basis->norm_squared();
297
298 // Compute PC expansion of phi_k in original basis
299 evaluateRecurrence(k, a, b, points, phi_vals);
300 for (ordinal_type j=0; j<npc; j++) {
301 value_type c = value_type(0);
302 for (ordinal_type i=0; i<nqp; i++)
303 c += weights[i]*phi_vals[i][k]*basis_values[i][j];
304 c /= norms[j];
305 phi_pce_coeffs[j] = c;
306 }
307
308 // Compute \int phi_k^2(\eta) d\eta
309 val1 = value_type(0);
310 for (ordinal_type j=0; j<npc; j++)
311 val1 += phi_pce_coeffs[j]*phi_pce_coeffs[j]*norms[j];
312
313 // Compute \int \eta phi_k^2(\eta) d\eta
314 val2 = value_type(0);
315 for (typename Cijk_type::k_iterator k_it = Cijk->k_begin();
316 k_it != Cijk->k_end(); ++k_it) {
317 ordinal_type l = index(k_it);
318 for (typename Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
319 j_it != Cijk->j_end(k_it); ++j_it) {
320 ordinal_type j = index(j_it);
321 for (typename Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
322 i_it != Cijk->i_end(j_it); ++i_it) {
323 ordinal_type i = index(i_it);
324 value_type c = value(i_it);
325 val2 += phi_pce_coeffs[i]*phi_pce_coeffs[j]*(*pce)[l]*c;
326 }
327 }
328 }
329}
330
331template <typename ordinal_type, typename value_type>
332void
334transformCoeffsFromStieltjes(const value_type *in, value_type *out) const
335{
336 Teuchos::BLAS<ordinal_type, value_type> blas;
337 blas.GEMV(Teuchos::TRANS, fromStieltjesMat.numRows(),
338 fromStieltjesMat.numCols(), 1.0, fromStieltjesMat.values(),
339 fromStieltjesMat.numRows(), in, 1, 0.0, out, 1);
340}
341
342template <typename ordinal_type, typename value_type>
343Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
345cloneWithOrder(ordinal_type p) const
346{
347 return Teuchos::rcp(new Stokhos::StieltjesPCEBasis<ordinal_type,value_type>(p,*this));
348}
349
350template <typename ordinal_type, typename value_type>
352StieltjesPCEBasis(ordinal_type p, const StieltjesPCEBasis& sbasis) :
353 RecurrenceBasis<ordinal_type, value_type>(p, sbasis),
354 pce(sbasis.pce),
355 quad(sbasis.quad),
356 pce_weights(quad->getQuadWeights()),
357 basis_values(quad->getBasisAtQuadPoints()),
358 pce_vals(sbasis.pce_vals),
359 phi_vals(),
360 use_pce_quad_points(sbasis.use_pce_quad_points),
361 fromStieltjesMat(p+1,pce->size()),
362 project_integrals(sbasis.project_integrals),
363 basis(pce->basis()),
364 Cijk(sbasis.Cijk),
365 phi_pce_coeffs()
366{
367 this->setup();
368}
Class to store coefficients of a projection onto an orthogonal polynomial basis.
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.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
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.
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.
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 .
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 integrateBasisSquaredProj(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 by projecting onto original PCE basis.
virtual void setup()
Setup basis after computing recurrence coefficients.
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.
StieltjesPCEBasis(ordinal_type p, const Teuchos::RCP< const Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, bool use_pce_quad_points, bool normalize=false, bool project_integrals=false, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk=Teuchos::null)
Constructor.
void transformCoeffsFromStieltjes(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
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.
Bi-directional iterator for traversing a sparse array.