Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Stokhos_ProductBasisUtilsUnitTest.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 "Teuchos_ArrayView.hpp"
49 
50 #include "Stokhos.hpp"
52 
53 #include <iterator>
54 
56 
57  // Common setup for unit tests
58  template <typename OrdinalType, typename ValueType>
59  struct UnitTestSetup {
60  ValueType rtol, atol;
61  OrdinalType sz,p,d;
62 
64  rtol = 1e-12;
65  atol = 1e-12;
66  d = 3;
67  p = 5;
68  }
69 
70  };
71 
72  typedef int ordinal_type;
73  typedef double value_type;
75 
76  // Utility function for computing factorials
77  template <typename ordinal_type>
79  ordinal_type res = 1;
80  for (ordinal_type i=1; i<=n; ++i)
81  res *= i;
82  return res;
83  }
84 
85  // Function for testing quadratures
86  template <typename scalar_type>
87  scalar_type quad_func1(const Teuchos::Array<scalar_type>& x) {
88  scalar_type val = 0.0;
89  for (int i=0; i<x.size(); i++)
90  val += x[i];
91  return std::exp(val);
92  }
93 
94  // Function for testing quadratures
95  template <typename scalar_type>
96  scalar_type quad_func2(const Teuchos::Array<scalar_type>& x) {
97  scalar_type val = 0.0;
98  for (int i=0; i<x.size(); i++)
99  val += x[i];
100  return std::sin(val);
101  }
102 
103  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, NChooseK ) {
104  ordinal_type n, k, v1, v2;
105 
106  success = true;
107 
108  n = 7; k = 3; // n-k > k
109  v1 = Stokhos::n_choose_k(n,k);
110  v2 = factorial(n)/(factorial(k)*factorial(n-k));
111  if (v1 != v2) {
112  out << "For n =" << n << ", k = " << k << ", n_choose_k = " << v1
113  << " != " << v2 << std::endl;
114  success = false;
115  }
116 
117  n = 7; k = 4; // n-k < k
118  v1 = Stokhos::n_choose_k(n,k);
119  v2 = factorial(n)/(factorial(k)*factorial(n-k));
120  if (v1 != v2) {
121  out << "For n =" << n << ", k = " << k << ", n_choose_k = " << v1
122  << " != " << v2 << std::endl;
123  success = false;
124  }
125 
126  }
127 
128  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, TotalOrderLess ) {
129  success = true;
130 
131  // Build sorted index set of dimension d and order p
132  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
133  typedef index_set_type::multiindex_type multiindex_type;
135  typedef std::set<multiindex_type, less_type> multiindex_set;
136  typedef multiindex_set::iterator iterator;
137  index_set_type indexSet(setup.d, 0, setup.p);
138  multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
139 
140  // Print sorted index set
141  std::ostream_iterator<multiindex_type> out_iterator(out, "\n");
142  out << std::endl << "Sorted total order index set (dimension = " << setup.d
143  << ", order = " << setup.p << "):" << std::endl;
144  std::copy(sortedIndexSet.begin(), sortedIndexSet.end(), out_iterator);
145 
146  // Ensure orders of each index are increasing
147  iterator prev = sortedIndexSet.begin();
148  iterator curr = prev; ++curr;
149  while (curr != sortedIndexSet.end()) {
150  ordinal_type order_prev = prev->order();
151  ordinal_type order_curr = curr->order();
152  ordinal_type i = 0;
153  while (i < setup.d && order_prev == order_curr) {
154  order_prev -= (*prev)[i];
155  order_curr -= (*curr)[i];
156  ++i;
157  }
158  if (order_prev >= order_curr) {
159  out << "previous index " << *prev << " and current index "
160  << *curr << " are out of order" << std::endl;
161  success = false;
162  }
163  prev = curr;
164  ++curr;
165  }
166  }
167 
168  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, LexographicLess ) {
169  success = true;
170 
171  // Build sorted index set of dimension d and order p
172  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
173  typedef index_set_type::multiindex_type multiindex_type;
175  typedef std::set<multiindex_type, less_type> multiindex_set;
176  typedef multiindex_set::iterator iterator;
177  index_set_type indexSet(setup.d, 0, setup.p);
178  multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
179 
180  // Print sorted index set
181  std::ostream_iterator<multiindex_type> out_iterator(out, "\n");
182  out << std::endl << "Sorted lexicographic index set (dimension = "
183  << setup.d << ", order = " << setup.p << "):" << std::endl;
184  std::copy(sortedIndexSet.begin(), sortedIndexSet.end(), out_iterator);
185 
186  // Ensure orders of each index are increasing
187  iterator prev = sortedIndexSet.begin();
188  iterator curr = prev; ++curr;
189  while (curr != sortedIndexSet.end()) {
190  ordinal_type i = 0;
191  while (i < setup.d && (*prev)[i] == (*curr)[i]) ++i;
192  if (i == setup.d || (*prev)[i] >= (*curr)[i]) {
193  out << "previous index " << *prev << " and current index "
194  << *curr << " are out of order" << std::endl;
195  success = false;
196  }
197  prev = curr;
198  ++curr;
199  }
200  }
201 
202  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, FloatingPointLess ) {
203  success = true;
204 
205  value_type tol=1e-12;
207 
208  TEUCHOS_TEST_EQUALITY(less(-0.774597,-0.774597), false, out, success);
209  TEUCHOS_TEST_EQUALITY(less(-0.774597+tol/2.0,-0.774597), false, out, success);
210  TEUCHOS_TEST_EQUALITY(less(-0.774597-tol/2.0,-0.774597), false, out, success);
211  TEUCHOS_TEST_EQUALITY(less(-0.774597,-0.774597+tol/2.0), false, out, success);
212  TEUCHOS_TEST_EQUALITY(less(-0.774597,-0.774597-tol/2.0), false, out, success);
213  TEUCHOS_TEST_EQUALITY(less(-0.774597,0.0), true, out, success);
214  TEUCHOS_TEST_EQUALITY(less(0.0,-0.774597), false, out, success);
215  }
216 
217  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, LexographicFloatingPointLess ) {
218  success = true;
219 
221  typedef Stokhos::FloatingPointLess<value_type> comp_type;
222  term_type a(2), b(2);
224  a[0] = -0.774597; a[1] = -0.774597;
225  b[0] = -0.774597; b[1] = 0.0;
226 
227  TEUCHOS_TEST_EQUALITY(less(a,b), true, out, success);
228  TEUCHOS_TEST_EQUALITY(less(b,a), false, out, success);
229  }
230 
231  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, TotalOrderIndexSet ) {
232  success = true;
233 
234  // Build index set of dimension d and order p
235  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
236  typedef index_set_type::multiindex_type multiindex_type;
237  typedef index_set_type::iterator iterator;
238  index_set_type indexSet(setup.d, 0, setup.p);
239 
240  // Print index set
241  out << std::endl << "Total order index set (dimension = " << setup.d
242  << ", order = " << setup.p << "):" << std::endl;
243  std::ostream_iterator<multiindex_type> out_iterator(out, "\n");
244  std::copy(indexSet.begin(), indexSet.end(), out_iterator);
245 
246  // Verify each index lies appropriatly in the set
247  for (iterator i=indexSet.begin(); i!=indexSet.end(); ++i) {
248  if (i->order() < 0 || i->order() > setup.p) {
249  out << "index " << *i << " does not lie in total order set! "
250  << std::endl;
251  success = false;
252  }
253  }
254 
255  // Put indices in sorted container -- this will ensure there are no
256  // duplicates, if we get the right size
258  typedef std::set<multiindex_type, less_type> multiindex_set;
259  multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
260 
261  out << "sorted index set size = " << sortedIndexSet.size() << std::endl;
262  out << "expected index set size = "
263  << Stokhos::n_choose_k(setup.p+setup.d,setup.d) << std::endl;
264  if (static_cast<ordinal_type>(sortedIndexSet.size()) !=
266  success = false;
267  }
268 
269  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils,
270  AnisotropicTotalOrderIndexSet ) {
271  success = true;
272 
273  // Build index set of dimension d and order p
275  typedef index_set_type::multiindex_type multiindex_type;
276  typedef index_set_type::iterator iterator;
277  multiindex_type upper(setup.d);
278  for (ordinal_type i=0; i<setup.d; ++i)
279  upper[i] = i+1;
280  index_set_type indexSet(setup.p, upper);
281 
282  // Print index set
283  out << std::endl << "Anisotropic total order index set (dimension = "
284  << setup.d << ", order = " << setup.p << ", and component orders = "
285  << upper << "):" << std::endl;
286  std::ostream_iterator<multiindex_type> out_iterator(out, "\n");
287  std::copy(indexSet.begin(), indexSet.end(), out_iterator);
288 
289  // Verify each index lies appropriatly in the set
290  for (iterator i=indexSet.begin(); i!=indexSet.end(); ++i) {
291  if (i->order() < 0 || i->order() > setup.p || !i->termWiseLEQ(upper)) {
292  out << "index " << *i << " does not lie in total order set! "
293  << std::endl;
294  success = false;
295  }
296  }
297 
298  // Need to figure out how to compute the size of such an index set
299  /*
300  // Put indices in sorted container -- this will ensure there are no
301  // duplicates, if we get the right size
302  typedef Stokhos::TotalOrderLess<multiindex_type> less_type;
303  typedef std::set<multiindex_type, less_type> multiindex_set;
304  multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
305 
306  out << "sorted index set size = " << sortedIndexSet.size() << std::endl;
307  out << "expected index set size = "
308  << Stokhos::n_choose_k(setup.p+setup.d,setup.d) << std::endl;
309  if (static_cast<ordinal_type>(sortedIndexSet.size()) !=
310  Stokhos::n_choose_k(setup.p+setup.d,setup.d))
311  success = false;
312  */
313  }
314 
315  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, TotalOrderBasis ) {
316  success = true;
317 
318  // Build index set of dimension d and order p
319  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
320  index_set_type indexSet(setup.d, 0, setup.p);
321 
322  // Build total-order basis from index set
324  typedef Stokhos::TotalOrderLess<coeff_type> less_type;
325  typedef std::map<coeff_type, ordinal_type, less_type> basis_set_type;
326  typedef Teuchos::Array<coeff_type> basis_map_type;
327  basis_set_type basis_set;
328  basis_map_type basis_map;
330  indexSet, basis_set, basis_map);
331 
332  // Build total-order basis directly
333  ordinal_type sz;
334  Teuchos::Array< Stokhos::MultiIndex<ordinal_type> > terms;
335  Teuchos::Array<ordinal_type> num_terms;
337  compute_terms(setup.p, setup.d, sz, terms, num_terms);
338 
339  // Check sizes
340  TEUCHOS_TEST_EQUALITY(static_cast<ordinal_type>(basis_set.size()),
341  static_cast<ordinal_type>(basis_map.size()),
342  out, success);
343  TEUCHOS_TEST_EQUALITY(static_cast<ordinal_type>(basis_set.size()),
344  static_cast<ordinal_type>(terms.size()),
345  out, success);
346  TEUCHOS_TEST_EQUALITY(sz, static_cast<ordinal_type>(terms.size()),
347  out, success);
348 
349  std::ostream_iterator<ordinal_type> out_iterator(out, " ");
350  for (ordinal_type i=0; i<sz; i++) {
351 
352  // Verify terms match
353  out << "term " << basis_map[i] << " == " << terms[i] << " : ";
354  bool is_equal = true;
355  for (ordinal_type j=0; j<setup.d; j++)
356  is_equal = is_equal && terms[i][j] == basis_map[i][j];
357  if (is_equal)
358  out << "passed" << std::endl;
359  else {
360  out << "failed" << std::endl;
361  success = false;
362  }
363 
364  // Verify global index mapping matches
365  TEUCHOS_TEST_EQUALITY(basis_set[basis_map[i]], i, out, success);
366  }
367 
368  }
369 
370  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, TensorProductBasis ) {
371  success = true;
372 
373  // Build index set of dimension d and order p
374  typedef Stokhos::TensorProductIndexSet<ordinal_type> index_set_type;
375  index_set_type indexSet(setup.d, 0, setup.p);
376 
377  // Build total-order basis from index set
379  typedef Stokhos::TotalOrderLess<coeff_type> less_type;
380  typedef std::map<coeff_type, ordinal_type, less_type> basis_set_type;
381  typedef Teuchos::Array<coeff_type> basis_map_type;
382  basis_set_type basis_set;
383  basis_map_type basis_map;
385  indexSet, basis_set, basis_map);
386 
387  // Compute expected size
388  ordinal_type sz = 1;
389  for (ordinal_type i=0; i<setup.d; ++i)
390  sz *= setup.p+1;
391 
392  // Check sizes
393  TEUCHOS_TEST_EQUALITY(static_cast<ordinal_type>(basis_set.size()),
394  static_cast<ordinal_type>(basis_map.size()),
395  out, success);
396  TEUCHOS_TEST_EQUALITY(sz, static_cast<ordinal_type>(basis_set.size()),
397  out, success);
398 
399  std::ostream_iterator<ordinal_type> out_iterator(out, " ");
400  for (ordinal_type i=0; i<sz; i++) {
401 
402  // Verify terms match
403  out << "term " << basis_map[i] << " <= " << setup.p << " : ";
404  bool is_less = true;
405  for (ordinal_type j=0; j<setup.d; j++)
406  is_less = is_less && basis_map[i][j] <= setup.p;
407  if (is_less)
408  out << "passed" << std::endl;
409  else {
410  out << "failed" << std::endl;
411  success = false;
412  }
413 
414  // Verify global index mapping matches
415  TEUCHOS_TEST_EQUALITY(basis_set[basis_map[i]], i, out, success);
416  }
417 
418  }
419 
420  template <typename ordinal_type>
423 
425  dim(dim_), order(order_) {}
426 
427  template <typename term_type>
428  bool operator() (const term_type& term) const {
429  ordinal_type sum = 0;
430  for (ordinal_type i=0; i<dim; ++i)
431  sum += term[i];
432  return sum <= order;
433  }
434 
435  };
436 
437  template <typename basis_set_type>
439  const basis_set_type& basis_set;
440 
441  general_predicate(const basis_set_type& basis_set_) :
442  basis_set(basis_set_) {}
443 
444  template <typename term_type>
445  bool operator() (const term_type& term) const {
446  return basis_set.find(term) != basis_set.end();
447  }
448 
449  };
450 
451 
452  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, TotalOrderSparse3Tensor ) {
453  success = true;
454  ordinal_type dim = setup.d;
455  ordinal_type order = setup.p;
456 
457  // Build index set of dimension d and order p
458  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
459  index_set_type indexSet(dim, 0, order);
460 
461  // Build total-order basis from index set
463  typedef Stokhos::TotalOrderLess<coeff_type> less_type;
464  typedef std::map<coeff_type, ordinal_type, less_type> basis_set_type;
465  typedef Teuchos::Array<coeff_type> basis_map_type;
466  basis_set_type basis_set;
467  basis_map_type basis_map;
469  indexSet, basis_set, basis_map);
470 
471  // 1-D bases
472  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
473  for (ordinal_type i=0; i<dim; i++)
474  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(order, true));
475 
476  // Build Cijk tensor
478  total_order_predicate<ordinal_type> pred(dim, order);
479  //general_predicate<basis_set_type> pred(basis_set);
480  Teuchos::RCP<Cijk_type> Cijk =
481  Stokhos::ProductBasisUtils::computeTripleProductTensor<ordinal_type,value_type>(bases, basis_set, basis_map, pred, pred);
482 
483  // Build Cijk tensor using original approach
484  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<ordinal_type,value_type> > basis = Teuchos::rcp(new Stokhos::CompletePolynomialBasis<ordinal_type,value_type>(bases));
485  Teuchos::RCP<Cijk_type> Cijk2 =
486  basis->computeTripleProductTensor();
487 
488  // Check sizes
489  TEUCHOS_TEST_EQUALITY(Cijk->num_k(), Cijk2->num_k(), out, success);
490  TEUCHOS_TEST_EQUALITY(Cijk->num_entries(), Cijk2->num_entries(), out, success);
491 
492  // Check tensors match
493  for (Cijk_type::k_iterator k_it=Cijk2->k_begin();
494  k_it!=Cijk2->k_end(); ++k_it) {
495  int k = Stokhos::index(k_it);
496  for (Cijk_type::kj_iterator j_it = Cijk2->j_begin(k_it);
497  j_it != Cijk2->j_end(k_it); ++j_it) {
498  int j = Stokhos::index(j_it);
499  for (Cijk_type::kji_iterator i_it = Cijk2->i_begin(j_it);
500  i_it != Cijk2->i_end(j_it); ++i_it) {
501  int i = Stokhos::index(i_it);
502  double c = Cijk->getValue(i,j,k);
503  double c2 = Stokhos::value(i_it);
504  double tol = setup.atol + c2*setup.rtol;
505  double err = std::abs(c-c2);
506  bool s = err < tol;
507  if (!s) {
508  out << std::endl
509  << "Check: rel_err( C(" << i << "," << j << "," << k << ") )"
510  << " = " << "rel_err( " << c << ", " << c2 << " ) = " << err
511  << " <= " << tol << " : ";
512  if (s) out << "Passed.";
513  else out << "Failed!";
514  out << std::endl;
515  }
516  success = success && s;
517  }
518  }
519  }
520  }
521 
522  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, TotalOrderSparse3TensorNew ) {
523  success = true;
524  // ordinal_type dim = setup.d;
525  // ordinal_type order = setup.p;
526  ordinal_type dim = 2;
527  ordinal_type order = 3;
528 
529  // Build index set of dimension d and order p
530  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
531  index_set_type indexSet(dim, 0, order);
532 
533  // Build total-order basis from index set
535  typedef Stokhos::TotalOrderLess<coeff_type> less_type;
536  typedef std::map<coeff_type, ordinal_type, less_type> basis_set_type;
537  typedef Teuchos::Array<coeff_type> basis_map_type;
538  basis_set_type basis_set;
539  basis_map_type basis_map;
541  indexSet, basis_set, basis_map);
542 
543  // 1-D bases
544  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
545  for (ordinal_type i=0; i<dim; i++)
546  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(order, true));
547 
548  // Build Cijk tensor
550  total_order_predicate<ordinal_type> pred(dim, order);
551  //general_predicate<basis_set_type> pred(basis_set);
552  Teuchos::RCP<Cijk_type> Cijk =
553  Stokhos::ProductBasisUtils::computeTripleProductTensorNew<ordinal_type,value_type>(bases, basis_set, basis_map, pred, pred, false);
554 
555  // Build Cijk tensor using original approach
556  Teuchos::RCP<const Stokhos::CompletePolynomialBasis<ordinal_type,value_type> > basis = Teuchos::rcp(new Stokhos::CompletePolynomialBasis<ordinal_type,value_type>(bases));
557  Teuchos::RCP<Cijk_type> Cijk2 =
558  basis->computeTripleProductTensor();
559 
560  // Check sizes
561  TEUCHOS_TEST_EQUALITY(Cijk->num_k(), Cijk2->num_k(), out, success);
562  TEUCHOS_TEST_EQUALITY(Cijk->num_entries(), Cijk2->num_entries(), out, success);
563 
564  // Check tensors match
565  for (Cijk_type::k_iterator k_it=Cijk2->k_begin();
566  k_it!=Cijk2->k_end(); ++k_it) {
567  int k = Stokhos::index(k_it);
568  for (Cijk_type::kj_iterator j_it = Cijk2->j_begin(k_it);
569  j_it != Cijk2->j_end(k_it); ++j_it) {
570  int j = Stokhos::index(j_it);
571  for (Cijk_type::kji_iterator i_it = Cijk2->i_begin(j_it);
572  i_it != Cijk2->i_end(j_it); ++i_it) {
573  int i = Stokhos::index(i_it);
574  double c = Cijk->getValue(i,j,k);
575  double c2 = Stokhos::value(i_it);
576  double tol = setup.atol + c2*setup.rtol;
577  double err = std::abs(c-c2);
578  bool s = err < tol;
579  if (!s) {
580  out << std::endl
581  << "Check: rel_err( C(" << i << "," << j << "," << k << ") )"
582  << " = " << "rel_err( " << c << ", " << c2 << " ) = " << err
583  << " <= " << tol << " : ";
584  if (s) out << "Passed.";
585  else out << "Failed!";
586  out << std::endl;
587  }
588  success = success && s;
589  }
590  }
591  }
592  }
593 
594  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, TotalOrderSparse3LTO ) {
595  success = true;
596  ordinal_type dim = setup.d;
597  ordinal_type order = setup.p;
598  // ordinal_type dim = 2;
599  // ordinal_type order = 3;
600 
601  // 1-D bases
602  Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> > > bases(dim);
603  for (ordinal_type i=0; i<dim; i++)
604  bases[i] = Teuchos::rcp(new Stokhos::LegendreBasis<ordinal_type,value_type>(order, true));
605 
606  // Product basis
607  typedef Stokhos::MultiIndex<ordinal_type> coeff_type;
608  typedef Stokhos::LexographicLess<coeff_type> less_type;
610  bases);
611 
612  // Build Cijk tensor
614  Teuchos::RCP<Cijk_type> Cijk =
616 
617  // Build Cijk tensor using original approach
618  Teuchos::RCP<Cijk_type> Cijk2 =
620 
621  // Check sizes
622  TEUCHOS_TEST_EQUALITY(Cijk->num_k(), Cijk2->num_k(), out, success);
623  TEUCHOS_TEST_EQUALITY(Cijk->num_entries(), Cijk2->num_entries(), out, success);
624 
625  // Check tensors match
626  for (Cijk_type::k_iterator k_it=Cijk2->k_begin();
627  k_it!=Cijk2->k_end(); ++k_it) {
628  int k = Stokhos::index(k_it);
629  for (Cijk_type::kj_iterator j_it = Cijk2->j_begin(k_it);
630  j_it != Cijk2->j_end(k_it); ++j_it) {
631  int j = Stokhos::index(j_it);
632  for (Cijk_type::kji_iterator i_it = Cijk2->i_begin(j_it);
633  i_it != Cijk2->i_end(j_it); ++i_it) {
634  int i = Stokhos::index(i_it);
635  double c = Cijk->getValue(i,j,k);
636  double c2 = Stokhos::value(i_it);
637  double tol = setup.atol + c2*setup.rtol;
638  double err = std::abs(c-c2);
639  bool s = err < tol;
640  if (!s) {
641  out << std::endl
642  << "Check: rel_err( C(" << i << "," << j << "," << k << ") )"
643  << " = " << "rel_err( " << c << ", " << c2 << " ) = " << err
644  << " <= " << tol << " : ";
645  if (s) out << "Passed.";
646  else out << "Failed!";
647  out << std::endl;
648  }
649  success = success && s;
650  }
651  }
652  }
653  }
654 
655  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, TotalOrderMapping ) {
656  success = true;
657 
658  // Build sorted index set of dimension d and order p
659  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
660  typedef index_set_type::multiindex_type multiindex_type;
662  typedef std::set<multiindex_type, less_type> multiindex_set;
663  typedef multiindex_set::iterator iterator;
664  index_set_type indexSet(setup.d, 0, setup.p);
665  multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
666 
667  // Loop over each index set element and test if mapping
668  // computes proper global index
669  iterator i = sortedIndexSet.begin();
670  ordinal_type idx = 0;
671  while (i != sortedIndexSet.end()) {
672  ordinal_type idx_mapping = totalOrderMapping(*i);
673  out << *i << ": index = " << idx << " mapped index = " << idx_mapping
674  << ": ";
675  if (idx == idx_mapping)
676  out << "passed";
677  else {
678  out << "failed";
679  success = false;
680  }
681  out << std::endl;
682  ++i;
683  ++idx;
684  }
685  }
686 
687  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, LexographicMapping ) {
688  success = true;
689 
690  // Build sorted index set of dimension d and order p
691  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
692  typedef index_set_type::multiindex_type multiindex_type;
694  typedef std::set<multiindex_type, less_type> multiindex_set;
695  typedef multiindex_set::iterator iterator;
696  index_set_type indexSet(setup.d, 0, setup.p);
697  multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
698 
699  // Loop over each index set element and test if mapping
700  // computes proper global index
701  iterator i = sortedIndexSet.begin();
702  ordinal_type idx = 0;
703  while (i != sortedIndexSet.end()) {
704  ordinal_type idx_mapping = lexicographicMapping(*i,setup.p);
705  out << *i << ": index = " << idx << " mapped index = " << idx_mapping
706  << ": ";
707  if (idx == idx_mapping)
708  out << "passed";
709  else {
710  out << "failed";
711  success = false;
712  }
713  out << std::endl;
714  ++i;
715  ++idx;
716  }
717  }
718 
719  TEUCHOS_UNIT_TEST( Stokhos_ProductBasisUtils, LexographicMapping2 ) {
720  success = true;
721 
722  // Build sorted index set of dimension d and order p
723  ordinal_type d = setup.d;
724  ordinal_type p = setup.p;
725  typedef Stokhos::TotalOrderIndexSet<ordinal_type> index_set_type;
726  typedef index_set_type::multiindex_type multiindex_type;
728  typedef std::set<multiindex_type, less_type> multiindex_set;
729  index_set_type indexSet(d, 0, p);
730  multiindex_set sortedIndexSet(indexSet.begin(), indexSet.end());
731 
732  // Loop over lexicographically sorted multi-indices, compute global index
733  // using combinatorial number system, and test if it is correct
734  multiindex_type index(d), num_terms(d), orders(d,p);
735  bool stop = false;
736  ordinal_type idx = 0;
737  index[d-1] = -1;
738  num_terms[d-1] = -1;
739  while (!stop) {
740  // Increment index to next lexicographic term
741  ordinal_type dim = d-1;
742  ++index[dim];
743  while (dim > 0 && index[dim] > orders[dim]) {
744  index[dim] = 0;
745  --dim;
746  ++index[dim];
747  }
748  for (ordinal_type i=dim+1; i<d; ++i)
749  orders[i] = orders[i-1] - index[i-1];
750 
751  if (index[dim] > orders[dim])
752  stop = true;
753  else {
754  // Update num_terms: num_terms[dim] = number of terms with
755  // order < index[dim] and dimension d-dim-1
756  num_terms[dim] += Stokhos::n_choose_k(orders[dim]-index[dim]+d-dim,
757  d-dim-1);
758  for (ordinal_type i=dim+1; i<d; ++i)
759  num_terms[i] = 0;
760 
761  // Compute global index
762  ordinal_type I = num_terms.order();
763  //ordinal_type I = lexicographicMapping(index, p);
764 
765  out << index << ": index = " << idx << " mapped index = " << I
766  << ": ";
767  if (idx == I)
768  out << "passed";
769  else {
770  out << "failed";
771  success = false;
772  }
773  out << std::endl;
774  }
775 
776  ++idx;
777  }
778  }
779 }
780 
781 int main( int argc, char* argv[] ) {
782  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
783  int res = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
784  Teuchos::TimeMonitor::summarize(std::cout);
785  return res;
786 }
A functor for comparing floating-point numbers to some tolerance.
ordinal_type factorial(const ordinal_type &n)
ordinal_type n_choose_k(const ordinal_type &n, const ordinal_type &k)
Compute bionomial coefficient (n ; k) = n!/( k! (n-k)! )
Teuchos::RCP< Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTO(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Container storing a term in a generalized tensor product.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
expr val()
TEUCHOS_UNIT_TEST(Stokhos_ProductBasisUtils, NChooseK)
static void buildProductBasis(const index_set_type &index_set, const growth_rule_type &growth_rule, basis_set_type &basis_set, basis_map_type &basis_map)
Generate a product basis from an index set.
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< RD, RP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< XD, XP... > >::value >::type sum(const Kokkos::View< RD, RP... > &r, const Kokkos::View< XD, XP... > &x)
int main(int argc, char *argv[])
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
static ordinal_type compute_terms(ordinal_type p, ordinal_type d, ordinal_type &sz, Teuchos::Array< MultiIndex< ordinal_type > > &terms, Teuchos::Array< ordinal_type > &num_terms)
Compute the 2-D array of basis terms which maps a basis index into the orders for each basis dimensio...
UnitTestSetup< ordinal_type, value_type > setup
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
scalar_type quad_func2(const Teuchos::Array< scalar_type > &x)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr expr expr expr j
An anisotropic total order index set.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
KOKKOS_INLINE_FUNCTION PCE< Storage > exp(const PCE< Storage > &a)
Legendre polynomial basis.
Stokhos::Sparse3Tensor< int, double > Cijk_type
scalar_type quad_func1(const Teuchos::Array< scalar_type > &x)
An isotropic total order index set.
A comparison functor implementing a strict weak ordering based lexographic ordering.
KOKKOS_INLINE_FUNCTION PCE< Storage > sin(const PCE< Storage > &a)
ordinal_type lexicographicMapping(const Stokhos::MultiIndex< ordinal_type > &index, ordinal_type max_order)
ordinal_type totalOrderMapping(const Stokhos::MultiIndex< ordinal_type > &index)