Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_LanczosUnitTest.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 "Teuchos_UnitTestHarness.hpp"
45#include "Teuchos_TestingHelpers.hpp"
46#include "Teuchos_UnitTestRepository.hpp"
47#include "Teuchos_GlobalMPISession.hpp"
48
49#include "Stokhos.hpp"
53
54// Quadrature functor to be passed into quadrature expansion for mapping
55// from Lanczos basis back to original PCE
59 pce(pce_), basis(basis_), vec(1) {}
60
61 double operator() (const double& a) const {
62 vec[0] = a;
63 return pce.evaluate(vec);
64 }
67 mutable Teuchos::Array<double> vec;
68};
69
70// Class encapsulating setup of the Lanczos basis for a given PCE
71// u = Func(x), where Func is specified by a template-parameter
72template <typename Func>
73struct Lanczos_PCE_Setup {
74 typedef typename Func::OrdinalType OrdinalType;
75 typedef typename Func::ValueType ValueType;
77 Func func;
79 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> > basis;
80 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType> > exp;
81 Teuchos::RCP<const Stokhos::LanczosProjPCEBasis<OrdinalType,ValueType> > st_1d_proj_basis;
82 Teuchos::RCP<const Stokhos::LanczosPCEBasis<OrdinalType,ValueType> > st_1d_basis;
83 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > st_bases;
84 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> > st_basis;
85 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > st_quad;
87
88 Lanczos_PCE_Setup(bool normalize, bool project) :
89 func()
90 {
91 rtol = 1e-8;
92 atol = 1e-10;
93 const OrdinalType d = 3;
94 const OrdinalType p = 5;
95
96 // Create product basis
97 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases(d);
98 for (OrdinalType i=0; i<d; i++)
99 bases[i] =
101 basis =
103
104 // Create approximation
105 sz = basis->size();
107 for (OrdinalType i=0; i<d; i++)
108 x.term(i, 1) = 1.0;
109
110 // Tensor product quadrature
111 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > quad =
113
114 // Triple product tensor
115 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk =
116 basis->computeTripleProductTensor();
117
118 // Quadrature expansion
120
121 // Compute PCE via quadrature expansion
122 u.reset(basis);
123 v.reset(basis);
124 func.eval(*exp, x, u);
125 exp->times(v,u,u);
126
127 // Compute Lanczos basis
128 st_bases.resize(1);
129 if (project) {
132 p, Teuchos::rcp(&u,false), Cijk, normalize));
134 }
135 else {
136 st_1d_basis =
138 p, Teuchos::rcp(&u,false), quad, normalize, false));
140 }
141
142 st_basis =
144 st_sz = st_basis->size();
147 if (project) {
148 u_st[0] = st_1d_proj_basis->getNewCoeffs(0);
149 u_st[1] = st_1d_proj_basis->getNewCoeffs(1);
150 }
151 else {
152 u_st[0] = st_1d_basis->getNewCoeffs(0);
153 u_st[1] = st_1d_basis->getNewCoeffs(1);
154 }
155
156 // Tensor product quadrature
157 st_quad =
159
160 // Triple product tensor
161 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > st_Cijk =
162 st_basis->computeTripleProductTensor();
163
164 // Quadrature expansion
166 st_Cijk,
167 st_quad);
168
169 st_exp.times(v_st, u_st, u_st);
170 }
171
172};
173
174#define LANCZOS_UNIT_TESTS(BASENAME, TAG, FUNC, NORMALIZE, PROJECT) \
175namespace BASENAME ## TAG { \
176 \
177 Lanczos_PCE_Setup< FUNC<int,double> > setup(NORMALIZE, PROJECT); \
178 \
179 TEUCHOS_UNIT_TEST( BASENAME, TAG ## Map ) { \
180 Stokhos::OrthogPolyApprox<int,double> u2(setup.basis); \
181 if (PROJECT) \
182 setup.st_1d_proj_basis->transformCoeffsFromLanczos( \
183 setup.u_st.coeff(), \
184 u2.coeff()); \
185 else \
186 setup.st_1d_basis->transformCoeffsFromLanczos( \
187 setup.u_st.coeff(), \
188 u2.coeff()); \
189 success = Stokhos::comparePCEs(setup.u, "u", u2, "u2", \
190 setup.rtol, setup.atol, out); \
191 } \
192 \
193 TEUCHOS_UNIT_TEST( BASENAME, TAG ## Orthog ) { \
194 const Teuchos::Array<double>& norms = \
195 setup.st_bases[0]->norm_squared(); \
196 const Teuchos::Array<double>& weights = \
197 setup.st_quad->getQuadWeights(); \
198 const Teuchos::Array< Teuchos::Array<double> >& values = \
199 setup.st_quad->getBasisAtQuadPoints(); \
200 Teuchos::SerialDenseMatrix<int,double> mat(setup.st_sz, \
201 setup.st_sz); \
202 for (int i=0; i<setup.st_sz; i++) { \
203 for (int j=0; j<setup.st_sz; j++) { \
204 for (unsigned int k=0; k<weights.size(); k++) \
205 mat(i,j) += weights[k]*values[k][i]*values[k][j]; \
206 mat(i,j) /= std::sqrt(norms[i]*norms[j]); \
207 } \
208 mat(i,i) -= 1.0; \
209 } \
210 success = mat.normInf() < setup.atol; \
211 if (!success) { \
212 out << "\n Error, mat.normInf() < atol = " << mat.normInf() \
213 << " < " << setup.atol << ": failed!\n"; \
214 out << "mat = " << printMat(mat) << std::endl; \
215 } \
216 } \
217 \
218 TEUCHOS_UNIT_TEST( BASENAME, TAG ## PCE ) { \
219 Stokhos::OrthogPolyApprox<int,double> v2(setup.basis); \
220 lanczos_pce_quad_func quad_func(setup.v_st, *setup.st_basis); \
221 setup.exp->unary_op(quad_func, v2, setup.u); \
222 success = comparePCEs(setup.v, "v", v2, "v2", setup.rtol, \
223 setup.atol, out); \
224 } \
225 \
226 TEUCHOS_UNIT_TEST( BASENAME, TAG ## Mean ) { \
227 success = Teuchos::testRelErr( \
228 "v.mean()", setup.v.mean(), \
229 "v_st.mean()", setup.v_st.mean(), \
230 "rtol", setup.rtol, \
231 "rtol", setup.rtol, \
232 Teuchos::Ptr<std::ostream>(out.getOStream().get())); \
233 \
234 } \
235 \
236 TEUCHOS_UNIT_TEST( BASENAME, TAG ## StandardDeviation ) { \
237 success = Teuchos::testRelErr( \
238 "v.standard_deviation()", \
239 setup.v.standard_deviation(), \
240 "v_st.standard_devaition()", \
241 setup.v_st.standard_deviation(), \
242 "rtol", 1e-1, \
243 "rtol", 1e-1, \
244 Teuchos::Ptr<std::ostream>(out.getOStream().get())); \
245 } \
246 \
247}
248
249//
250// Lanczos tests based on expansion of u = cos(x) where x is a U([-1,1])
251// random variable
252//
253template <typename Ordinal_Type, typename Value_Type>
254struct Lanczos_Cos_Func {
255 typedef Ordinal_Type OrdinalType;
256 typedef Value_Type ValueType;
257 static const bool is_even = true;
258 void
262 exp.cos(u,x);
263 }
264};
265LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_Proj, Cos, Lanczos_Cos_Func,
266 false, true)
267LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_ProjNorm, Cos, Lanczos_Cos_Func,
269LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis, Cos, Lanczos_Cos_Func,
270 false, false)
271LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_Norm, Cos, Lanczos_Cos_Func,
272 true, false)
273
274//
275// Lanczos tests based on expansion of u = sin(x) where x is a U([-1,1])
276// random variable
277//
278template <typename Ordinal_Type, typename Value_Type>
279struct Lanczos_Sin_Func {
280 typedef Ordinal_Type OrdinalType;
281 typedef Value_Type ValueType;
282 static const bool is_even = false;
283 void
287 exp.sin(u,x);
288 }
289};
290LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_Proj, Sin, Lanczos_Sin_Func,
291 false, true)
292LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_ProjNorm, Sin, Lanczos_Sin_Func,
293 true, true)
294LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis, Sin, Lanczos_Sin_Func,
295 false, false)
296LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_Norm, Sin, Lanczos_Sin_Func,
297 true, false)
298
299//
300// Lanczos tests based on expansion of u = exp(x) where x is a U([-1,1])
301// random variable. For this test we don't use the PCE quad points and
302// instead use those generated for the Lanczos basis
303//
304template <typename Ordinal_Type, typename Value_Type>
305struct Lanczos_Exp_Func {
306 typedef Ordinal_Type OrdinalType;
307 typedef Value_Type ValueType;
308 static const bool is_even = false;
309 void
313 exp.exp(u,x);
314 }
315};
316LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_Proj, Exp, Lanczos_Exp_Func,
317 false, true)
318LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_ProjNorm, Exp, Lanczos_Exp_Func,
319 true, true)
320LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis, Exp, Lanczos_Exp_Func,
321 false, false)
322LANCZOS_UNIT_TESTS(Stokhos_LanczosPCEBasis_Norm, Exp, Lanczos_Exp_Func,
323 true, false)
324
325int main( int argc, char* argv[] ) {
326 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
327 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
328}
#define LANCZOS_UNIT_TESTS(BASENAME, TAG, FUNC, NORMALIZE, PROJECT)
true false int main(int argc, char *argv[])
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...
Legendre polynomial basis.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void reset(const Teuchos::RCP< const Stokhos::OrthogPolyBasis< ordinal_type, value_type > > &new_basis, ordinal_type sz=0)
Reset to a new basis.
value_type evaluate(const Teuchos::Array< value_type > &point) const
Evaluate polynomial approximation at a point.
Abstract base class for multivariate orthogonal polynomials.
Orthogonal polynomial expansions based on numerical quadrature.
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)
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > st_quad
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v
Teuchos::RCP< const Stokhos::LanczosPCEBasis< OrdinalType, ValueType > > st_1d_basis
Teuchos::RCP< const Stokhos::HouseTriDiagPCEBasis< OrdinalType, ValueType > > st_1d_proj_basis
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > st_basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > v_st
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp
Teuchos::RCP< const Stokhos::LanczosProjPCEBasis< OrdinalType, ValueType > > st_1d_proj_basis
Func::OrdinalType OrdinalType
Lanczos_PCE_Setup(bool normalize, bool project)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u_st
Teuchos::Array< Teuchos::RCP< const Stokhos::OneDOrthogPolyBasis< OrdinalType, ValueType > > > st_bases
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
void eval(Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > &exp, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &x, Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &u)
double operator()(const double &a) const
const Stokhos::OrthogPolyBasis< int, double > & basis
lanczos_pce_quad_func(const Stokhos::OrthogPolyApprox< int, double > &pce_, const Stokhos::OrthogPolyBasis< int, double > &basis_)
const Stokhos::OrthogPolyApprox< int, double > & pce