Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_Sparse3TensorUnitTest.cpp
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 "Teuchos_UnitTestHarness.hpp"
43#include "Teuchos_TestingHelpers.hpp"
44#include "Teuchos_UnitTestRepository.hpp"
45#include "Teuchos_GlobalMPISession.hpp"
46
47#include "Stokhos.hpp"
49
51
53
54 // Common setup for unit tests
55 template <typename OrdinalType, typename ValueType>
57 ValueType rtol, atol, sparse_tol;
58 OrdinalType d, p, sz;
59 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases;
60 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<OrdinalType,ValueType> > basis;
61 Teuchos::RCP<const Stokhos::Quadrature<OrdinalType,ValueType> > quad;
62 Teuchos::RCP<Cijk_type> Cijk;
63
64 UnitTestSetup(): d(3), p(4), bases(d) {
65 rtol = 1e-14;
66 atol = 1e-14;
67 sparse_tol = 1e-14;
68
69 // Create product basis
70 for (OrdinalType i=0; i<d; i++)
71 bases[i] =
73 basis =
75 sz = basis->size();
76
77 // Tensor product quadrature
78 quad =
80
81
82 // Triple product tensor
83 Cijk = basis->computeTripleProductTensor();
84 }
85
86 };
87
88 UnitTestSetup<int,double> setup;
89
90 TEUCHOS_UNIT_TEST( Stokhos_Sparse3Tensor, Values ) {
91 const Teuchos::Array<double>& weights = setup.quad->getQuadWeights();
92 const Teuchos::Array< Teuchos::Array<double> > & values =
93 setup.quad->getBasisAtQuadPoints();
94
95 success = true;
96 for (Cijk_type::k_iterator k_it=setup.Cijk->k_begin();
97 k_it!=setup.Cijk->k_end(); ++k_it) {
98 int k = Stokhos::index(k_it);
99 for (Cijk_type::kj_iterator j_it = setup.Cijk->j_begin(k_it);
100 j_it != setup.Cijk->j_end(k_it); ++j_it) {
101 int j = Stokhos::index(j_it);
102 for (Cijk_type::kji_iterator i_it = setup.Cijk->i_begin(j_it);
103 i_it != setup.Cijk->i_end(j_it); ++i_it) {
104 int i = Stokhos::index(i_it);
105 double c = Stokhos::value(i_it);
106
107 double c2 = 0.0;
108 int nqp = weights.size();
109 for (int qp=0; qp<nqp; qp++)
110 c2 += weights[qp]*values[qp][i]*values[qp][j]*values[qp][k];
111
112 double tol = setup.atol + c*setup.rtol;
113 double err = std::abs(c-c2);
114 bool s = err < tol;
115 if (!s) {
116 out << std::endl
117 << "Check: rel_err( C(" << i << "," << j << "," << k << ") )"
118 << " = " << "rel_err( " << c << ", " << c2 << " ) = " << err
119 << " <= " << tol << " : ";
120 if (s) out << "Passed.";
121 else
122 out << "Failed!";
123 out << std::endl;
124 }
125 success = success && s;
126 }
127 }
128 }
129 }
130
131 TEUCHOS_UNIT_TEST( Stokhos_Sparse3Tensor, Sparsity ) {
132 Teuchos::RCP< Stokhos::Sparse3Tensor<int, double> > Cijk_quad =
133 Teuchos::rcp(new Stokhos::Sparse3Tensor<int,double>);
134
135 // Create 1-D triple products
136 Teuchos::Array< Teuchos::RCP<Stokhos::Dense3Tensor<int,double> > > Cijk_1d(setup.d);
137 for (int i=0; i<setup.d; i++)
138 Cijk_1d[i] = setup.bases[i]->computeTripleProductTensor();
139
140 for (int j=0; j<setup.sz; j++) {
141 Stokhos::MultiIndex<int> terms_j = setup.basis->term(j);
142 for (int i=0; i<setup.sz; i++) {
143 Stokhos::MultiIndex<int> terms_i = setup.basis->term(i);
144 for (int k=0; k<setup.sz; k++) {
145 Stokhos::MultiIndex<int> terms_k = setup.basis->term(k);
146 double c = 1.0;
147 for (int l=0; l<setup.d; l++)
148 c *= (*Cijk_1d[l])(terms_i[l],terms_j[l],terms_k[l]);
149 if (std::abs(c/setup.basis->norm_squared(i)) > setup.sparse_tol)
150 Cijk_quad->add_term(i,j,k,c);
151 }
152 }
153 }
154 Cijk_quad->fillComplete();
155
156 // Check number of nonzeros
157 int nnz = setup.Cijk->num_entries();
158 int nnz_quad = Cijk_quad->num_entries();
159 success = (nnz == nnz_quad);
160 if (!success)
161 out << std::endl
162 << "Check: nnz(C) = " << nnz << " == nnz(C_quad) = " << nnz_quad
163 << ": Failed";
164 for (Cijk_type::k_iterator k_it=Cijk_quad->k_begin();
165 k_it!=Cijk_quad->k_end(); ++k_it) {
166 int k = Stokhos::index(k_it);
167 for (Cijk_type::kj_iterator j_it = Cijk_quad->j_begin(k_it);
168 j_it != Cijk_quad->j_end(k_it); ++j_it) {
169 int j = Stokhos::index(j_it);
170 for (Cijk_type::kji_iterator i_it = Cijk_quad->i_begin(j_it);
171 i_it != Cijk_quad->i_end(j_it); ++i_it) {
172 int i = Stokhos::index(i_it);
173 double c = setup.Cijk->getValue(i,j,k);
174 double c2 = Stokhos::value(i_it);
175
176 double tol = setup.atol + c2*setup.rtol;
177 double err = std::abs(c-c2);
178 bool s = err < tol;
179 if (!s) {
180 out << std::endl
181 << "Check: rel_err( C(" << i << "," << j << "," << k << ") )"
182 << " = " << "rel_err( " << c << ", " << c2 << " ) = " << err
183 << " <= " << tol << " : ";
184 if (s) out << "Passed.";
185 else
186 out << "Failed!";
187 out << std::endl;
188 }
189 success = success && s;
190 }
191 }
192 }
193 }
194
195 TEUCHOS_UNIT_TEST( Stokhos_Sparse3Tensor, GetValue ) {
196 success = true;
197 bool s;
198 double c, c_true;
199
200 // Check getValue() for a few different indices
201
202 c = setup.Cijk->getValue(0, 0, 0);
203 c_true = 1.0;
204 s = Stokhos::compareValues(c, "c", c_true, "c_true", setup.rtol, setup.atol,
205 out);
206 success = success && s;
207
208 c = setup.Cijk->getValue(9, 25, 4);
209 c_true = 0.04;
210 s = Stokhos::compareValues(c, "c", c_true, "c_true", setup.rtol, setup.atol,
211 out);
212 success = success && s;
213
214 c = setup.Cijk->getValue(8, 25, 4);
215 c_true = 0.0;
216 s = Stokhos::compareValues(c, "c", c_true, "c_true", setup.rtol, setup.atol,
217 out);
218 success = success && s;
219 }
220
221}
222
223int main( int argc, char* argv[] ) {
224 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
225 return Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
226}
int main(int argc, char *argv[])
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
A multidimensional index.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
j_sparse_array::const_iterator kji_iterator
Iterator for looping over i entries given k and j.
ji_sparse_array::const_iterator kj_iterator
Iterator for looping over j entries given k.
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
TEUCHOS_UNIT_TEST(Stokhos_Sparse3Tensor, Values)
Stokhos::Sparse3Tensor< int, double > Cijk_type
UnitTestSetup< int, double > setup
bool compareValues(const ValueType &a1, const std::string &a1_name, const ValueType &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
Teuchos::Array< Teuchos::RCP< const Stokhos::OneDOrthogPolyBasis< OrdinalType, ValueType > > > bases
Teuchos::RCP< const Stokhos::CompletePolynomialBasis< OrdinalType, ValueType > > basis
Teuchos::RCP< const Stokhos::Quadrature< OrdinalType, ValueType > > quad