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