Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
gram_schmidt_example.cpp
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 <iostream>
45#include <iomanip>
46
47#include "Stokhos.hpp"
51
53
58 pce(pce_), basis(basis_), vec(2) {}
59
60 double operator() (const double& a, const double& b) const {
61 vec[0] = a;
62 vec[1] = b;
63 return pce.evaluate(vec);
64 }
67 mutable Teuchos::Array<double> vec;
68};
69
70int main(int argc, char **argv)
71{
72 try {
73
74 const unsigned int d = 2;
75 const unsigned int p = 5;
76
77 // Create product basis
78 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
79 for (unsigned int i=0; i<d; i++)
80 bases[i] = Teuchos::rcp(new basis_type(p));
81 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
82 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
83
84 // Create approximation
85 Stokhos::OrthogPolyApprox<int,double> x(basis), u(basis), v(basis),
86 w(basis), w2(basis), w3(basis);
87 for (unsigned int i=0; i<d; i++) {
88 x.term(i, 1) = 1.0;
89 }
90
91 // Tensor product quadrature
92 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
93 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
94
95 // Triple product tensor
96 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
97 basis->computeTripleProductTensor();
98
99 // Quadrature expansion
100 Stokhos::QuadOrthogPolyExpansion<int,double> quad_exp(basis, Cijk, quad);
101
102 // Compute PCE via quadrature expansion
103 quad_exp.sin(u,x);
104 quad_exp.exp(v,x);
105 quad_exp.times(w,v,u);
106
107 // Compute tensor product Stieltjes basis for u and v
108 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > st_bases(2);
109 st_bases[0] =
111 p, Teuchos::rcp(&u,false), quad, true));
112 st_bases[1] =
114 p, Teuchos::rcp(&v,false), quad, true));
115 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > st_basis =
116 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(st_bases));
117 Stokhos::OrthogPolyApprox<int,double> u_st(st_basis), v_st(st_basis),
118 w_st(st_basis);
119 u_st.term(0, 0) = u.mean();
120 u_st.term(0, 1) = 1.0;
121 v_st.term(0, 0) = v.mean();
122 v_st.term(1, 1) = 1.0;
123
124 // Use Gram-Schmidt to orthogonalize Stieltjes basis of u and v
125 Teuchos::Array<double> st_points_0;
126 Teuchos::Array<double> st_weights_0;
127 Teuchos::Array< Teuchos::Array<double> > st_values_0;
128 st_bases[0]->getQuadPoints(p+1, st_points_0, st_weights_0, st_values_0);
129 Teuchos::Array<double> st_points_1;
130 Teuchos::Array<double> st_weights_1;
131 Teuchos::Array< Teuchos::Array<double> > st_values_1;
132 st_bases[1]->getQuadPoints(p+1, st_points_1, st_weights_1, st_values_1);
133 Teuchos::Array< Teuchos::Array<double> > st_points(st_points_0.size());
134 for (int i=0; i<st_points_0.size(); i++) {
135 st_points[i].resize(2);
136 st_points[i][0] = st_points_0[i];
137 st_points[i][1] = st_points_1[i];
138 }
139 Teuchos::Array<double> st_weights = st_weights_0;
140
141 Teuchos::RCP< Stokhos::GramSchmidtBasis<int,double> > gs_basis =
142 Teuchos::rcp(new Stokhos::GramSchmidtBasis<int,double>(st_basis,
143 st_points,
144 st_weights,
145 1e-15));
146
147 // Create quadrature for Gram-Schmidt basis using quad points and
148 // and weights from original basis mapped to Stieljtes basis
149 Teuchos::RCP< const Teuchos::Array< Teuchos::Array<double> > > points =
150 Teuchos::rcp(&st_points,false);
151 Teuchos::RCP< const Teuchos::Array<double> > weights =
152 Teuchos::rcp(&st_weights,false);
153 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gs_quad =
154 Teuchos::rcp(new Stokhos::UserDefinedQuadrature<int,double>(gs_basis,
155 points,
156 weights));
157
158 // Triple product tensor
159 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk =
160 gs_basis->computeTripleProductTensor();
161
162 // Gram-Schmidt quadrature expansion
164 gs_Cijk,
165 gs_quad);
166
167 Stokhos::OrthogPolyApprox<int,double> u_gs(gs_basis), v_gs(gs_basis),
168 w_gs(gs_basis);
169
170 // Map expansion in Stieltjes basis to Gram-Schmidt basis
171 gs_basis->transformCoeffs(u_st.coeff(), u_gs.coeff());
172 gs_basis->transformCoeffs(v_st.coeff(), v_gs.coeff());
173
174 // Compute w_gs = u_gs*v_gs in Gram-Schmidt basis
175 gs_quad_exp.times(w_gs, u_gs, v_gs);
176
177 // Project w_gs back to original basis
178 pce_quad_func gs_func(w_gs, *gs_basis);
179 quad_exp.binary_op(gs_func, w2, u, v);
180
181 // Triple product tensor
182 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
183 st_basis->computeTripleProductTensor();
184
185 // Stieltjes quadrature expansion
186 Teuchos::RCP<const Stokhos::Quadrature<int,double> > st_quad =
187 Teuchos::rcp(new Stokhos::UserDefinedQuadrature<int,double>(st_basis,
188 points,
189 weights));
191 st_Cijk,
192 st_quad);
193
194 // Compute w_st = u_st*v_st in Stieltjes basis
195 st_quad_exp.times(w_st, u_st, v_st);
196
197 // Project w_st back to original basis
198 pce_quad_func st_func(w_st, *st_basis);
199 quad_exp.binary_op(st_func, w3, u, v);
200
201 std::cout.precision(12);
202 std::cout << "w = " << std::endl << w;
203 std::cout << "w2 = " << std::endl << w2;
204 std::cout << "w3 = " << std::endl << w3;
205 std::cout << "w_gs = " << std::endl << w_gs;
206 std::cout << "w_st = " << std::endl << w_st;
207
208 std::cout.setf(std::ios::scientific);
209 std::cout << "w.mean() = " << w.mean() << std::endl
210 << "w2.mean() = " << w2.mean() << std::endl
211 << "w3.mean() = " << w3.mean() << std::endl
212 << "w_gs.mean() = " << w_gs.mean() << std::endl
213 << "w_st.mean() = " << w_st.mean() << std::endl
214 << "w.std_dev() = " << w.standard_deviation() << std::endl
215 << "w2.std_dev() = " << w2.standard_deviation() << std::endl
216 << "w3.std_dev() = " << w3.standard_deviation() << std::endl
217 << "w_gs.std_dev() = " << w_gs.standard_deviation() << std::endl
218 << "w_st.std_dev() = " << w_st.standard_deviation() << std::endl;
219 }
220 catch (std::exception& e) {
221 std::cout << e.what() << std::endl;
222 }
223}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Transforms a non-orthogonal multivariate basis to an orthogonal one using the Gram-Schmit procedure.
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
value_type mean() const
Compute mean of expansion.
value_type standard_deviation() const
Compute standard deviation of expansion.
Abstract base class for multivariate orthogonal polynomials.
Orthogonal polynomial expansions based on numerical quadrature.
Generates three-term recurrence using the Discretized Stieltjes procedure applied to a polynomial cha...
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
int main(int argc, char **argv)
Stokhos::LegendreBasis< int, double > basis_type
const Stokhos::OrthogPolyApprox< int, double > & pce
double operator()(const double &a, const double &b) const
const Stokhos::OrthogPolyBasis< int, double > & basis
pce_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
Teuchos::Array< double > vec