Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_LexicographicTreeBasisUnitTest.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"
48 #include "Stokhos_UnitTestHelpers.hpp"
49 #include "Stokhos_LTBSparse3Tensor.hpp"
50 
51 namespace TotalOrderBasisUnitTest {
52 
53  // Common setup for unit tests
54  template <typename OrdinalType, typename ValueType>
55  struct UnitTestSetup {
56  ValueType rtol, atol, sparse_tol;
57  OrdinalType p,d;
58 
59  UnitTestSetup() {
60  rtol = 1e-12;
61  atol = 1e-12;
62  sparse_tol = 1e-12;
63  d = 3;
64  p = 5;
65  }
66 
67  };
68 
69  typedef int ordinal_type;
70  typedef double value_type;
71  UnitTestSetup<ordinal_type,value_type> setup;
72 
73  template <typename ordinal_type>
74  bool test_lexicographic_tree_coeffs(
75  const Teuchos::Array<ordinal_type>& basis_orders,
76  const ordinal_type p,
77  const bool isotropic,
78  Teuchos::FancyOStream& out) {
79  bool success = true;
80 
81  // Build total order basis of dimension d
82  ordinal_type d = basis_orders.size();
83  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(d);
84  for (ordinal_type i=0; i<d; i++)
85  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(basis_orders[i], true));
86  typedef Stokhos::LexographicLess< Stokhos::MultiIndex<int> > less_type;
87  typedef Stokhos::TotalOrderBasis<ordinal_type,value_type,less_type> basis_type;
88  Teuchos::RCP<const basis_type> basis = Teuchos::rcp(new basis_type(bases));
89 
90  // Build tree basis
91  typedef Stokhos::LexicographicTreeBasisNode<ordinal_type> node_type;
92  Teuchos::RCP<node_type> head =
93  Stokhos::build_lexicographic_basis_tree(basis_orders(), p);
94 
95  // Check head block size is equal to the whole basis size
96  TEUCHOS_TEST_EQUALITY(head->block_size, basis->size(), out, success);
97 
98  // Check tree is consistent
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;
104  ordinal_type child_index;
105  ordinal_type block_size_expected;
106  while (node_stack.size() > 0) {
107  node = node_stack.back();
108  child_index = index_stack.back();
109 
110  // Check block size is the sum of each child's block size
111  // or the number of terms for leaves
112  if (node->children.size() > 0) {
113  block_size_expected = 0;
114  for (ordinal_type i=0; i<node->children.size(); ++i)
115  block_size_expected += node->children[i]->block_size;
116  }
117  else {
118  block_size_expected = node->terms.size();
119  }
120  TEUCHOS_TEST_EQUALITY(node->block_size, block_size_expected,
121  out, success);
122 
123  // Check block size based on formula (only for isotropic)
124  if (isotropic) {
125  ordinal_type sum_prev = 0;
126  Stokhos::MultiIndex<ordinal_type> term_prefix = node->terms[0];
127  for (ordinal_type i=0; i<term_prefix.dimension()-1; ++i)
128  sum_prev += term_prefix[i];
129  ordinal_type d_prev = term_prefix.dimension()-1;
130  ordinal_type my_p = std::min(p-sum_prev,basis_orders[d_prev]);
131  ordinal_type d_left = d - d_prev;
132  ordinal_type block_size_expected2 =
133  Stokhos::n_choose_k(my_p+d_left,d_left);
134  TEUCHOS_TEST_EQUALITY(node->block_size, block_size_expected2,
135  out, success);
136  }
137 
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;
141 
142  // Leaf -- check global indices
143  if (node->children.size() == 0) {
144  TEUCHOS_TEST_EQUALITY(node->block_size, node->terms.size(),
145  out, success);
146  for (ordinal_type i=0; i<node->terms.size(); ++i) {
147  Stokhos::MultiIndex<ordinal_type> term = node->terms[i];
148  ordinal_type index = node->index_begin + i;
149  ordinal_type index_expected = basis->index(term);
150 
151  out << term << " : index = " << index
152  << " index expected = " << index_expected << std::endl;
153  TEUCHOS_TEST_EQUALITY(index, index_expected, out, success);
154  }
155  node_stack.pop_back();
156  index_stack.pop_back();
157  }
158 
159  // More children to process -- process them first
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);
165  }
166 
167  // No more children
168  else {
169  node_stack.pop_back();
170  index_stack.pop_back();
171  }
172 
173  }
174  return success;
175  }
176 
177  template <typename ordinal_type, typename value_type>
178  bool test_lexicographic_tree_sparse_3_tensor(
179  const Teuchos::Array<ordinal_type>& basis_orders,
180  const ordinal_type p,
181  const bool symmetric,
182  const value_type sparse_tol,
183  const value_type atol,
184  const value_type rtol,
185  Teuchos::FancyOStream& out) {
186 
187  bool success = true;
188  ordinal_type d = basis_orders.size();
189 
190  // Build tensor product basis
191  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(d);
192  if (symmetric)
193  for (ordinal_type i=0; i<d; i++)
194  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(basis_orders[i], true));
195  else
196  for (ordinal_type i=0; i<d; i++)
197  bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<ordinal_type,value_type>(basis_orders[i], 1.0, 2.0, true));
198  typedef Stokhos::LexographicLess< Stokhos::MultiIndex<int> > less_type;
199  typedef Stokhos::TotalOrderBasis<ordinal_type,value_type,less_type> basis_type;
200  Teuchos::RCP<const basis_type> basis =
201  Teuchos::rcp(new basis_type(bases, sparse_tol));
202 
203  // Build "standard" Cijk
204  Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
205  basis->computeTripleProductTensor();
206 
207  // Build LTB Cijk
208  typedef Stokhos::LTBSparse3Tensor<ordinal_type, value_type> Cijk_LTB_type;
209  Teuchos::RCP<Cijk_LTB_type> Cijk_LTB =
210  computeTripleProductTensorLTB(*basis, symmetric);
211 
212  // Check sizes
213  TEUCHOS_TEST_EQUALITY(Cijk->num_entries(), Cijk_LTB->num_entries(),
214  out, success);
215 
216  // Check non-zero structure
217  typedef typename Cijk_LTB_type::CijkNode node_type;
218  typedef Stokhos::ProductBasisUtils::Cijk_1D_Iterator<ordinal_type> Cijk_Iterator;
219 
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;
225  ordinal_type child_index;
226  while (node_stack.size() > 0) {
227  node = node_stack.back();
228  child_index = index_stack.back();
229 
230  // Leaf -- check Cijk indices and values
231  if (node->is_leaf) {
232  TEUCHOS_TEST_EQUALITY(node->my_num_entries, node->values.size(),
233  out, success);
234  Cijk_Iterator cijk_iterator(node->p_i, node->p_j, node->p_k, symmetric);
235  ordinal_type idx = 0;
236  bool more = true;
237  while (more) {
238  value_type cijk = node->values[idx];
239  ordinal_type I = node->i_begin + cijk_iterator.i;
240  ordinal_type J = node->j_begin + cijk_iterator.j;
241  ordinal_type K = node->k_begin + cijk_iterator.k;
242  value_type cijk2 = Cijk->getValue(I,J,K);
243 
244  value_type tol = atol + std::abs(cijk2)*rtol;
245  value_type err = std::abs(cijk-cijk2);
246  bool s = err < tol;
247  if (!s) {
248  out << std::endl
249  << "Check: rel_err( C(" << I << "," << J << "," << K << ") )"
250  << " = " << "rel_err( " << cijk << ", " << cijk2 << " ) = "
251  << err << " <= " << tol << " : ";
252  if (s) out << "Passed.";
253  else
254  out << "Failed!";
255  out << std::endl;
256  }
257  success = success && s;
258  more = cijk_iterator.increment();
259  ++idx;
260  }
261  TEUCHOS_TEST_EQUALITY(node->my_num_entries, idx, out, success);
262  node_stack.pop_back();
263  index_stack.pop_back();
264  }
265 
266  // More children to process -- process them first
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);
272  }
273 
274  // No more children
275  else {
276  node_stack.pop_back();
277  index_stack.pop_back();
278  }
279 
280  }
281 
282  return success;
283  }
284 
285  template <typename ordinal_type, typename value_type>
286  bool test_lexicographic_tree_sparse_3_tensor_block(
287  const Teuchos::Array<ordinal_type>& basis_orders,
288  const ordinal_type p,
289  const bool symmetric,
290  const value_type sparse_tol,
291  const value_type atol,
292  const value_type rtol,
293  Teuchos::FancyOStream& out) {
294 
295  bool success = true;
296  ordinal_type d = basis_orders.size();
297 
298  // Build tensor product basis
299  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(d);
300  if (symmetric)
301  for (ordinal_type i=0; i<d; i++)
302  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(basis_orders[i], true));
303  else
304  for (ordinal_type i=0; i<d; i++)
305  bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<ordinal_type,value_type>(basis_orders[i], 1.0, 2.0, true));
306  typedef Stokhos::LexographicLess< Stokhos::MultiIndex<int> > less_type;
307  typedef Stokhos::TotalOrderBasis<ordinal_type,value_type,less_type> basis_type;
308  Teuchos::RCP<const basis_type> basis =
309  Teuchos::rcp(new basis_type(bases, sparse_tol));
310 
311  // Build "standard" Cijk
312  Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
313  basis->computeTripleProductTensor();
314 
315  // Build LTB Cijk
316  typedef Stokhos::LTBSparse3Tensor<ordinal_type, value_type> Cijk_LTB_type;
317  Teuchos::RCP<Cijk_LTB_type> Cijk_LTB =
318  computeTripleProductTensorLTBBlockLeaf(*basis, symmetric);
319 
320  // Check non-zero structure
321  typedef typename Cijk_LTB_type::CijkNode node_type;
322 
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;
328  ordinal_type child_index;
329  while (node_stack.size() > 0) {
330  node = node_stack.back();
331  child_index = index_stack.back();
332 
333  // Leaf -- check Cijk indices and values
334  if (node->is_leaf) {
335  TEUCHOS_TEST_EQUALITY(node->my_num_entries, node->values.size(),
336  out, success);
337  ordinal_type idx = 0;
338  for (ordinal_type i=0; i<=node->p_i; ++i) {
339  for (ordinal_type j=0; j<=node->p_j; ++j) {
340  // ordinal_type k0 = node->parent_j_equals_k ?
341  // std::max(j,std::abs(i-j)) : std::abs(i-j);
342  // if (symmetric && (k0%2 != (i+j)%2)) ++k0;
343  // const ordinal_type k_end = std::min(node->p_k,i+j);
344  // const ordinal_type k_inc = symmetric ? 2 : 1;
345 
346  ordinal_type k0 = node->parent_j_equals_k ? j : 0;
347  if (symmetric && (k0%2 != (i+j)%2)) ++k0;
348  const ordinal_type k_end = node->p_k;
349  const ordinal_type k_inc = symmetric ? 2 : 1;
350  for (ordinal_type k=k0; k<=k_end; k+=k_inc) {
351  value_type cijk = node->values[idx++];
352  ordinal_type I = node->i_begin + i;
353  ordinal_type J = node->j_begin + j;
354  ordinal_type K = node->k_begin + k;
355  if (J == K) cijk *= 2.0;
356  value_type cijk2 = Cijk->getValue(I,J,K);
357 
358  value_type tol = atol + std::abs(cijk2)*rtol;
359  value_type err = std::abs(cijk-cijk2);
360  bool s = err < tol;
361  if (!s) {
362  out << std::endl
363  << "Check: rel_err( C(" << I << "," << J << ","
364  << K << ") )"
365  << " = " << "rel_err( " << cijk << ", " << cijk2 << " ) = "
366  << err << " <= " << tol << " : ";
367  if (s) out << "Passed.";
368  else
369  out << "Failed!";
370  out << std::endl;
371  }
372  success = success && s;
373 
374  value_type cijk3 = Cijk->getValue(I,K,J);
375  value_type tol2 = atol + std::abs(cijk3)*rtol;
376  value_type err2 = std::abs(cijk-cijk3);
377  bool s2 = err2 < tol2;
378  if (!s2) {
379  out << std::endl
380  << "Check: rel_err( C(" << I << "," << J << ","
381  << K << ") )"
382  << " = " << "rel_err( " << cijk << ", " << cijk3 << " ) = "
383  << err << " <= " << tol << " : ";
384  if (s2) out << "Passed.";
385  else
386  out << "Failed!";
387  out << std::endl;
388  }
389  success = success && s2;
390  }
391  }
392  }
393  TEUCHOS_TEST_EQUALITY(node->my_num_entries, idx, out, success);
394  node_stack.pop_back();
395  index_stack.pop_back();
396  }
397 
398  // More children to process -- process them first
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);
404  }
405 
406  // No more children
407  else {
408  node_stack.pop_back();
409  index_stack.pop_back();
410  }
411 
412  }
413 
414  return success;
415  }
416 
417  template <typename ordinal_type, typename value_type>
418  bool test_lexicographic_tree_sparse_3_tensor_multiply(
419  const Teuchos::Array<ordinal_type>& basis_orders,
420  const ordinal_type p,
421  const bool symmetric,
422  const value_type sparse_tol,
423  const value_type atol,
424  const value_type rtol,
425  Teuchos::FancyOStream& out) {
426 
427  bool success = true;
428  ordinal_type d = basis_orders.size();
429 
430  // Build tensor product basis
431  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(d);
432  if (symmetric)
433  for (ordinal_type i=0; i<d; i++)
434  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(basis_orders[i], true));
435  else
436  for (ordinal_type i=0; i<d; i++)
437  bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<ordinal_type,value_type>(basis_orders[i], 1.0, 2.0, true));
438  typedef Stokhos::LexographicLess< Stokhos::MultiIndex<int> > less_type;
439  typedef Stokhos::TotalOrderBasis<ordinal_type,value_type,less_type> basis_type;
440  Teuchos::RCP<const basis_type> basis =
441  Teuchos::rcp(new basis_type(bases, sparse_tol));
442 
443  // Build "standard" Cijk
444  Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
445  basis->computeTripleProductTensor();
446 
447  // Build "standard" algebraic expansion
448  Stokhos::AlgebraicOrthogPolyExpansion<ordinal_type, value_type> expn(
449  basis, Cijk);
450 
451  // Build quadrature expansion for sin/cos
452  Teuchos::RCP<const Stokhos::Quadrature<ordinal_type,value_type> > quad =
453  Teuchos::rcp(new Stokhos::TensorProductQuadrature<ordinal_type,value_type>(basis));
454  Stokhos::QuadOrthogPolyExpansion<ordinal_type, value_type> quad_expn(
455  basis, Cijk, quad);
456 
457  // Build flat LTB 3 tensor
458  Teuchos::RCP< Stokhos::FlatLTBSparse3Tensor<ordinal_type,value_type> >
459  flat_Cijk =
460  Stokhos::computeFlatTripleProductTensorLTB(*basis, symmetric);
461 
462  // Expansions
463  Stokhos::OrthogPolyApprox<ordinal_type,value_type>
464  x(basis), a(basis), b(basis), c1(basis), c2(basis);
465 
466  // Initialize a and b to reasonable values
467  x.term(0, 0) = 1.0;
468  for (ordinal_type i=0; i<d; i++)
469  x.term(i, 1) = 0.1;
470  quad_expn.sin(a,x);
471  quad_expn.cos(b,x);
472 
473  // Do multiplications
474  expn.times(c1,a,b);
475  Stokhos::flatLTB3TensorMultiply<10>(c2, a, b, *flat_Cijk);
476 
477  // Test c1 == c2
478  success = Stokhos::comparePCEs(c1, "c1", c2, "c2", rtol, atol, out);
479 
480  return success;
481  }
482 
483  TEUCHOS_UNIT_TEST( LexicographicTreeCoefficients, Isotropic ) {
484  Teuchos::Array<ordinal_type> basis_orders(setup.d, setup.p);
485  success = test_lexicographic_tree_coeffs(basis_orders, setup.p, true, out);
486  }
487 
488  TEUCHOS_UNIT_TEST( LexicographicTreeCoefficients, Anisotropic ) {
489  Teuchos::Array<ordinal_type> basis_orders(setup.d);
490  for (ordinal_type i=0; i<setup.d; ++i)
491  basis_orders[i] = i+1;
492  success = test_lexicographic_tree_coeffs(basis_orders, setup.d, false, out);
493  }
494 
495  TEUCHOS_UNIT_TEST( LTBSparse3Tensor, Isotropic_Symmetric ) {
496  Teuchos::Array<ordinal_type> basis_orders(setup.d, setup.p);
497  success = test_lexicographic_tree_sparse_3_tensor(
498  basis_orders, setup.p, true, setup.sparse_tol, setup.atol, setup.rtol,
499  out);
500  }
501 
502  TEUCHOS_UNIT_TEST( LTBSparse3Tensor, Anisotropic_Symmetric ) {
503  Teuchos::Array<ordinal_type> basis_orders(setup.d);
504  for (ordinal_type i=0; i<setup.d; ++i)
505  basis_orders[i] = i+1;
506  success = test_lexicographic_tree_sparse_3_tensor(
507  basis_orders, setup.p, true, setup.sparse_tol, setup.atol, setup.rtol,
508  out);
509  }
510 
511  TEUCHOS_UNIT_TEST( LTBSparse3Tensor, Isotropic_Asymmetric ) {
512  Teuchos::Array<ordinal_type> basis_orders(setup.d, setup.p);
513  success = test_lexicographic_tree_sparse_3_tensor(
514  basis_orders, setup.p, false, setup.sparse_tol, setup.atol, setup.rtol,
515  out);
516  }
517 
518  TEUCHOS_UNIT_TEST( LTBSparse3Tensor, Anisotropic_Asymmetric ) {
519  Teuchos::Array<ordinal_type> basis_orders(setup.d);
520  for (ordinal_type i=0; i<setup.d; ++i)
521  basis_orders[i] = i+1;
522  success = test_lexicographic_tree_sparse_3_tensor(
523  basis_orders, setup.p, false, setup.sparse_tol, setup.atol, setup.rtol,
524  out);
525  }
526 
527  TEUCHOS_UNIT_TEST( LTBSparse3TensorBlock, Isotropic_Symmetric ) {
528  Teuchos::Array<ordinal_type> basis_orders(setup.d, setup.p);
529  success = test_lexicographic_tree_sparse_3_tensor_block(
530  basis_orders, setup.p, true, setup.sparse_tol, setup.atol, setup.rtol,
531  out);
532  }
533 
534  TEUCHOS_UNIT_TEST( LTBSparse3TensorBlockMultiply, Isotropic_Symmetric ) {
535  Teuchos::Array<ordinal_type> basis_orders(setup.d, setup.p);
536  success = test_lexicographic_tree_sparse_3_tensor_multiply(
537  basis_orders, setup.p, true, setup.sparse_tol, setup.atol, setup.rtol,
538  out);
539  }
540 
541  TEUCHOS_UNIT_TEST( LTBSparse3TensorBlockMultiply, Isotropic_Asymmetric ) {
542  Teuchos::Array<ordinal_type> basis_orders(setup.d, setup.p);
543  success = test_lexicographic_tree_sparse_3_tensor_multiply(
544  basis_orders, setup.p, false, setup.sparse_tol, setup.atol, setup.rtol,
545  out);
546  }
547 
548  TEUCHOS_UNIT_TEST( LTBSparse3TensorBlockMultiply, Anisotropic_Symmetric ) {
549  Teuchos::Array<ordinal_type> basis_orders(setup.d);
550  for (ordinal_type i=0; i<setup.d; ++i)
551  basis_orders[i] = i+1;
552  success = test_lexicographic_tree_sparse_3_tensor_multiply(
553  basis_orders, setup.p, true, setup.sparse_tol, setup.atol, setup.rtol,
554  out);
555  }
556 
557  TEUCHOS_UNIT_TEST( LTBSparse3TensorBlockMultiply, Anisotropic_Asymmetric ) {
558  Teuchos::Array<ordinal_type> basis_orders(setup.d);
559  for (ordinal_type i=0; i<setup.d; ++i)
560  basis_orders[i] = i+1;
561  success = test_lexicographic_tree_sparse_3_tensor_multiply(
562  basis_orders, setup.p, false, setup.sparse_tol, setup.atol, setup.rtol,
563  out);
564  }
565 }
566 
567 int main( int argc, char* argv[] ) {
568  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
569  int res = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
570  return res;
571 }
TotalOrderBasisUnitTest::test_lexicographic_tree_sparse_3_tensor_block
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)
Definition: Stokhos_LexicographicTreeBasisUnitTest.cpp:286
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:286
Stokhos::QuadOrthogPolyExpansion::sin
void sin(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition: Stokhos_QuadOrthogPolyExpansionImp.hpp:635
Stokhos::computeTripleProductTensorLTB
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Definition: Stokhos_LTBSparse3Tensor.hpp:226
Stokhos::n_choose_k
ordinal_type n_choose_k(const ordinal_type &n, const ordinal_type &k)
Compute bionomial coefficient (n ; k) = n!/( k! (n-k)! )
Definition: Stokhos_ProductBasisUtils.hpp:60
Stokhos::TotalOrderBasis
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Definition: Stokhos_TotalOrderBasis.hpp:68
Stokhos::ProductBasisUtils::Cijk_1D_Iterator
Definition: Stokhos_ProductBasisUtils.hpp:1293
Stokhos::QuadOrthogPolyExpansion::cos
void cos(OrthogPolyApprox< ordinal_type, value_type, node_type > &c, const OrthogPolyApprox< ordinal_type, value_type, node_type > &a)
Definition: Stokhos_QuadOrthogPolyExpansionImp.hpp:644
Stokhos::MultiIndex::dimension
ordinal_type dimension() const
Dimension.
Definition: Stokhos_ProductBasisUtils.hpp:96
Stokhos::comparePCEs
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)
Definition: Stokhos_SacadoUQPCEUnitTest.cpp:58
TotalOrderBasisUnitTest::test_lexicographic_tree_sparse_3_tensor
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)
Definition: Stokhos_LexicographicTreeBasisUnitTest.cpp:178
Stokhos::RecurrenceBasis::computeTripleProductTensor
virtual Teuchos::RCP< Stokhos::Dense3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
Definition: Stokhos_RecurrenceBasisImp.hpp:147
Stokhos::AlgebraicOrthogPolyExpansion
Orthogonal polynomial expansions limited to algebraic operations.
Definition: Stokhos_AlgebraicOrthogPolyExpansion.hpp:62
Stokhos::computeTripleProductTensorLTBBlockLeaf
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTBBlockLeaf(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Definition: Stokhos_LTBSparse3Tensor.hpp:368
Stokhos_LTBSparse3Tensor.hpp
basis_type
Stokhos::LegendreBasis< int, double > basis_type
Definition: gram_schmidt_example.cpp:52
Sacado::UQ::min
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:1247
node_type
Kokkos::Serial node_type
Definition: Stokhos_SacadoTraitsUnitTest.cpp:59
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Teuchos_UnitTestRepository.hpp
Stokhos::LTBSparse3Tensor
Data structure storing a sparse 3-tensor C(i,j,k) in a a tree-based format for lexicographically orde...
Definition: Stokhos_LTBSparse3Tensor.hpp:57
TotalOrderBasisUnitTest::TEUCHOS_UNIT_TEST
TEUCHOS_UNIT_TEST(LexicographicTreeCoefficients, Isotropic)
Definition: Stokhos_LexicographicTreeBasisUnitTest.cpp:483
TotalOrderBasisUnitTest::test_lexicographic_tree_coeffs
bool test_lexicographic_tree_coeffs(const Teuchos::Array< ordinal_type > &basis_orders, const ordinal_type p, const bool isotropic, Teuchos::FancyOStream &out)
Definition: Stokhos_LexicographicTreeBasisUnitTest.cpp:74
TotalOrderBasisUnitTest::test_lexicographic_tree_sparse_3_tensor_multiply
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)
Definition: Stokhos_LexicographicTreeBasisUnitTest.cpp:418
Stokhos::build_lexicographic_basis_tree
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 >())
Definition: Stokhos_LTBSparse3Tensor.hpp:172
Teuchos::UnitTestRepository::runUnitTestsFromMain
static int runUnitTestsFromMain(int argc, char *argv[])
Stokhos::computeFlatTripleProductTensorLTB
Teuchos::RCP< FlatLTBSparse3Tensor< ordinal_type, value_type > > computeFlatTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Definition: Stokhos_LTBSparse3Tensor.hpp:545
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::basic_FancyOStream
Teuchos_UnitTestHarness.hpp
Stokhos::JacobiBasis
Jacobi polynomial basis.
Definition: Stokhos_JacobiBasis.hpp:94
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1181
Stokhos::OrthogPolyExpansionBase< ordinal_type, value_type, Stokhos::StandardStorage< ordinal_type, value_type > >::times
void times(OrthogPolyApprox< ordinal_type, value_type, Stokhos::StandardStorage< ordinal_type, value_type > > &c, const OrthogPolyApprox< ordinal_type, value_type, Stokhos::StandardStorage< ordinal_type, value_type > > &a, const OrthogPolyApprox< ordinal_type, value_type, Stokhos::StandardStorage< ordinal_type, value_type > > &b)
Stokhos_UnitTestHelpers.hpp
Kokkos::cijk
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Definition: KokkosExp_View_UQ_PCE_Contiguous.hpp:217
Stokhos::LegendreBasis
Legendre polynomial basis.
Definition: Stokhos_LegendreBasis.hpp:67
main
int main(int argc, char **argv)
Definition: cijk_ltb_partition.cpp:57
Teuchos::Array::back
reference back()
Teuchos_TestingHelpers.hpp
Teuchos::Array::push_back
void push_back(const value_type &x)
TEUCHOS_TEST_EQUALITY
#define TEUCHOS_TEST_EQUALITY(v1, v2, out, success)
Teuchos_GlobalMPISession.hpp
Teuchos::Array::size
size_type size() const
Stokhos::LexographicLess
A comparison functor implementing a strict weak ordering based lexographic ordering.
Definition: Stokhos_ProductBasisUtils.hpp:799
Stokhos::RecurrenceBasis::size
virtual ordinal_type size() const
Return total size of basis (given by order() + 1).
Definition: Stokhos_RecurrenceBasisImp.hpp:123
Stokhos::QuadOrthogPolyExpansion
Orthogonal polynomial expansions based on numerical quadrature.
Definition: Stokhos_QuadOrthogPolyExpansion.hpp:62
Stokhos::TensorProductQuadrature
Defines quadrature for a tensor product basis by tensor products of 1-D quadrature rules...
Definition: Stokhos_TensorProductQuadrature.hpp:58
TotalOrderBasisUnitTest::setup
UnitTestSetup< ordinal_type, value_type > setup
Definition: Stokhos_LexicographicTreeBasisUnitTest.cpp:71
Generated by   doxygen 1.8.5