Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_MonoProjPCEBasisImp.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_LAPACK.hpp"
45#include "Teuchos_TimeMonitor.hpp"
46#include "Teuchos_SerialDenseHelpers.hpp"
47
48template <typename ordinal_type, typename value_type>
51 ordinal_type p,
55 bool limit_integration_order_) :
56 RecurrenceBasis<ordinal_type, value_type>("Monomial Projection", p, true),
57 limit_integration_order(limit_integration_order_),
58 pce_sz(pce.basis()->size()),
59 pce_norms(pce.basis()->norm_squared()),
60 a(pce_sz),
61 b(pce_sz),
62 basis_vecs(pce_sz, p+1),
63 new_pce(p+1)
64{
65 // If the original basis is normalized, we can use the standard QR
66 // factorization. For simplicity, we renormalize the PCE coefficients
67 // for a normalized basis
69 for (ordinal_type i=0; i<pce_sz; i++) {
70 pce_norms[i] = std::sqrt(pce_norms[i]);
71 normalized_pce[i] *= pce_norms[i];
72 }
73
74 // Evaluate PCE at quad points
75 ordinal_type nqp = quad.size();
76 Teuchos::Array<value_type> pce_vals(nqp);
77 const Teuchos::Array<value_type>& weights = quad.getQuadWeights();
78 const Teuchos::Array< Teuchos::Array<value_type> >& quad_points =
79 quad.getQuadPoints();
80 const Teuchos::Array< Teuchos::Array<value_type> >& basis_values =
82 for (ordinal_type i=0; i<nqp; i++) {
83 pce_vals[i] = normalized_pce.evaluate(quad_points[i], basis_values[i]);
84 }
85
86 // Form Kylov matrix up to order pce_sz
88
89 // Compute matrix
92 for (typename Cijk_type::k_iterator k_it = Cijk.k_begin();
93 k_it != Cijk.k_end(); ++k_it) {
94 ordinal_type k = index(k_it);
95 for (typename Cijk_type::kj_iterator j_it = Cijk.j_begin(k_it);
96 j_it != Cijk.j_end(k_it); ++j_it) {
97 ordinal_type j = index(j_it);
98 value_type val = 0;
99 for (typename Cijk_type::kji_iterator i_it = Cijk.i_begin(j_it);
100 i_it != Cijk.i_end(j_it); ++i_it) {
101 ordinal_type i = index(i_it);
102 value_type c = value(i_it) / (pce_norms[j]*pce_norms[k]);
103 val += pce[i]*c;
104 }
105 A(k,j) = val;
106 }
107 }
108
109 // Each column i is given by projection of the i-th order monomial
110 // onto original basis
111 vector_type u0 = Teuchos::getCol(Teuchos::View, K, 0);
112 u0(0) = 1.0;
113 for (ordinal_type i=1; i<pce_sz; i++)
114 u0(i) = 0.0;
115 for (ordinal_type k=1; k<pce_sz; k++) {
116 vector_type u = Teuchos::getCol(Teuchos::View, K, k);
117 vector_type up = Teuchos::getCol(Teuchos::View, K, k-1);
118 u.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, A, up, 0.0);
119 }
120 /*
121 for (ordinal_type j=0; j<pce_sz; j++) {
122 for (ordinal_type i=0; i<pce_sz; i++) {
123 value_type val = 0.0;
124 for (ordinal_type k=0; k<nqp; k++)
125 val += weights[k]*std::pow(pce_vals[k],j)*basis_values[k][i];
126 K(i,j) = val;
127 }
128 }
129 */
130
131 std::cout << K << std::endl << std::endl;
132
133 // Compute QR factorization of K
134 ordinal_type ws_size, info;
135 value_type ws_size_query;
136 Teuchos::Array<value_type> tau(pce_sz);
138 lapack.GEQRF(pce_sz, pce_sz, K.values(), K.stride(), &tau[0],
139 &ws_size_query, -1, &info);
140 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
141 "GEQRF returned value " << info);
142 ws_size = static_cast<ordinal_type>(ws_size_query);
143 Teuchos::Array<value_type> work(ws_size);
144 lapack.GEQRF(pce_sz, pce_sz, K.values(), K.stride(), &tau[0],
145 &work[0], ws_size, &info);
146 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
147 "GEQRF returned value " << info);
148
149 // Get Q
150 lapack.ORGQR(pce_sz, pce_sz, pce_sz, K.values(), K.stride(), &tau[0],
151 &ws_size_query, -1, &info);
152 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
153 "ORGQR returned value " << info);
154 ws_size = static_cast<ordinal_type>(ws_size_query);
155 work.resize(ws_size);
156 lapack.ORGQR(pce_sz, pce_sz, pce_sz, K.values(), K.stride(), &tau[0],
157 &work[0], ws_size, &info);
158 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
159 "ORGQR returned value " << info);
160
161 // Get basis vectors
162 for (ordinal_type j=0; j<p+1; j++)
163 for (ordinal_type i=0; i<pce_sz; i++)
164 basis_vecs(i,j) = K(i,j);
165
166
167 // Compute T = Q'*A*Q
169 prod.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, K, A, 0.0);
171 T.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, prod, K, 0.0);
172
173 //std::cout << T << std::endl;
174
175 // Recursion coefficients are diagonal and super diagonal
176 b[0] = 1.0;
177 for (ordinal_type i=0; i<pce_sz-1; i++) {
178 a[i] = T(i,i);
179 b[i+1] = T(i,i+1);
180 }
181 a[pce_sz-1] = T(pce_sz-1,pce_sz-1);
182
183 // Setup rest of basis
184 this->setup();
185
186 // Project original PCE into the new basis
188 for (ordinal_type i=0; i<pce_sz; i++)
189 u[i] = normalized_pce[i];
190 new_pce.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, basis_vecs, u,
191 0.0);
192 for (ordinal_type i=0; i<=p; i++)
193 new_pce[i] /= this->norms[i];
194}
195
196template <typename ordinal_type, typename value_type>
199{
200}
201
202template <typename ordinal_type, typename value_type>
203void
205getQuadPoints(ordinal_type quad_order,
206 Teuchos::Array<value_type>& quad_points,
207 Teuchos::Array<value_type>& quad_weights,
208 Teuchos::Array< Teuchos::Array<value_type> >& quad_values) const
209{
210#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
211 TEUCHOS_FUNC_TIME_MONITOR("Stokhos::MonoProjPCEBasis -- compute Gauss points");
212#endif
213
214 // Call base class
215 ordinal_type num_points =
216 static_cast<ordinal_type>(std::ceil((quad_order+1)/2.0));
217
218 // We can't always reliably generate quadrature points of order > 2*p
219 // when using sparse grids for the underlying quadrature
220 if (limit_integration_order && quad_order > 2*this->p)
221 quad_order = 2*this->p;
223 quad_points,
224 quad_weights,
225 quad_values);
226
227 // Fill in the rest of the points with zero weight
228 if (quad_weights.size() < num_points) {
229 ordinal_type old_size = quad_weights.size();
230 quad_weights.resize(num_points);
231 quad_points.resize(num_points);
232 quad_values.resize(num_points);
233 for (ordinal_type i=old_size; i<num_points; i++) {
234 quad_weights[i] = value_type(0);
235 quad_points[i] = quad_points[0];
236 quad_values[i].resize(this->p+1);
237 evaluateBases(quad_points[i], quad_values[i]);
238 }
239 }
240}
241
242template <typename ordinal_type, typename value_type>
243Teuchos::RCP<Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> >
245cloneWithOrder(ordinal_type p) const
246{
248 p,*this));
249}
250
251template <typename ordinal_type, typename value_type>
252value_type
254getNewCoeffs(ordinal_type i) const
255{
256 return new_pce[i];
257}
258
259template <typename ordinal_type, typename value_type>
260void
262transformCoeffs(const value_type *in, value_type *out) const
263{
264 // Transform coefficients to normalized basis
265 Teuchos::BLAS<ordinal_type, value_type> blas;
266 blas.GEMV(Teuchos::NO_TRANS, pce_sz, this->p+1,
267 value_type(1.0), basis_vecs.values(), pce_sz,
268 in, ordinal_type(1), value_type(0.0), out, ordinal_type(1));
269
270 // Transform from normalized to original
271 for (ordinal_type i=0; i<pce_sz; i++)
272 out[i] /= pce_norms[i];
273}
274
275template <typename ordinal_type, typename value_type>
276bool
278computeRecurrenceCoefficients(ordinal_type n,
279 Teuchos::Array<value_type>& alpha,
280 Teuchos::Array<value_type>& beta,
281 Teuchos::Array<value_type>& delta,
282 Teuchos::Array<value_type>& gamma) const
283{
284 // Get recurrence coefficients from the full set we already computed
285 for (ordinal_type i=0; i<n; i++) {
286 alpha[i] = a[i];
287 beta[i] = b[i];
288 delta[i] = value_type(1.0);
289 gamma[i] = b[i];
290
291 std::cout << "i = " << i << " alpha = " << alpha[i] << " beta = " << beta[i]
292 << " gamma = " << gamma[i] << std::endl;
293 }
294
295 return true;
296}
297
298template <typename ordinal_type, typename value_type>
300MonoProjPCEBasis(ordinal_type p, const MonoProjPCEBasis& basis) :
301 RecurrenceBasis<ordinal_type, value_type>("Lanczos PCE", p, false),
302 limit_integration_order(basis.limit_integration_order),
303 pce_sz(basis.pce_sz),
304 pce_norms(basis.pce_norms),
305 a(basis.a),
306 b(basis.b),
307 basis_vecs(basis.basis_vecs),
308 new_pce(basis.new_pce)
309{
310 this->setup();
311}
expr val()
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
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.
Teuchos::SerialDenseVector< ordinal_type, value_type > vector_type
vector_type new_pce
Projection of pce in new basis.
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.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > matrix_type
Teuchos::Array< value_type > a
Stores full set of alpha coefficients.
matrix_type basis_vecs
Basis vectors.
MonoProjPCEBasis(ordinal_type p, const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce, const Stokhos::Quadrature< ordinal_type, value_type > &quad, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, bool limit_integration_order=false)
Constructor.
Teuchos::Array< value_type > pce_norms
Basis norms.
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::Array< value_type > b
Stores full set of beta coefficients.
ordinal_type pce_sz
Size of PC expansion.
void transformCoeffs(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.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Abstract base class for quadrature methods.
virtual ordinal_type size() const =0
Get number of quadrature points.
virtual const Teuchos::Array< value_type > & getQuadWeights() const =0
Get quadrature weights.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getBasisAtQuadPoints() const =0
Get values of basis at quadrature points.
virtual const Teuchos::Array< Teuchos::Array< value_type > > & getQuadPoints() const =0
Get quadrature points.
Implementation of OneDOrthogPolyBasis based on the general three-term recurrence relationship:
Teuchos::Array< value_type > norms
Norms.
ordinal_type p
Order of basis.
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.
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
k_iterator k_begin() const
Iterator pointing to first k entry.
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
k_iterator k_end() const
Iterator pointing to last k entry.
Specialization for Sacado::UQ::PCE< Storage<...> >