Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_LanczosPCEBasisImp.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_BLAS.hpp"
44#include "Teuchos_TimeMonitor.hpp"
45
46template <typename ordinal_type, typename value_type>
49 ordinal_type p,
50 const Teuchos::RCP< const Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce_,
51 const Teuchos::RCP< const Stokhos::Quadrature<ordinal_type, value_type> >& quad_,
52 bool normalize,
53 bool limit_integration_order_) :
54 RecurrenceBasis<ordinal_type, value_type>("Lanczos PCE", p, normalize),
55 pce(pce_),
56 quad(quad_),
57 limit_integration_order(limit_integration_order_),
58 nqp(quad->size()),
59 pce_weights(Teuchos::Copy,
60 const_cast<value_type*>(quad->getQuadWeights().getRawPtr()),
61 nqp),
62 pce_vals(nqp),
63 u0(nqp),
64 lanczos_vecs(nqp, p+1),
65 fromStieltjesMat(),
66 new_pce()
67{
68 // Evaluate PCE at quad points
69 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
70 quad->getQuadPoints();
71 const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
72 quad->getBasisAtQuadPoints();
73 for (ordinal_type i=0; i<nqp; i++) {
74 pce_vals[i] = pce->evaluate(quad_points[i], basis_values[i]);
75 u0[i] = value_type(1);
76 }
77
78 // Setup rest of basis
79 this->setup();
80}
81
82template <typename ordinal_type, typename value_type>
85{
86}
87
88template <typename ordinal_type, typename value_type>
89void
91getQuadPoints(ordinal_type quad_order,
92 Teuchos::Array<value_type>& quad_points,
93 Teuchos::Array<value_type>& quad_weights,
94 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
95{
96#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
97 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::LanczosPCEBasis -- compute Gauss points");
98#endif
99
100 // Call base class
101 ordinal_type num_points =
102 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
103
104 // We can't always reliably generate quadrature points of order > 2*p
105 // when using sparse grids for the underlying quadrature
106 if (limit_integration_order && quad_order > 2*this->p)
107 quad_order = 2*this->p;
109 quad_points,
110 quad_weights,
111 quad_values);
112
113 // Fill in the rest of the points with zero weight
114 if (quad_weights.size() < num_points) {
115 ordinal_type old_size = quad_weights.size();
116 quad_weights.resize(num_points);
117 quad_points.resize(num_points);
118 quad_values.resize(num_points);
119 for (ordinal_type i=old_size; i<num_points; i++) {
120 quad_weights[i] = value_type(0);
121 quad_points[i] = quad_points[0];
122 quad_values[i].resize(this->p+1);
123 this->evaluateBases(quad_points[i], quad_values[i]);
124 }
125 }
126}
127
128template <typename ordinal_type, typename value_type>
129Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
131cloneWithOrder(ordinal_type p) const
132{
134 p,*this));
135}
136
137template <typename ordinal_type, typename value_type>
138value_type
140getNewCoeffs(ordinal_type i) const
141{
142 return new_pce[i];
143}
144
145template <typename ordinal_type, typename value_type>
146void
148transformCoeffsFromLanczos(const value_type *in, value_type *out) const
149{
150 Teuchos::BLAS<ordinal_type, value_type> blas;
151 ordinal_type sz = fromStieltjesMat.numRows();
152 blas.GEMV(Teuchos::NO_TRANS, sz, this->p+1,
153 value_type(1.0), fromStieltjesMat.values(), sz,
154 in, ordinal_type(1), value_type(0.0), out, ordinal_type(1));
155}
156
157template <typename ordinal_type, typename value_type>
158bool
160computeRecurrenceCoefficients(ordinal_type n,
161 Teuchos::Array<value_type>& alpha,
162 Teuchos::Array<value_type>& beta,
163 Teuchos::Array<value_type>& delta,
164 Teuchos::Array<value_type>& gamma) const
165{
166 Teuchos::Array<value_type> nrm(n);
167 vectorspace_type vs(pce_weights);
168 operator_type A(pce_vals);
169
170 // Create space to store lanczos vectors -- use lanczos_vecs if
171 // we are requesting p+1 vectors
172 Teuchos::RCP<matrix_type> lv;
173 if (n == this->p+1)
174 lv = Teuchos::rcp(&lanczos_vecs, false);
175 else
176 lv = Teuchos::rcp(new matrix_type(nqp,n));
177
178 if (this->normalize)
179 lanczos_type::computeNormalized(n, vs, A, u0, *lv, alpha, beta, nrm);
180 else
181 lanczos_type::compute(n, vs, A, u0, *lv, alpha, beta, nrm);
182
183 for (ordinal_type i=0; i<n; i++) {
184 delta[i] = value_type(1.0);
185 }
186 if (this->normalize)
187 gamma = beta;
188 else
189 for (ordinal_type i=0; i<n; i++)
190 gamma[i] = value_type(1.0);
191
192 return this->normalize;
193}
194
195template <typename ordinal_type, typename value_type>
196void
198setup()
199{
201
202 // Compute transformation matrix back to original basis
203 ordinal_type sz = pce->size();
204 fromStieltjesMat.shape(sz, this->p+1);
205 fromStieltjesMat.putScalar(0.0);
206 const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
207 quad->getBasisAtQuadPoints();
208 for (ordinal_type i=0; i<sz; i++) {
209 for (ordinal_type j=0; j<=this->p; j++) {
210 for (ordinal_type k=0; k<nqp; k++)
211 fromStieltjesMat(i,j) +=
212 pce_weights[k]*lanczos_vecs(k,j)*basis_values[k][i];
213 fromStieltjesMat(i,j) /= pce->basis()->norm_squared(i);
214 }
215 }
216
217 // Project original PCE into the new basis
218 new_pce.resize(this->p+1);
219 vector_type u(sz);
220 for (ordinal_type i=0; i<sz; i++)
221 u[i] = (*pce)[i]*pce->basis()->norm_squared(i);
222 new_pce.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, fromStieltjesMat, u,
223 0.0);
224 for (ordinal_type i=0; i<=this->p; i++)
225 new_pce[i] /= this->norms[i];
226}
227
228template <typename ordinal_type, typename value_type>
230LanczosPCEBasis(ordinal_type p, const LanczosPCEBasis& basis) :
231 RecurrenceBasis<ordinal_type, value_type>(p, basis),
232 pce(basis.pce),
233 quad(basis.quad),
234 limit_integration_order(basis.limit_integration_order),
235 nqp(basis.nqp),
236 pce_weights(basis.pce_weights),
237 pce_vals(basis.pce_vals),
238 u0(basis.u0),
239 lanczos_vecs(nqp, p+1),
240 fromStieltjesMat(),
241 new_pce()
242{
243 this->setup();
244}
Copy
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
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.
lanczos_type::matrix_type matrix_type
Teuchos::RCP< const Stokhos::OrthogPolyApprox< ordinal_type, value_type > > pce
PCE Lanczos procedure is based on.
virtual void setup()
Setup basis after computing recurrence coefficients.
vector_type pce_vals
Values of PCE at quadrature points.
ordinal_type nqp
Number of quadrature points.
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.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > quad
Quadrature object.
void transformCoeffsFromLanczos(const value_type *in, value_type *out) const
Map expansion coefficients from this basis to original.
value_type getNewCoeffs(ordinal_type i) const
Get new coefficients in this new 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
Get Gauss quadrature points, weights, and values of basis at points.
vector_type u0
Initial Lanczos vector.
lanczos_type::vector_type vector_type
LanczosPCEBasis(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 normalize, bool limit_integration_order)
Constructor.
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.