Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_DivisionOperatorUnitTest.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#include "Stokhos.hpp"
55
57
58 // Common setup for unit tests
59 template <typename OrdinalType, typename ValueType>
61 ValueType rtol, atol;
62 ValueType crtol, catol;
63 OrdinalType sz;
64 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> > basis;
65 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > quad;
66 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk, Cijk_linear;
67 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType> > exp, exp_linear;
68 Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion<OrdinalType,ValueType> > qexp;
70 ValueType a;
71 Teuchos::RCP< Stokhos::DenseDirectDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > direct_division_strategy;
72
73
75 rtol = 1e-10;//4
76 atol = 1e-10;//5
77 crtol = 1e-12;
78 catol = 1e-12;
79 a = 3.1;
80 const OrdinalType d = 2;
81 const OrdinalType p = 7;
82
83 // Create product basis
84 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases(d);
85 for (OrdinalType i=0; i<d; i++)
86 bases[i] =
88 basis =
90
91 // Tensor product quadrature
92 quad =
94
95 // Triple product tensor
96 Cijk = basis->computeTripleProductTensor();
97 Cijk_linear = basis->computeLinearTripleProductTensor();
98
99 // Algebraic expansion
100 exp =
102 exp_linear =
104
105 // Quadrature expansion
106 qexp =
108 //Dense Direct Division Operator
111
112
113 // Create approximation
114// sz = basis->size();
115 x.reset(basis);
116 y.reset(basis);
117 u.reset(basis);
118 u2.reset(basis);
119 cx.reset(basis, 1);
120 x.term(0, 0) = 1.0;
121// y.term(0, 0) = 1.0:
122 cx.term(0, 0) = a;
123 cu.reset(basis);
124// cu2.reset(basis, 1);
125// sx.reset(basis, d+1);
126// su.reset(basis, d+1);
127// su2.reset(basis, d+1);
128 for (OrdinalType i=0; i<d; i++) {
129 x.term(i, 1) = 1.0;
130// y.term(i, 1) = 0.1;
131 }
132// y.term(0, 0) = 2.0;
133// for (OrdinalType i=0; i<d; i++)
134// y.term(i, 1) = 0.25;
135
136 c1.reset(basis);
137 c1.term(0,0)=1;
138 exp->exp(cu, x);
139
140 }
141};
142
144
145
146
147 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_Divide ) {
148 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_division_strategy =
149 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 0, 100, 0, 0,1));
150 cg_division_strategy->divide(setup.u, 1.0, setup.c1, setup.x, 0.0);
151 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.x, 0.0);
152 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
153 setup.rtol, setup.atol, out);
154 }
155 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_Jacobi_Divide ) {
156 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_diag_division_strategy =
157 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 2, 100, 0, 0,1));
158 cg_diag_division_strategy->divide(setup.u, 1.0, setup.c1, setup.x, 0.0);
159 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.x, 0.0);
160 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
161 setup.rtol, setup.atol, out);
162 }
163 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_SymGaussSeidel_Divide ) {
164 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_jacobi_division_strategy =
165 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 3, 100, 0, 0,1));
166 cg_jacobi_division_strategy->divide(setup.u, 1.0, setup.c1, setup.x, 0.0);
167 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.x, 0.0);
168 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
169 setup.rtol, setup.atol, out);
170 }
171 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_Schur_Divide ) {
172 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_schur_division_strategy =
173 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 0, 1e-12, 4, 100, 0, 0,1));
174 cg_schur_division_strategy->divide(setup.u, 1.0, setup.c1, setup.x, 0.0);
175 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.x, 0.0);
176 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
177 setup.rtol, setup.atol, out);
178 }
179
180 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_Nonlin_Divide ) {
181 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_nonlin_division_strategy =
182 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 0, 100, 0, 0,1));
183 cg_nonlin_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
184 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
185 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
186 setup.rtol, setup.atol, out);
187 }
188 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_Nonlin_Jacobi_Divide ) {
189 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_nonlin_division_strategy =
190 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 2, 100, 0, 0,1));
191 cg_nonlin_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
192 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
193 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
194 setup.rtol, setup.atol, out);
195 }
196 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_Nonlin_SymGaussSeidel_Divide ) {
197 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_nonlin_division_strategy =
198 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 3, 100, 0, 0,1));
199 cg_nonlin_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
200 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
201 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
202 setup.rtol, setup.atol, out);
203 }
204
205 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_Nonlin_Schur_Divide ) {
206 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_nonlin_division_strategy =
207 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 0, 1e-12, 4, 100, 0, 0,1));
208 cg_nonlin_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
209 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
210 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
211 setup.rtol, setup.atol, out);
212 }
213 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, CG_Nonlin_Schur_linearprec_Divide ) {
214 Teuchos::RCP< Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > cg_nonlin_division_strategy =
215 Teuchos::rcp(new Stokhos::CGDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 0, 1e-12, 4, 100, 1, 0,1));
216 cg_nonlin_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
217 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
218 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
219 setup.rtol, setup.atol, out);
220 }
221 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_Divide ) {
222 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
223 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 0, 100, 0, 0,1));
224 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.x, 0.0);
225 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.x, 0.0);
226 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
227 setup.rtol, setup.atol, out);
228 }
229
230 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_Jacobi_Divide ) {
231 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
232 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 2, 100, 0, 0,1));
233 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.x, 0.0);
234 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.x, 0.0);
235 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
236 setup.rtol, setup.atol, out);
237 }
238 TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_GaussSeidel_Divide ) {
239 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
240 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 3, 100, 0, 0,1));
241 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.x, 0.0);
242 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.x, 0.0);
243 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
244 setup.rtol, setup.atol, out);
245 }
246
247TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_Schur_Divide ) {
248 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
249 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 0, 1e-12, 4, 100, 0, 0,1));
250 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.x, 0.0);
251 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.x, 0.0);
252 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
253 setup.rtol, setup.atol, out);
254 }
255TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_Nonlin_Divide ) {
256 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
257 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 0, 100, 0, 0,1));
258 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
259 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
260 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
261 setup.rtol, setup.atol, out);
262 }
263TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_Nonlin_Jacobi_Divide ) {
264 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
265 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 2, 100, 0, 0,1));
266 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
267 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
268 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
269 setup.rtol, setup.atol, out);
270
271 }
272TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_Nonlin_GaussSeidel_Divide ) {
273 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
274 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 1, 1e-12, 3, 100, 0, 0,1));
275 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
276 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
277 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
278 setup.rtol, setup.atol, out);
279
280 }
281
282
283TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_Nonlin_Schur_Divide ) {
284 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
285 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 0, 1e-12, 4, 100, 0, 0,1));
286 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
287 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
288 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
289 setup.rtol, setup.atol, out);
290 }
291TEUCHOS_UNIT_TEST( Stokhos_DivisionOperator, GMRES_Nonlin_Schur_linearprec_Divide ) {
292 Teuchos::RCP< Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> > > gmres_division_strategy =
293 Teuchos::rcp(new Stokhos::GMRESDivisionExpansionStrategy<int,double,Stokhos::StandardStorage<int, double> >(setup.basis, setup.Cijk, 0, 1e-12, 4, 100, 1, 0,1));
294 gmres_division_strategy->divide(setup.u, 1.0, setup.c1, setup.cu, 0.0);
295 setup.direct_division_strategy->divide(setup.u2, 1.0, setup.c1, setup.cu, 0.0);
296 success = Stokhos::comparePCEs(setup.u, "u", setup.u2, "u2",
297 setup.rtol, setup.atol, out);
298 }
299
300
301
302
303
304}
305int main( int argc, char* argv[] ) {
306 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
307 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
308}
309
310
int main(int argc, char *argv[])
Strategy interface for computing PCE of a/b using only b[0].
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Strategy interface for computing PCE of a/b using only b[0].
Strategy interface for computing PCE of a/b using only b[0].
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.
reference term(ordinal_type dimension, ordinal_type order)
Get coefficient term for given dimension and order.
Orthogonal polynomial expansions based on numerical quadrature.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
TEUCHOS_UNIT_TEST(Stokhos_DivisionOperator, CG_Divide)
bool comparePCEs(const PCEType &a1, const std::string &a1_name, const Stokhos::OrthogPolyApprox< OrdinalType, ValueType > &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > c1
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp
Teuchos::RCP< Stokhos::Sparse3Tensor< int, double > > Cijk
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > x
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u2
Teuchos::RCP< Stokhos::DenseDirectDivisionExpansionStrategy< int, double, Stokhos::StandardStorage< int, double > > > direct_division_strategy
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > u
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > y
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > qexp
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > sx
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > quad
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cu2
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su2
Teuchos::RCP< Stokhos::Sparse3Tensor< int, double > > Cijk_linear
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > cx
Stokhos::OrthogPolyApprox< OrdinalType, ValueType > su
Teuchos::RCP< Stokhos::QuadOrthogPolyExpansion< OrdinalType, ValueType > > exp_linear