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,
127 for (
ordinal_type i=0; i<term_prefix.dimension()-1; ++i)
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();
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();
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();
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);