Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
gram_schmidt_example2.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"
49
50#include "Teuchos_CommandLineProcessor.hpp"
51
53
54int main(int argc, char **argv)
55{
56 try {
57
58 // Setup command line options
59 Teuchos::CommandLineProcessor CLP;
60 CLP.setDocString(
61 "This example runs a Gram-Schmidt-based dimension reduction example.\n");
62
63 int d = 2;
64 CLP.setOption("d", &d, "Stochastic dimension");
65
66 int p = 5;
67 CLP.setOption("p", &p, "Polynomial order");
68
69 CLP.parse( argc, argv );
70
71 // Create product basis
72 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
73 for (int i=0; i<d; i++)
74 bases[i] = Teuchos::rcp(new basis_type(p, true));
75 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
76 Teuchos::rcp(new Stokhos::CompletePolynomialBasis<int,double>(bases));
77
78 std::cout << "original basis size = " << basis->size() << std::endl;
79
80 // Create approximation
81 Stokhos::OrthogPolyApprox<int,double> x(basis), u(basis), v(basis),
82 w(basis), w2(basis);
83 for (int i=0; i<d; i++) {
84 x.term(i, 1) = 1.0;
85 }
86
87 // Tensor product quadrature
88 Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
89 Teuchos::rcp(new Stokhos::TensorProductQuadrature<int,double>(basis));
90 // Teuchos::RCP<const Stokhos::Quadrature<int,double> > quad =
91 // Teuchos::rcp(new Stokhos::SparseGridQuadrature<int,double>(basis,
92 // p+1));
93
94 std::cout << "original quadrature size = " << quad->size() << std::endl;
95
96 // Triple product tensor
97 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
98 basis->computeTripleProductTensor();
99
100 // Quadrature expansion
101 Stokhos::QuadOrthogPolyExpansion<int,double> quad_exp(basis, Cijk, quad);
102
103 // Compute PCE via quadrature expansion
104 quad_exp.sin(u,x);
105 quad_exp.exp(v,x);
106 quad_exp.times(w,v,u);
107
108 // Create new basis from u and v
109 Teuchos::Array< Stokhos::OrthogPolyApprox<int,double> > pces(2);
110 pces[0] = u;
111 pces[1] = v;
112 Teuchos::ParameterList params;
113 params.set("Reduced Basis Method", "Monomial Proj Gram-Schmidt");
115 Teuchos::RCP< Stokhos::ReducedPCEBasis<int,double> > gs_basis =
116 factory.createReducedBasis(p, pces, quad, Cijk);
117 Teuchos::RCP<const Stokhos::Quadrature<int,double> > gs_quad =
118 gs_basis->getReducedQuadrature();
119 Stokhos::OrthogPolyApprox<int,double> u_gs(gs_basis), v_gs(gs_basis),
120 w_gs(gs_basis);
121 gs_basis->transformFromOriginalBasis(u.coeff(), u_gs.coeff());
122 gs_basis->transformFromOriginalBasis(v.coeff(), v_gs.coeff());
123
124 std::cout << "reduced basis size = " << gs_basis->size() << std::endl;
125 std::cout << "reduced quadrature size = " << gs_quad->size() << std::endl;
126
127 Stokhos::OrthogPolyApprox<int,double> u2(basis), v2(basis);
128 gs_basis->transformToOriginalBasis(u_gs.coeff(), u2.coeff());
129 gs_basis->transformToOriginalBasis(v_gs.coeff(), v2.coeff());
130 double err_u = 0.0;
131 double err_v = 0.0;
132 for (int i=0; i<basis->size(); i++) {
133 double eu = std::abs(u[i]-u2[i]);
134 double ev = std::abs(v[i]-v2[i]);
135 if (eu > err_u) err_u = eu;
136 if (ev > err_v) err_v = ev;
137 }
138 std::cout << "error in u transformation = " << err_u << std::endl;
139 std::cout << "error in v transformation = " << err_v << std::endl;
140
141 // Triple product tensor
142 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > gs_Cijk =
143 Teuchos::null;
144
145 // Gram-Schmidt quadrature expansion
146 Teuchos::RCP< Teuchos::ParameterList > gs_exp_params =
147 Teuchos::rcp(new Teuchos::ParameterList);
148 gs_exp_params->set("Use Quadrature for Times", true);
150 gs_Cijk,
151 gs_quad,
152 gs_exp_params);
153
154 // Compute w_gs = u_gs*v_gs in Gram-Schmidt basis
155 gs_quad_exp.times(w_gs, u_gs, v_gs);
156
157 // Project w_gs back to original basis
158 gs_basis->transformToOriginalBasis(w_gs.coeff(), w2.coeff());
159
160 std::cout.precision(12);
161 std::cout << "w = " << std::endl << w;
162 std::cout << "w2 = " << std::endl << w2;
163 std::cout << "w_gs = " << std::endl << w_gs;
164
165 double err_w = 0.0;
166 for (int i=0; i<basis->size(); i++) {
167 double ew = std::abs(w[i]-w2[i]);
168 if (ew > err_w) err_w = ew;
169 }
170
171 std::cout.setf(std::ios::scientific);
172 std::cout << "w.mean() = " << w.mean() << std::endl
173 << "w2.mean() = " << w2.mean() << std::endl
174 << "mean error = "
175 << std::abs(w.mean()-w2.mean()) << std::endl
176 << "w.std_dev() = " << w.standard_deviation() << std::endl
177 << "w2.std_dev() = " << w2.standard_deviation() << std::endl
178 << "std_dev error = "
179 << std::abs(w.standard_deviation()-w2.standard_deviation())
180 << std::endl
181 << "w coeff error = " << err_w << std::endl;
182
183 }
184 catch (std::exception& e) {
185 std::cout << e.what() << std::endl;
186 }
187}
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
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.
pointer coeff()
Return coefficient array.
Orthogonal polynomial expansions based on numerical quadrature.
void exp(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void times(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a, const OrthogPolyApprox< ordinal_type, value_type, node_type > &b)
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Generate a basis from a given set of PCE expansions that is orthogonal with respect to the product me...
virtual Teuchos::RCP< Stokhos::ReducedPCEBasis< ordinal_type, value_type > > createReducedBasis(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
Get reduced quadrature object.
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