42#include "Teuchos_UnitTestHarness.hpp"
43#include "Teuchos_TestingHelpers.hpp"
44#include "Teuchos_UnitTestRepository.hpp"
45#include "Teuchos_GlobalMPISession.hpp"
54 template <
typename OrdinalType,
typename ValueType>
73 template <
typename ordinal_type>
75 const Teuchos::Array<ordinal_type>& basis_orders,
78 Teuchos::FancyOStream& out) {
83 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(d);
88 Teuchos::RCP<const basis_type> basis = Teuchos::rcp(
new basis_type(bases));
92 Teuchos::RCP<node_type> head =
96 TEUCHOS_TEST_EQUALITY(head->block_size, basis->size(), out, success);
99 Teuchos::Array< Teuchos::RCP<node_type> > node_stack;
100 Teuchos::Array< ordinal_type > index_stack;
101 node_stack.push_back(head);
102 index_stack.push_back(0);
103 Teuchos::RCP<node_type> node;
106 while (node_stack.size() > 0) {
107 node = node_stack.back();
108 child_index = index_stack.back();
112 if (node->children.size() > 0) {
113 block_size_expected = 0;
115 block_size_expected += node->children[i]->block_size;
118 block_size_expected = node->terms.size();
120 TEUCHOS_TEST_EQUALITY(node->block_size, block_size_expected,
128 sum_prev += term_prefix[i];
130 ordinal_type my_p = std::min(p-sum_prev,basis_orders[d_prev]);
134 TEUCHOS_TEST_EQUALITY(node->block_size, block_size_expected2,
138 if (child_index < node->terms.size() && node->children.size() > 0)
139 out << node->terms[child_index] <<
" : block_size = "
140 << node->children[child_index]->block_size << std::endl;
143 if (node->children.size() == 0) {
144 TEUCHOS_TEST_EQUALITY(node->block_size, node->terms.size(),
151 out << term <<
" : index = " << index
152 <<
" index expected = " << index_expected << std::endl;
153 TEUCHOS_TEST_EQUALITY(index, index_expected, out, success);
155 node_stack.pop_back();
156 index_stack.pop_back();
160 else if (child_index < node->children.size()) {
161 ++index_stack.back();
162 node = node->children[child_index];
163 node_stack.push_back(node);
164 index_stack.push_back(0);
169 node_stack.pop_back();
170 index_stack.pop_back();
177 template <
typename ordinal_type,
typename value_type>
179 const Teuchos::Array<ordinal_type>& basis_orders,
181 const bool symmetric,
185 Teuchos::FancyOStream& out) {
191 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(d);
200 Teuchos::RCP<const basis_type> basis =
201 Teuchos::rcp(
new basis_type(bases, sparse_tol));
204 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
205 basis->computeTripleProductTensor();
209 Teuchos::RCP<Cijk_LTB_type> Cijk_LTB =
210 computeTripleProductTensorLTB(*basis, symmetric);
213 TEUCHOS_TEST_EQUALITY(Cijk->num_entries(), Cijk_LTB->num_entries(),
217 typedef typename Cijk_LTB_type::CijkNode
node_type;
220 Teuchos::Array< Teuchos::RCP<const node_type> > node_stack;
221 Teuchos::Array< ordinal_type > index_stack;
222 node_stack.push_back(Cijk_LTB->getHeadNode());
223 index_stack.push_back(0);
224 Teuchos::RCP<const node_type> node;
226 while (node_stack.size() > 0) {
227 node = node_stack.back();
228 child_index = index_stack.back();
232 TEUCHOS_TEST_EQUALITY(node->my_num_entries, node->values.size(),
234 Cijk_Iterator cijk_iterator(node->p_i, node->p_j, node->p_k, symmetric);
249 <<
"Check: rel_err( C(" << I <<
"," << J <<
"," << K <<
") )"
250 <<
" = " <<
"rel_err( " << cijk <<
", " << cijk2 <<
" ) = "
251 << err <<
" <= " << tol <<
" : ";
252 if (s) out <<
"Passed.";
257 success = success && s;
258 more = cijk_iterator.increment();
261 TEUCHOS_TEST_EQUALITY(node->my_num_entries, idx, out, success);
262 node_stack.pop_back();
263 index_stack.pop_back();
267 else if (child_index < node->children.size()) {
268 ++index_stack.back();
269 node = node->children[child_index];
270 node_stack.push_back(node);
271 index_stack.push_back(0);
276 node_stack.pop_back();
277 index_stack.pop_back();
285 template <
typename ordinal_type,
typename value_type>
287 const Teuchos::Array<ordinal_type>& basis_orders,
289 const bool symmetric,
293 Teuchos::FancyOStream& out) {
299 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(d);
308 Teuchos::RCP<const basis_type> basis =
309 Teuchos::rcp(
new basis_type(bases, sparse_tol));
312 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
313 basis->computeTripleProductTensor();
317 Teuchos::RCP<Cijk_LTB_type> Cijk_LTB =
318 computeTripleProductTensorLTBBlockLeaf(*basis, symmetric);
321 typedef typename Cijk_LTB_type::CijkNode
node_type;
323 Teuchos::Array< Teuchos::RCP<const node_type> > node_stack;
324 Teuchos::Array< ordinal_type > index_stack;
325 node_stack.push_back(Cijk_LTB->getHeadNode());
326 index_stack.push_back(0);
327 Teuchos::RCP<const node_type> node;
329 while (node_stack.size() > 0) {
330 node = node_stack.back();
331 child_index = index_stack.back();
335 TEUCHOS_TEST_EQUALITY(node->my_num_entries, node->values.size(),
347 if (symmetric && (k0%2 != (i+
j)%2)) ++k0;
355 if (J == K) cijk *= 2.0;
363 <<
"Check: rel_err( C(" << I <<
"," << J <<
","
365 <<
" = " <<
"rel_err( " << cijk <<
", " << cijk2 <<
" ) = "
366 << err <<
" <= " << tol <<
" : ";
367 if (s) out <<
"Passed.";
372 success = success && s;
375 value_type tol2 = atol + std::abs(cijk3)*rtol;
377 bool s2 = err2 < tol2;
380 <<
"Check: rel_err( C(" << I <<
"," << J <<
","
382 <<
" = " <<
"rel_err( " << cijk <<
", " << cijk3 <<
" ) = "
383 << err <<
" <= " << tol <<
" : ";
384 if (s2) out <<
"Passed.";
389 success = success && s2;
393 TEUCHOS_TEST_EQUALITY(node->my_num_entries, idx, out, success);
394 node_stack.pop_back();
395 index_stack.pop_back();
399 else if (child_index < node->children.size()) {
400 ++index_stack.back();
401 node = node->children[child_index];
402 node_stack.push_back(node);
403 index_stack.push_back(0);
408 node_stack.pop_back();
409 index_stack.pop_back();
417 template <
typename ordinal_type,
typename value_type>
419 const Teuchos::Array<ordinal_type>& basis_orders,
421 const bool symmetric,
425 Teuchos::FancyOStream& out) {
431 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(d);
440 Teuchos::RCP<const basis_type> basis =
441 Teuchos::rcp(
new basis_type(bases, sparse_tol));
444 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
445 basis->computeTripleProductTensor();
452 Teuchos::RCP<const Stokhos::Quadrature<ordinal_type,value_type> > quad =
458 Teuchos::RCP< Stokhos::FlatLTBSparse3Tensor<ordinal_type,value_type> >
464 x(basis), a(basis), b(basis), c1(basis), c2(basis);
475 Stokhos::flatLTB3TensorMultiply<10>(c2, a, b, *flat_Cijk);
484 Teuchos::Array<ordinal_type> basis_orders(
setup.d,
setup.p);
489 Teuchos::Array<ordinal_type> basis_orders(
setup.d);
491 basis_orders[i] = i+1;
496 Teuchos::Array<ordinal_type> basis_orders(
setup.d,
setup.p);
503 Teuchos::Array<ordinal_type> basis_orders(
setup.d);
505 basis_orders[i] = i+1;
512 Teuchos::Array<ordinal_type> basis_orders(
setup.d,
setup.p);
519 Teuchos::Array<ordinal_type> basis_orders(
setup.d);
521 basis_orders[i] = i+1;
528 Teuchos::Array<ordinal_type> basis_orders(
setup.d,
setup.p);
535 Teuchos::Array<ordinal_type> basis_orders(
setup.d,
setup.p);
542 Teuchos::Array<ordinal_type> basis_orders(
setup.d,
setup.p);
549 Teuchos::Array<ordinal_type> basis_orders(
setup.d);
551 basis_orders[i] = i+1;
558 Teuchos::Array<ordinal_type> basis_orders(
setup.d);
560 basis_orders[i] = i+1;
568 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
569 int res = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
int main(int argc, char *argv[])
Orthogonal polynomial expansions limited to algebraic operations.
Data structure storing a sparse 3-tensor C(i,j,k) in a a tree-based format for lexicographically orde...
Legendre polynomial basis.
A comparison functor implementing a strict weak ordering based lexographic ordering.
A multidimensional index.
Teuchos::Array< ordinal_type > index
index terms
ordinal_type dimension() const
Dimension.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
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)
Orthogonal polynomial expansions based on numerical quadrature.
void cos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Stokhos::LegendreBasis< int, double > basis_type
Teuchos::RCP< FlatLTBSparse3Tensor< ordinal_type, value_type > > computeFlatTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Teuchos::RCP< LexicographicTreeBasisNode< ordinal_type > > build_lexicographic_basis_tree(const Teuchos::ArrayView< const ordinal_type > &basis_orders, const ordinal_type total_order, const ordinal_type index_begin=ordinal_type(0), const ordinal_type order_sum=ordinal_type(0), const Stokhos::MultiIndex< ordinal_type > &term_prefix=Stokhos::MultiIndex< ordinal_type >())
ordinal_type n_choose_k(const ordinal_type &n, const ordinal_type &k)
Compute bionomial coefficient (n ; k) = n!/( k! (n-k)! )
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)
bool test_lexicographic_tree_sparse_3_tensor_block(const Teuchos::Array< ordinal_type > &basis_orders, const ordinal_type p, const bool symmetric, const value_type sparse_tol, const value_type atol, const value_type rtol, Teuchos::FancyOStream &out)
bool test_lexicographic_tree_sparse_3_tensor_multiply(const Teuchos::Array< ordinal_type > &basis_orders, const ordinal_type p, const bool symmetric, const value_type sparse_tol, const value_type atol, const value_type rtol, Teuchos::FancyOStream &out)
bool test_lexicographic_tree_sparse_3_tensor(const Teuchos::Array< ordinal_type > &basis_orders, const ordinal_type p, const bool symmetric, const value_type sparse_tol, const value_type atol, const value_type rtol, Teuchos::FancyOStream &out)
bool test_lexicographic_tree_coeffs(const Teuchos::Array< ordinal_type > &basis_orders, const ordinal_type p, const bool isotropic, Teuchos::FancyOStream &out)
UnitTestSetup< ordinal_type, value_type > setup
TEUCHOS_UNIT_TEST(LexicographicTreeCoefficients, Isotropic)