Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_ProductLanczosGramSchmidtPCEBasisImp.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 "Stokhos_SDMUtils.hpp"
49
50template <typename ordinal_type, typename value_type>
53 ordinal_type max_p,
54 const Teuchos::Array< Stokhos::OrthogPolyApprox<ordinal_type, value_type> >& pce,
55 const Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >& quad,
56 const Teuchos::RCP< const Stokhos::Sparse3Tensor<ordinal_type, value_type> >& Cijk,
57 const Teuchos::ParameterList& params_) :
58 name("Product Lanczos Gram-Schmidt PCE Basis"),
59 params(params_),
60 p(max_p)
61{
62 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > pce_basis = pce[0].basis();
63 pce_sz = pce_basis->size();
64
65 // Check if basis is a product basis
66 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(pce_basis);
67 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type,value_type> > > coord_bases;
68 if (prod_basis != Teuchos::null)
69 coord_bases = prod_basis->getCoordinateBases();
70
71 // Build Lanczos basis for each pce
72 bool project = params.get("Project", true);
73 bool normalize = params.get("Normalize", true);
74 bool limit_integration_order = params.get("Limit Integration Order", false);
75 bool use_stieltjes = params.get("Use Old Stieltjes Method", false);
76 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double > > > coordinate_bases;
77 Teuchos::Array<int> is_invariant(pce.size(),-2);
78 for (ordinal_type i=0; i<pce.size(); i++) {
79
80 // Check for pce's lying in invariant subspaces, which are pce's that
81 // depend on only a single dimension. In this case use the corresponding
82 // original coordinate basis. Convention is: -2 -- not invariant, -1 --
83 // constant, i >= 0 pce depends only on dimension i
84 if (prod_basis != Teuchos::null)
85 is_invariant[i] = isInvariant(pce[i]);
86 if (is_invariant[i] >= 0) {
87 coordinate_bases.push_back(coord_bases[is_invariant[i]]);
88 }
89
90 // Exclude constant pce's from the basis since they don't represent
91 // stochastic dimensions
92 else if (is_invariant[i] != -1) {
93 if (use_stieltjes) {
94 coordinate_bases.push_back(
95 Teuchos::rcp(
97 p, Teuchos::rcp(&(pce[i]),false), quad, false,
98 normalize, project, Cijk)));
99 }
100 else {
101 if (project)
102 coordinate_bases.push_back(
103 Teuchos::rcp(
105 p, Teuchos::rcp(&(pce[i]),false), Cijk,
106 normalize, limit_integration_order)));
107 else
108 coordinate_bases.push_back(
109 Teuchos::rcp(
111 p, Teuchos::rcp(&(pce[i]),false), quad,
112 normalize, limit_integration_order)));
113 }
114 }
115 }
116 d = coordinate_bases.size();
117
118 // Build tensor product basis
120 Teuchos::rcp(
122 coordinate_bases,
123 params.get("Cijk Drop Tolerance", 1.0e-15),
124 params.get("Use Old Cijk Algorithm", false)));
125
126 // Build Psi matrix -- Psi_ij = Psi_i(x^j)*w_j/<Psi_i^2>
127 const Teuchos::Array<value_type>& weights = quad->getQuadWeights();
128 const Teuchos::Array< Teuchos::Array<value_type> >& points =
129 quad->getQuadPoints();
130 const Teuchos::Array< Teuchos::Array<value_type> >& basis_vals =
131 quad->getBasisAtQuadPoints();
132 ordinal_type nqp = weights.size();
133 SDM Psi(pce_sz, nqp);
134 for (ordinal_type i=0; i<pce_sz; i++)
135 for (ordinal_type k=0; k<nqp; k++)
136 Psi(i,k) = basis_vals[k][i]*weights[k]/pce_basis->norm_squared(i);
137
138 // Build Phi matrix -- Phi_ij = Phi_i(y(x^j))
139 sz = tensor_lanczos_basis->size();
140 Teuchos::Array<value_type> red_basis_vals(sz);
141 Teuchos::Array<value_type> pce_vals(d);
142 SDM Phi(sz, nqp);
143 SDM F(nqp, d);
144 for (int k=0; k<nqp; k++) {
145 ordinal_type jdx = 0;
146 for (int j=0; j<pce.size(); j++) {
147
148 // Exclude constant pce's
149 if (is_invariant[j] != -1) {
150
151 // Use the identity mapping for invariant subspaces
152 if (is_invariant[j] >= 0)
153 pce_vals[jdx] = points[k][is_invariant[j]];
154 else
155 pce_vals[jdx] = pce[j].evaluate(points[k], basis_vals[k]);
156 F(k,jdx) = pce_vals[jdx];
157 jdx++;
158
159 }
160
161 }
162 tensor_lanczos_basis->evaluateBases(pce_vals, red_basis_vals);
163 for (int i=0; i<sz; i++)
164 Phi(i,k) = red_basis_vals[i];
165 }
166
167 bool verbose = params.get("Verbose", false);
168
169 // Compute matrix A mapping reduced space to original
170 SDM A(pce_sz, sz);
171 ordinal_type ret =
172 A.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, Psi, Phi, 0.0);
173 TEUCHOS_ASSERT(ret == 0);
174
175 // Rescale columns of A to have unit norm
176 const Teuchos::Array<value_type>& basis_norms = pce_basis->norm_squared();
177 for (ordinal_type j=0; j<sz; j++) {
178 value_type nrm = 0.0;
179 for (ordinal_type i=0; i<pce_sz; i++)
180 nrm += A(i,j)*A(i,j)*basis_norms[i];
181 nrm = std::sqrt(nrm);
182 for (ordinal_type i=0; i<pce_sz; i++)
183 A(i,j) /= nrm;
184 }
185
186 // Compute our new basis -- each column of Qp is the coefficients of the
187 // new basis in the original basis. Constraint pivoting so first d+1
188 // columns and included in Qp.
189 value_type rank_threshold = params.get("Rank Threshold", 1.0e-12);
190 std::string orthogonalization_method =
191 params.get("Orthogonalization Method", "Householder");
192 Teuchos::Array<value_type> w(pce_sz, 1.0);
193 SDM R;
194 Teuchos::Array<ordinal_type> piv(sz);
195 for (int i=0; i<d+1; i++)
196 piv[i] = 1;
198 sz = SOF::createOrthogonalBasis(
199 orthogonalization_method, rank_threshold, verbose, A, w, Qp, R, piv);
200
201 // Original basis at quadrature points -- needed to transform expansions
202 // in this basis back to original
203 SDM B(nqp, pce_sz);
204 for (ordinal_type i=0; i<nqp; i++)
205 for (ordinal_type j=0; j<pce_sz; j++)
206 B(i,j) = basis_vals[i][j];
207
208 // Evaluate new basis at original quadrature points
209 Q.reshape(nqp, sz);
210 ret = Q.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, B, Qp, 0.0);
211 TEUCHOS_ASSERT(ret == 0);
212
213 // Compute reduced quadrature rule
215 params.sublist("Reduced Quadrature"));
216 Teuchos::SerialDenseMatrix<ordinal_type, value_type> Q2;
217 reduced_quad = quad_factory.createReducedQuadrature(Q, Q2, F, weights);
218
219 // Basis is orthonormal by construction
220 norms.resize(sz, 1.0);
221}
222
223template <typename ordinal_type, typename value_type>
226{
227}
228
229template <typename ordinal_type, typename value_type>
230ordinal_type
232order() const
233{
234 return p;
235}
236
237template <typename ordinal_type, typename value_type>
238ordinal_type
240dimension() const
241{
242 return d;
243}
244
245template <typename ordinal_type, typename value_type>
246ordinal_type
248size() const
249{
250 return sz;
251}
252
253template <typename ordinal_type, typename value_type>
254const Teuchos::Array<value_type>&
256norm_squared() const
257{
258 return norms;
259}
260
261template <typename ordinal_type, typename value_type>
262const value_type&
264norm_squared(ordinal_type i) const
265{
266 return norms[i];
267}
268
269template <typename ordinal_type, typename value_type>
270Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
273
274{
275 return Teuchos::null;
276}
277
278template <typename ordinal_type, typename value_type>
279Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
282
283{
284 return Teuchos::null;
285}
286
287template <typename ordinal_type, typename value_type>
288value_type
290evaluateZero(ordinal_type i) const
291{
292 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
293}
294
295template <typename ordinal_type, typename value_type>
296void
298evaluateBases(const Teuchos::ArrayView<const value_type>& point,
299 Teuchos::Array<value_type>& basis_vals) const
300{
301 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Not implemented!");
302}
303
304template <typename ordinal_type, typename value_type>
305void
307print(std::ostream& os) const
308{
309 os << "Product Lanczos Gram-Schmidt basis of order " << p << ", dimension "
310 << d
311 << ", and size " << sz << ". Matrix coefficients:\n";
312 os << printMat(Qp) << std::endl;
313 os << "Basis vector norms (squared):\n\t";
314 for (ordinal_type i=0; i<sz; i++)
315 os << norms[i] << " ";
316 os << "\n";
317}
318
319template <typename ordinal_type, typename value_type>
320const std::string&
322getName() const
323{
324 return name;
325}
326
327template <typename ordinal_type, typename value_type>
328void
330transformToOriginalBasis(const value_type *in, value_type *out,
331 ordinal_type ncol, bool transpose) const
332{
333 if (transpose) {
334 SDM zbar(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, sz);
335 SDM z(Teuchos::View, out, ncol, ncol, pce_sz);
336 ordinal_type ret =
337 z.multiply(Teuchos::NO_TRANS, Teuchos::TRANS, 1.0, zbar, Qp, 0.0);
338 TEUCHOS_ASSERT(ret == 0);
339 }
340 else {
341 SDM zbar(Teuchos::View, const_cast<value_type*>(in), sz, sz, ncol);
342 SDM z(Teuchos::View, out, pce_sz, pce_sz, ncol);
343 ordinal_type ret =
344 z.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, Qp, zbar, 0.0);
345 TEUCHOS_ASSERT(ret == 0);
346 }
347}
348
349template <typename ordinal_type, typename value_type>
350void
352transformFromOriginalBasis(const value_type *in, value_type *out,
353 ordinal_type ncol, bool transpose) const
354{
355 if (transpose) {
356 SDM z(Teuchos::View, const_cast<value_type*>(in), ncol, ncol, pce_sz);
357 SDM zbar(Teuchos::View, out, ncol, ncol, sz);
358 ordinal_type ret =
359 zbar.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, 1.0, z, Qp, 0.0);
360 TEUCHOS_ASSERT(ret == 0);
361 }
362 else {
363 SDM z(Teuchos::View, const_cast<value_type*>(in), pce_sz, pce_sz, ncol);
364 SDM zbar(Teuchos::View, out, sz, sz, ncol);
365 ordinal_type ret =
366 zbar.multiply(Teuchos::TRANS, Teuchos::NO_TRANS, 1.0, Qp, z, 0.0);
367 TEUCHOS_ASSERT(ret == 0);
368 }
369}
370
371template <typename ordinal_type, typename value_type>
372Teuchos::RCP<const Stokhos::Quadrature<ordinal_type, value_type> >
375{
376 return reduced_quad;
377}
378
379template <typename ordinal_type, typename value_type>
380ordinal_type
383{
384 Teuchos::RCP<const Stokhos::OrthogPolyBasis<ordinal_type,value_type> > basis =
385 pce.basis();
386 ordinal_type dim = basis->dimension();
387 value_type tol = 1.0e-15;
388
389 // Check if basis is a product basis
390 Teuchos::RCP<const Stokhos::ProductBasis<ordinal_type,value_type> > prod_basis = Teuchos::rcp_dynamic_cast<const Stokhos::ProductBasis<ordinal_type,value_type> >(basis);
391 if (prod_basis == Teuchos::null)
392 return -2;
393
394 // Build list of dimensions pce depends on by looping over each dimension,
395 // computing norm of pce with just that dimension -- note we don't include
396 // the constant term
397 Teuchos::Array<ordinal_type> dependent_dims;
398 tmp_pce.reset(basis);
399 for (ordinal_type i=0; i<dim; i++) {
400 ordinal_type p = prod_basis->getCoordinateBases()[i]->order();
401 tmp_pce.init(0.0);
402 for (ordinal_type j=1; j<=p; j++)
403 tmp_pce.term(i,j) = pce.term(i,j);
404 value_type nrm = tmp_pce.two_norm();
405 if (nrm > tol) dependent_dims.push_back(i);
406 }
407
408 // If dependent_dims has length 1, pce a function of a single variable,
409 // which is an invariant subspace
410 if (dependent_dims.size() == 1)
411 return dependent_dims[0];
412
413 // If dependent_dims has length 0, pce is constant
414 else if (dependent_dims.size() == 0)
415 return -1;
416
417 // Otherwise pce depends on more than one variable
418 return -2;
419}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Generates three-term recurrence using the Lanczos procedure applied to a polynomial chaos expansion i...
Class to store coefficients of a projection onto an orthogonal polynomial basis.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > basis() const
Return basis.
Encapsulate various orthogonalization (ie QR) methods.
SDM Q
Values of transformed basis at quadrature points.
ordinal_type isInvariant(const Stokhos::OrthogPolyApprox< ordinal_type, value_type > &pce) const
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
Teuchos::SerialDenseMatrix< ordinal_type, value_type > SDM
virtual const std::string & getName() const
Return string name of basis.
virtual void transformFromOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients from original basis to this basis.
virtual ordinal_type size() const
Return total size of basis.
virtual void transformToOriginalBasis(const value_type *in, value_type *out, ordinal_type ncol=1, bool transpose=false) const
Transform coefficients to original basis from this basis.
Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > reduced_quad
Reduced quadrature object.
virtual Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > getReducedQuadrature() const
Get reduced quadrature object.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
Teuchos::RCP< Stokhos::CompletePolynomialBasis< ordinal_type, value_type > > tensor_lanczos_basis
Product Lanczos basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
virtual void print(std::ostream &os) const
Print basis to stream os.
SDM Qp
Coefficients of transformed basis in original basis.
ProductLanczosGramSchmidtPCEBasis(ordinal_type p, const Teuchos::Array< Stokhos::OrthogPolyApprox< ordinal_type, value_type > > &pce, const Teuchos::RCP< const Stokhos::Quadrature< ordinal_type, value_type > > &quad, const Teuchos::RCP< const Stokhos::Sparse3Tensor< ordinal_type, value_type > > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())
Constructor.
Abstract base class for quadrature methods.
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual Teuchos::RCP< const Stokhos::UserDefinedQuadrature< ordinal_type, value_type > > createReducedQuadrature(const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &Q2, const Teuchos::SerialDenseMatrix< ordinal_type, value_type > &F, const Teuchos::Array< value_type > &weights) const
Get reduced quadrature object.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
void printMat(const char *name, Epetra_IntSerialDenseMatrix &matrix)