Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_ProductBasisUtils.hpp
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#ifndef STOKHOS_PRODUCT_BASIS_UTILS_HPP
43#define STOKHOS_PRODUCT_BASIS_UTILS_HPP
44
45#include "Teuchos_Array.hpp"
46#include "Teuchos_RCP.hpp"
47#include "Teuchos_SerialDenseVector.hpp"
48#include "Teuchos_TimeMonitor.hpp"
49
50#include "Stokhos_SDMUtils.hpp"
53
54namespace Stokhos {
55
59 template <typename ordinal_type>
60 ordinal_type n_choose_k(const ordinal_type& n, const ordinal_type& k) {
61 // Use formula
62 // n!/(k!(n-k)!) = n(n-1)...(k+1) / ( (n-k)(n-k-1)...1 ) ( n-k terms )
63 // = n(n-1)...(n-k+1) / ( k(k-1)...1 ) ( k terms )
64 // which ever has fewer terms
65 if (k > n)
66 return 0;
67 ordinal_type num = 1;
68 ordinal_type den = 1;
69 ordinal_type l = std::min(n-k,k);
70 for (ordinal_type i=0; i<l; i++) {
71 num *= n-i;
72 den *= i+1;
73 }
74 return num / den;
75 }
76
78 template <typename ordinal_t>
79 class MultiIndex {
80 public:
81
82 typedef ordinal_t ordinal_type;
83 typedef ordinal_t element_type;
84
87
90 index(dim,v) {}
91
94
96 ordinal_type dimension() const { return index.size(); }
97
99 ordinal_type size() const { return index.size(); }
100
102 const ordinal_type& operator[] (ordinal_type i) const { return index[i]; }
103
106
108 const Teuchos::Array<element_type>& getTerm() const { return index; }
109
111 Teuchos::Array<element_type>& getTerm() { return index; }
112
115 for (ordinal_type i=0; i<dimension(); i++)
116 index[i] = v;
117 }
118
121 index.resize(d,v);
122 }
123
126 ordinal_type my_order = 0;
127 for (ordinal_type i=0; i<dimension(); ++i) my_order += index[i];
128 return my_order;
129 }
130
132 bool operator==(const MultiIndex& idx) const {
133 if (dimension() != idx.dimension())
134 return false;
135 for (ordinal_type i=0; i<dimension(); i++) {
136 if (index[i] != idx.index[i])
137 return false;
138 }
139 return true;
140 }
141
143 bool operator!=(const MultiIndex& idx) const { return !(*this == idx); }
144
146 bool termWiseLEQ(const MultiIndex& idx) const {
147 for (ordinal_type i=0; i<dimension(); i++) {
148 if (index[i] > idx.index[i])
149 return false;
150 }
151 return true;
152 }
153
155 std::ostream& print(std::ostream& os) const {
156 os << "[ ";
157 for (ordinal_type i=0; i<dimension(); i++)
158 os << index[i] << " ";
159 os << "]";
160 return os;
161 }
162
165 for (ordinal_type i=0; i<dimension(); i++)
166 index[i] = index[i] <= idx[i] ? index[i] : idx[i];
167 return *this;
168 }
169
172 for (ordinal_type i=0; i<dimension(); i++)
173 index[i] = index[i] <= idx ? index[i] : idx;
174 return *this;
175 }
176
179 for (ordinal_type i=0; i<dimension(); i++)
180 index[i] = index[i] >= idx[i] ? index[i] : idx[i];
181 return *this;
182 }
183
186 for (ordinal_type i=0; i<dimension(); i++)
187 index[i] = index[i] >= idx ? index[i] : idx;
188 return *this;
189 }
190
191 protected:
192
194 Teuchos::Array<ordinal_type> index;
195
196 };
197
198 template <typename ordinal_type>
199 std::ostream& operator << (std::ostream& os,
200 const MultiIndex<ordinal_type>& m) {
201 return m.print(os);
202 }
203
205
214 template <typename ordinal_t>
216 public:
217
218 // Forward declaration of our iterator
219 class Iterator;
220
221 typedef ordinal_t ordinal_type;
225
227
232 ordinal_type lower_,
233 ordinal_type upper_) :
234 dim(dim_), lower(lower_), upper(upper_) {}
235
237
242 ordinal_type upper_) :
243 dim(dim_), lower(0), upper(upper_) {}
244
246 ordinal_type dimension() const { return dim; }
247
250 return multiindex_type(dim, upper);
251 }
252
255 multiindex_type index(dim);
256 index[0] = lower;
257 return Iterator(upper, index);
258 }
259
262 multiindex_type index(dim);
263 index[dim-1] = upper+1;
264 return Iterator(upper, index);
265 }
266
267 protected:
268
271
274
277
278 public:
279
281 class Iterator : public std::iterator<std::input_iterator_tag,
282 multiindex_type> {
283 public:
284
285 typedef std::iterator<std::input_iterator_tag,multiindex_type> base_type;
286 typedef typename base_type::iterator_category iterator_category;
287 typedef typename base_type::value_type value_type;
288 typedef typename base_type::difference_type difference_type;
289 typedef typename base_type::reference reference;
290 typedef typename base_type::pointer pointer;
291
294
296
300 Iterator(ordinal_type max_order_, const multiindex_type& index_) :
301 max_order(max_order_), index(index_), dim(index.dimension()),
302 orders(dim) {
303 orders[dim-1] = max_order;
304 for (ordinal_type i=dim-2; i>=0; --i)
305 orders[i] = orders[i+1] - index[i+1];
306 }
307
309 bool operator==(const Iterator& it) const { return index == it.index; }
310
312 bool operator!=(const Iterator& it) const { return index != it.index; }
313
315 const_reference operator*() const { return index; }
316
318 const_pointer operator->() const { return &index; }
319
321
331 ++index[0];
332 ordinal_type i=0;
333 while (i<dim-1 && index[i] > orders[i]) {
334 index[i] = 0;
335 ++i;
336 ++index[i];
337 }
338 for (ordinal_type ii=dim-2; ii>=0; --ii)
339 orders[ii] = orders[ii+1] - index[ii+1];
340
341 return *this;
342 }
343
346 Iterator tmp(*this);
347 ++(*this);
348 return tmp;
349 }
350
351 protected:
352
355
358
361
363 Teuchos::Array<ordinal_type> orders;
364 };
365 };
366
368
377 template <typename ordinal_t>
379 public:
380
381 // Forward declaration of our iterator
382 class Iterator;
383
384 typedef ordinal_t ordinal_type;
388
390
395 const multiindex_type& lower_,
396 const multiindex_type& upper_) :
397 dim(lower_.dimension()),
398 upper_order(upper_order_),
399 lower(lower_),
400 upper(upper_) {}
401
403
408 const multiindex_type& upper_) :
409 dim(upper_.dimension()),
410 upper_order(upper_order_),
411 lower(dim,0),
412 upper(upper_) {}
413
415 ordinal_type dimension() const { return dim; }
416
418 multiindex_type max_orders() const { return upper; }
419
423 }
424
427 multiindex_type index(dim);
428 index[dim-1] = std::min(upper_order, upper[dim-1]) + 1;
429 return Iterator(upper_order, upper, index);
430 }
431
432 protected:
433
436
439
442
445
448
449 public:
450
452 class Iterator : public std::iterator<std::input_iterator_tag,
453 multiindex_type> {
454 public:
455
456 typedef std::iterator<std::input_iterator_tag,multiindex_type> base_type;
457 typedef typename base_type::iterator_category iterator_category;
458 typedef typename base_type::value_type value_type;
459 typedef typename base_type::difference_type difference_type;
460 typedef typename base_type::reference reference;
461 typedef typename base_type::pointer pointer;
462
465
467
472 const multiindex_type& component_max_order_,
473 const multiindex_type& index_) :
474 max_order(max_order_),
475 component_max_order(component_max_order_),
476 index(index_),
478 orders(dim)
479 {
480 orders[dim-1] = max_order;
481 for (ordinal_type i=dim-2; i>=0; --i)
482 orders[i] = orders[i+1] - index[i+1];
483 }
484
486 bool operator==(const Iterator& it) const { return index == it.index; }
487
489 bool operator!=(const Iterator& it) const { return index != it.index; }
490
492 const_reference operator*() const { return index; }
493
495 const_pointer operator->() const { return &index; }
496
498
508 ++index[0];
509 ordinal_type i=0;
510 while (i<dim-1 && (index[i] > orders[i] || index[i] > component_max_order[i])) {
511 index[i] = 0;
512 ++i;
513 ++index[i];
514 }
515 for (ordinal_type ii=dim-2; ii>=0; --ii)
516 orders[ii] = orders[ii+1] - index[ii+1];
517
518 return *this;
519 }
520
523 Iterator tmp(*this);
524 ++(*this);
525 return tmp;
526 }
527
528 protected:
529
532
535
538
541
543 Teuchos::Array<ordinal_type> orders;
544 };
545 };
546
548
557 template <typename ordinal_t>
559 public:
560
561 // Forward declaration of our iterator
562 class Iterator;
563
564 typedef ordinal_t ordinal_type;
568
570
575 ordinal_type lower_,
576 ordinal_type upper_) :
577 dim(dim_), lower(dim_,lower_), upper(dim_,upper_) {}
578
580
585 ordinal_type upper_) :
586 dim(dim_), lower(dim_,ordinal_type(0)), upper(dim_,upper_) {}
587
589
594 const multiindex_type& upper_) :
595 dim(lower_.dimension()), lower(lower_), upper(upper_) {}
596
598
603 dim(upper_.dimension()), lower(dim,ordinal_type(0)), upper(upper_) {}
604
606 ordinal_type dimension() const { return dim; }
607
610 return upper;
611 }
612
615 return Iterator(upper, lower);
616 }
617
620 multiindex_type index(dim);
621 index[dim-1] = upper[dim-1]+1;
622 return Iterator(upper, index);
623 }
624
625 protected:
626
629
632
635
636 public:
637
639 class Iterator : public std::iterator<std::input_iterator_tag,
640 multiindex_type> {
641 public:
642
643 typedef std::iterator<std::input_iterator_tag,multiindex_type> base_type;
644 typedef typename base_type::iterator_category iterator_category;
645 typedef typename base_type::value_type value_type;
646 typedef typename base_type::difference_type difference_type;
647 typedef typename base_type::reference reference;
648 typedef typename base_type::pointer pointer;
649
652
654
658 Iterator(const multiindex_type& upper_, const multiindex_type& index_) :
659 upper(upper_), index(index_), dim(index.dimension()) {}
660
662 bool operator==(const Iterator& it) const { return index == it.index; }
663
665 bool operator!=(const Iterator& it) const { return index != it.index; }
666
668 const_reference operator*() const { return index; }
669
671 const_pointer operator->() const { return &index; }
672
674
684 ++index[0];
685 ordinal_type i=0;
686 while (i<dim-1 && index[i] > upper[i]) {
687 index[i] = 0;
688 ++i;
689 ++index[i];
690 }
691 return *this;
692 }
693
696 Iterator tmp(*this);
697 ++(*this);
698 return tmp;
699 }
700
701 protected:
702
705
708
711
712 };
713 };
714
716 template <typename ordinal_t, typename element_t>
718 public:
719
720 typedef ordinal_t ordinal_type;
721 typedef element_t element_type;
722
725
728 const element_type& val = element_type(0)) :
729 term(dim,val) {}
730
733
735 ordinal_type dimension() const { return term.size(); }
736
738 ordinal_type size() const { return term.size(); }
739
741 const element_type& operator[] (ordinal_type i) const { return term[i]; }
742
745
747 const Teuchos::Array<element_type>& getTerm() const { return term; }
748
750 Teuchos::Array<element_type>& getTerm() { return term; }
751
753 operator Teuchos::ArrayView<element_type>() { return term; }
754
756 operator Teuchos::ArrayView<const element_type>() const { return term; }
757
760 element_type my_order = 0;
761 for (ordinal_type i=0; i<dimension(); ++i) my_order += term[i];
762 return my_order;
763 }
764
766 std::ostream& print(std::ostream& os) const {
767 os << "[ ";
768 for (ordinal_type i=0; i<dimension(); i++)
769 os << term[i] << " ";
770 os << "]";
771 return os;
772 }
773
774 protected:
775
777 Teuchos::Array<element_type> term;
778
779 };
780
781 template <typename ordinal_type, typename element_type>
782 std::ostream& operator << (
783 std::ostream& os,
785 return m.print(os);
786 }
787
792 /*
793 * Objects of type \c term_type must implement \c dimension() and
794 * \c operator[] methods, as well as contain ordinal_type and element_type
795 * nested types.
796 */
797 template <typename term_type,
798 typename compare_type = std::less<typename term_type::element_type> >
800 public:
801
802 typedef term_type product_element_type;
803 typedef typename term_type::ordinal_type ordinal_type;
804 typedef typename term_type::element_type element_type;
805
807 LexographicLess(const compare_type& cmp_ = compare_type()) : cmp(cmp_) {}
808
810 bool operator()(const term_type& a, const term_type& b) const {
811 ordinal_type i=0;
812 while(i < a.dimension() && i < b.dimension() && equal(a[i],b[i])) i++;
813
814 // if a is shorter than b and the first a.dimension() elements agree
815 // then a is always less than b
816 if (i == a.dimension()) return i != b.dimension();
817
818 // if a is longer than b and the first b.dimension() elements agree
819 // then b is always less than a
820 if (i == b.dimension()) return false;
821
822 // a and b different at element i, a is less than b if a[i] < b[i]
823 return cmp(a[i],b[i]);
824 }
825
826 protected:
827
829 compare_type cmp;
830
832 bool equal(const element_type& a, const element_type& b) const {
833 return !(cmp(a,b) || cmp(b,a));
834 }
835
836 };
837
842 /*
843 * Objects of type \c term_type must implement \c dimension(), \c order, and
844 * \c operator[] methods, as well as contain ordinal_type and element_type
845 * nested types.
846 */
847 template <typename term_type,
848 typename compare_type = std::less<typename term_type::element_type> >
850 public:
851
852 typedef term_type product_element_type;
853 typedef typename term_type::ordinal_type ordinal_type;
854 typedef typename term_type::element_type element_type;
855
857 TotalOrderLess(const compare_type& cmp_ = compare_type()) : cmp(cmp_) {}
858
859 bool operator()(const term_type& a, const term_type& b) const {
860 element_type a_order = a.order();
861 element_type b_order = b.order();
862 ordinal_type i=0;
863 while (i < a.dimension() && i < b.dimension() && equal(a_order,b_order)) {
864 a_order -= a[i];
865 b_order -= b[i];
866 ++i;
867 }
868 return cmp(a_order,b_order);
869 }
870
871 protected:
872
874 compare_type cmp;
875
877 bool equal(const element_type& a, const element_type& b) const {
878 return !(cmp(a,b) || cmp(b,a));
879 }
880
881 };
882
887 /*
888 * A Morton Z-ordering is based on thinking of terms of dimension d as
889 * coordinates in a d-dimensional space, and forming a linear ordering of
890 * the points in that space by interleaving the bits from each coordinate.
891 * As implemented, this ordering only works with integral types and
892 * directly uses <. The code was taken from here:
893 *
894 * http://en.wikipedia.org/wiki/Z-order_curve
895 *
896 * With the idea behind it published here:
897 *
898 * Chan, T. (2002), "Closest-point problems simplified on the RAM",
899 * ACM-SIAM Symposium on Discrete Algorithms.
900 */
901 template <typename term_type>
903 public:
904
905 typedef term_type product_element_type;
906 typedef typename term_type::ordinal_type ordinal_type;
907 typedef typename term_type::element_type element_type;
908
911
912 bool operator()(const term_type& a, const term_type& b) const {
913 ordinal_type d = a.dimension();
915 //ordinal_type k = ordinal_type(0); // causes shadowing warning
917 for (ordinal_type k=0; k<d; ++k) {
918 element_type y = a[k] ^ b[k];
919 if ( (x < y) && (x < (x^y)) ) {
920 j = k;
921 x = y;
922 }
923 }
924 return a[j] < b[j];
925 }
926
927 };
928
930
934 template <typename value_type>
936 public:
937
939 FloatingPointLess(const value_type& tol_ = 1.0e-12) : tol(tol_) {}
940
943
945 bool operator() (const value_type& a, const value_type& b) const {
946 return a < b - tol;
947 }
948
949 protected:
950
952 value_type tol;
953
954 };
955
957 template <typename ordinal_type>
961
962 TensorProductPredicate(const coeff_type& orders_) : orders(orders_) {}
963
964 bool operator() (const coeff_type& term) const {
965 return term.termWiseLEQ(orders);
966 }
967
968 };
969
971 template <typename ordinal_type>
974 ordinal_type max_order;
976
977 TotalOrderPredicate(ordinal_type max_order_, const coeff_type& orders_)
978 : max_order(max_order_), orders(orders_) {}
979
980 bool operator() (const coeff_type& term) const {
981 return term.termWiseLEQ(orders) && term.order() <= max_order;
982 }
983
984 };
985
986 // Compute global index from a total-order-sorted multi-index
987 /*
988 * The strategy is based on the fact we can compute the number of terms
989 * of a total order expansion for any given order and dimension. Thus
990 * starting from the last dimension of the multi-index:
991 * -- compute the order p of the multi-index from the current dimension
992 * dim-d to the last
993 * -- compute the number of terms in an order p-1 expansion of dimension
994 * d
995 * -- add this number of terms to the global index
996 * The logic here is independent of the values of the multi-index and
997 * requires exactly 8*d + (3/2)*d*(d+1) integer operations and comparisons
998 */
999 template <typename ordinal_type>
1000 ordinal_type
1002 ordinal_type dim = index.dimension();
1003 ordinal_type idx = ordinal_type(0);
1004 ordinal_type p = ordinal_type(0);
1005 ordinal_type den = ordinal_type(1);
1006 for (ordinal_type d=1; d<=dim; ++d) {
1007 p += index[dim-d];
1008
1009 // offset = n_choose_k(p-1+d,d) = (p-1+d)! / ( (p-1)!*d! ) =
1010 // ( p*(p+1)*...*(p+d-1) )/( 1*2*...*d )
1011 den *= d;
1012 ordinal_type num = 1;
1013 for (ordinal_type i=p; i<p+d; ++i)
1014 num *= i;
1015
1016 idx += num / den;
1017 }
1018 return idx;
1019 }
1020
1021 // Compute global index from a lexicographic-sorted multi-index
1022 /*
1023 * For a given dimension d, let p be the sum of the orders for dimensions
1024 * <= d. Then the number of terms that agree with the first d entries
1025 * is (max_order-p+dim-d)!/(max_order-p)!*(dim-d)! where max_order is
1026 * the maximum order of the multi-index and dim is the dimension. Thus
1027 * we loop over each dimension d and compute the number of terms that come
1028 * before that term that agree with the first d dimensions.
1029 */
1030 template <typename ordinal_type>
1031 ordinal_type
1033 ordinal_type max_order) {
1034 ordinal_type dim = index.dimension();
1035 ordinal_type idx = ordinal_type(0);
1036 for (ordinal_type d=1; d<=dim; ++d) {
1037 ordinal_type p = index[d-1];
1038
1039 ordinal_type offset = ordinal_type(0);
1040 for (ordinal_type pp=0; pp<p; ++pp)
1041 offset += Stokhos::n_choose_k(max_order-pp+dim-d,dim-d);
1042
1043 idx += offset;
1044 max_order -= p;
1045 }
1046 return idx;
1047 }
1048
1058 public:
1059
1063 template <typename index_set_type,
1064 typename growth_rule_type,
1065 typename basis_set_type,
1066 typename basis_map_type>
1067 static void
1068 buildProductBasis(const index_set_type& index_set,
1069 const growth_rule_type& growth_rule,
1070 basis_set_type& basis_set,
1071 basis_map_type& basis_map) {
1072
1073 typedef typename index_set_type::ordinal_type ordinal_type;
1074 typedef typename index_set_type::iterator index_iterator_type;
1075 typedef typename basis_set_type::iterator basis_iterator_type;
1076 typedef typename basis_set_type::key_type coeff_type;
1077
1078 ordinal_type dim = index_set.dimension();
1079
1080 // Iterator over elements of index set
1081 index_iterator_type index_iterator = index_set.begin();
1082 index_iterator_type index_iterator_end = index_set.end();
1083 for (; index_iterator != index_iterator_end; ++index_iterator) {
1084
1085 // Generate product coefficient
1086 coeff_type coeff(dim);
1087 for (ordinal_type i=0; i<dim; i++)
1088 coeff[i] = growth_rule[i]((*index_iterator)[i]);
1089
1090 // Insert coefficient into set
1091 basis_set[coeff] = ordinal_type(0);
1092 }
1093
1094 // Generate linear ordering of basis_set elements
1095 basis_map.resize(basis_set.size());
1096 ordinal_type idx = 0;
1097 basis_iterator_type basis_iterator = basis_set.begin();
1098 basis_iterator_type basis_iterator_end = basis_set.end();
1099 for (; basis_iterator != basis_iterator_end; ++basis_iterator) {
1100 basis_iterator->second = idx;
1101 basis_map[idx] = basis_iterator->first;
1102 ++idx;
1103 }
1104 }
1105
1110 template <typename index_set_type,
1111 typename basis_set_type,
1112 typename basis_map_type>
1113 static void
1114 buildProductBasis(const index_set_type& index_set,
1115 basis_set_type& basis_set,
1116 basis_map_type& basis_map) {
1117 typedef typename index_set_type::ordinal_type ordinal_type;
1118 ordinal_type dim = index_set.dimension();
1119 Teuchos::Array< IdentityGrowthRule<ordinal_type> > growth_rule(dim);
1120 buildProductBasis(index_set, growth_rule, basis_set, basis_map);
1121 }
1122
1123 template <typename ordinal_type,
1124 typename value_type,
1125 typename basis_set_type,
1126 typename basis_map_type,
1127 typename coeff_predicate_type,
1128 typename k_coeff_predicate_type>
1129 static Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
1131 const Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases,
1132 const basis_set_type& basis_set,
1133 const basis_map_type& basis_map,
1134 const coeff_predicate_type& coeff_pred,
1135 const k_coeff_predicate_type& k_coeff_pred,
1136 const value_type sparse_tol = 1.0e-12)
1137 {
1138#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
1139 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Time");
1140#endif
1141 typedef typename basis_map_type::value_type coeff_type;
1142 ordinal_type d = bases.size();
1143
1144 // The algorithm for computing Cijk = < \Psi_i \Psi_j \Psi_k > here
1145 // works by factoring
1146 // < \Psi_i \Psi_j \Psi_k > =
1147 // < \psi^1_{i_1}\psi^1_{j_1}\psi^1_{k_1} >_1 x ... x
1148 // < \psi^d_{i_d}\psi^d_{j_d}\psi^d_{k_d} >_d
1149 // We compute the sparse triple product < \psi^l_i\psi^l_j\psi^l_k >_l
1150 // for each dimension l, and then compute all non-zero products of these
1151 // terms. The complexity arises from iterating through all possible
1152 // combinations, throwing out terms that aren't in the basis and are
1153 // beyond the k-order limit provided
1154 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
1156
1157 // Create 1-D triple products
1158 Teuchos::Array< Teuchos::RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
1159 for (ordinal_type i=0; i<d; i++) {
1160 Cijk_1d[i] =
1161 bases[i]->computeSparseTripleProductTensor(bases[i]->order()+1);
1162 }
1163
1164 // Create i, j, k iterators for each dimension
1165 // Note: we have to supply an initializer in the arrays of iterators
1166 // to avoid checked-stl errors about singular iterators
1168 typedef typename Cijk_type::k_iterator k_iterator;
1169 typedef typename Cijk_type::kj_iterator kj_iterator;
1170 typedef typename Cijk_type::kji_iterator kji_iterator;
1171 Teuchos::Array<k_iterator> k_iterators(d, Cijk_1d[0]->k_begin());
1172 Teuchos::Array<kj_iterator > j_iterators(d, Cijk_1d[0]->j_begin(k_iterators[0]));
1173 Teuchos::Array<kji_iterator > i_iterators(d, Cijk_1d[0]->i_begin(j_iterators[0]));
1174 coeff_type terms_i(d), terms_j(d), terms_k(d);
1175 for (ordinal_type dim=0; dim<d; dim++) {
1176 k_iterators[dim] = Cijk_1d[dim]->k_begin();
1177 j_iterators[dim] = Cijk_1d[dim]->j_begin(k_iterators[dim]);
1178 i_iterators[dim] = Cijk_1d[dim]->i_begin(j_iterators[dim]);
1179 terms_i[dim] = index(i_iterators[dim]);
1180 terms_j[dim] = index(j_iterators[dim]);
1181 terms_k[dim] = index(k_iterators[dim]);
1182 }
1183
1184 ordinal_type I = 0;
1185 ordinal_type J = 0;
1186 ordinal_type K = 0;
1187 bool valid_i = coeff_pred(terms_i);
1188 bool valid_j = coeff_pred(terms_j);
1189 bool valid_k = k_coeff_pred(terms_k);
1190 bool inc_i = true;
1191 bool inc_j = true;
1192 bool inc_k = true;
1193 bool stop = false;
1194 ordinal_type cnt = 0;
1195 while (!stop) {
1196
1197 // Add term if it is in the basis
1198 if (valid_i && valid_j && valid_k) {
1199 if (inc_k) {
1200 typename basis_set_type::const_iterator k =
1201 basis_set.find(terms_k);
1202 K = k->second;
1203 }
1204 if (inc_i) {
1205 typename basis_set_type::const_iterator i =
1206 basis_set.find(terms_i);
1207 I = i->second;
1208 }
1209 if (inc_j) {
1210 typename basis_set_type::const_iterator j =
1211 basis_set.find(terms_j);
1212 J = j->second;
1213 }
1214 value_type c = value_type(1.0);
1215 value_type nrm = value_type(1.0);
1216 for (ordinal_type dim=0; dim<d; dim++) {
1217 c *= value(i_iterators[dim]);
1218 nrm *= bases[dim]->norm_squared(terms_i[dim]);
1219 }
1220 if (std::abs(c/nrm) > sparse_tol)
1221 Cijk->add_term(I,J,K,c);
1222 }
1223
1224 // Increment iterators to the next valid term
1225 ordinal_type cdim = 0;
1226 bool inc = true;
1227 inc_i = false;
1228 inc_j = false;
1229 inc_k = false;
1230 while (inc && cdim < d) {
1231 ordinal_type cur_dim = cdim;
1232 ++i_iterators[cdim];
1233 inc_i = true;
1234 if (i_iterators[cdim] != Cijk_1d[cdim]->i_end(j_iterators[cdim])) {
1235 terms_i[cdim] = index(i_iterators[cdim]);
1236 valid_i = coeff_pred(terms_i);
1237 }
1238 if (i_iterators[cdim] == Cijk_1d[cdim]->i_end(j_iterators[cdim]) ||
1239 !valid_i) {
1240 ++j_iterators[cdim];
1241 inc_j = true;
1242 if (j_iterators[cdim] != Cijk_1d[cdim]->j_end(k_iterators[cdim])) {
1243 terms_j[cdim] = index(j_iterators[cdim]);
1244 valid_j = coeff_pred(terms_j);
1245 }
1246 if (j_iterators[cdim] == Cijk_1d[cdim]->j_end(k_iterators[cdim]) ||
1247 !valid_j) {
1248 ++k_iterators[cdim];
1249 inc_k = true;
1250 if (k_iterators[cdim] != Cijk_1d[cdim]->k_end()) {
1251 terms_k[cdim] = index(k_iterators[cdim]);
1252 valid_k = k_coeff_pred(terms_k);
1253 }
1254 if (k_iterators[cdim] == Cijk_1d[cdim]->k_end() || !valid_k) {
1255 k_iterators[cdim] = Cijk_1d[cdim]->k_begin();
1256 ++cdim;
1257 terms_k[cur_dim] = index(k_iterators[cur_dim]);
1258 valid_k = k_coeff_pred(terms_k);
1259 }
1260 else
1261 inc = false;
1262 j_iterators[cur_dim] =
1263 Cijk_1d[cur_dim]->j_begin(k_iterators[cur_dim]);
1264 terms_j[cur_dim] = index(j_iterators[cur_dim]);
1265 valid_j = coeff_pred(terms_j);
1266 }
1267 else
1268 inc = false;
1269 i_iterators[cur_dim] =
1270 Cijk_1d[cur_dim]->i_begin(j_iterators[cur_dim]);
1271 terms_i[cur_dim] = index(i_iterators[cur_dim]);
1272 valid_i = coeff_pred(terms_i);
1273 }
1274 else
1275 inc = false;
1276
1277 if (!valid_i || !valid_j || !valid_k)
1278 inc = true;
1279 }
1280
1281 if (cdim == d)
1282 stop = true;
1283
1284 cnt++;
1285 }
1286
1287 Cijk->fillComplete();
1288
1289 return Cijk;
1290 }
1291
1292 template <typename ordinal_type>
1294 ordinal_type i_order, j_order, k_order;
1296 ordinal_type i, j, k;
1297
1298 Cijk_1D_Iterator(ordinal_type p = 0, bool sym = false) :
1299 i_order(p), j_order(p), k_order(p),
1300 symmetric(sym),
1301 i(0), j(0), k(0) {}
1302
1303 Cijk_1D_Iterator(ordinal_type p_i, ordinal_type p_j, ordinal_type p_k,
1304 bool sym) :
1305 i_order(p_i), j_order(p_j), k_order(p_k),
1306 symmetric(sym),
1307 i(0), j(0), k(0) {}
1308
1309 // Reset i,j,k to first non-zero
1310 void reset() { i = 0; j = 0; k = 0; }
1311
1312 // Increment i,j,k to next non-zero with constraint i >= j >= k
1313 // Return false if no more non-zeros
1314 bool increment() {
1315 bool zero = true;
1316
1317 // Increment terms to next non-zero
1318 while (zero) {
1319 bool more = increment_once();
1320 if (!more) return false;
1321 zero = is_zero();
1322 }
1323
1324 return true;
1325 }
1326
1327 // Increment i,j,k to next non-zero with constraint i >= j >= k
1328 // Return false if no more non-zeros
1329 bool increment(ordinal_type& delta_i,
1330 ordinal_type& delta_j,
1331 ordinal_type& delta_k) {
1332 bool zero = true;
1333 bool more = true;
1334 ordinal_type i0 = i;
1335 ordinal_type j0 = j;
1336 ordinal_type k0 = k;
1337
1338 // Increment terms to next non-zero
1339 while (more && zero) {
1340 more = increment_once();
1341 if (more)
1342 zero = is_zero();
1343 }
1344
1345 delta_i = i-i0;
1346 delta_j = j-j0;
1347 delta_k = k-k0;
1348
1349 if (!more)
1350 return false;
1351
1352 return true;
1353 }
1354
1355 private:
1356
1357 // Increment i,j,k to next term with constraint i >= j >= k
1358 // If no more terms, reset to first term and return false
1360 ++i;
1361 if (i > i_order) {
1362 ++j;
1363 if (j > j_order) {
1364 ++k;
1365 if (k > k_order) {
1366 i = 0;
1367 j = 0;
1368 k = 0;
1369 return false;
1370 }
1371 j = 0;
1372 }
1373 i = 0;
1374 }
1375 return true;
1376 }
1377
1378 // Determine if term is zero
1379 bool is_zero() const {
1380 if (k+j < i || i+j < k || i+k < j) return true;
1381 if (symmetric && ((k+j) % 2) != (i % 2) ) return true;
1382 return false;
1383 }
1384 };
1385
1386 template <typename ordinal_type,
1387 typename value_type,
1388 typename basis_set_type,
1389 typename basis_map_type,
1390 typename coeff_predicate_type,
1391 typename k_coeff_predicate_type>
1392 static Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
1394 const Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases,
1395 const basis_set_type& basis_set,
1396 const basis_map_type& basis_map,
1397 const coeff_predicate_type& coeff_pred,
1398 const k_coeff_predicate_type& k_coeff_pred,
1399 bool symmetric = false,
1400 const value_type sparse_tol = 1.0e-12) {
1401#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
1402 TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Time");
1403#endif
1404 typedef typename basis_map_type::value_type coeff_type;
1405 ordinal_type d = bases.size();
1406
1407 // The algorithm for computing Cijk = < \Psi_i \Psi_j \Psi_k > here
1408 // works by factoring
1409 // < \Psi_i \Psi_j \Psi_k > =
1410 // < \psi^1_{i_1}\psi^1_{j_1}\psi^1_{k_1} >_1 x ... x
1411 // < \psi^d_{i_d}\psi^d_{j_d}\psi^d_{k_d} >_d
1412 // We compute the sparse triple product < \psi^l_i\psi^l_j\psi^l_k >_l
1413 // for each dimension l, and then compute all non-zero products of these
1414 // terms. The complexity arises from iterating through all possible
1415 // combinations, throwing out terms that aren't in the basis and are
1416 // beyond the k-order limit provided
1417 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> > Cijk =
1419
1420 // Create 1-D triple products
1421 Teuchos::Array< Teuchos::RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
1422 for (ordinal_type i=0; i<d; i++) {
1423 Cijk_1d[i] =
1424 bases[i]->computeSparseTripleProductTensor(bases[i]->order()+1);
1425 }
1426
1427 // Create i, j, k iterators for each dimension
1428 typedef Cijk_1D_Iterator<ordinal_type> Cijk_Iterator;
1429 Teuchos::Array< Cijk_1D_Iterator<ordinal_type> > Cijk_1d_iterators(d);
1430 coeff_type terms_i(d), terms_j(d), terms_k(d);
1431 for (ordinal_type dim=0; dim<d; dim++) {
1432 Cijk_1d_iterators[dim] = Cijk_Iterator(bases[dim]->order(), symmetric);
1433 terms_i[dim] = 0;
1434 terms_j[dim] = 0;
1435 terms_k[dim] = 0;
1436 }
1437
1438 ordinal_type I = 0;
1439 ordinal_type J = 0;
1440 ordinal_type K = 0;
1441 bool valid_i = coeff_pred(terms_i);
1442 bool valid_j = coeff_pred(terms_j);
1443 bool valid_k = k_coeff_pred(terms_k);
1444 bool stop = !valid_i || !valid_j || !valid_k;
1445 ordinal_type cnt = 0;
1446 while (!stop) {
1447
1448 typename basis_set_type::const_iterator k = basis_set.find(terms_k);
1449 typename basis_set_type::const_iterator i = basis_set.find(terms_i);
1450 typename basis_set_type::const_iterator j = basis_set.find(terms_j);
1451 K = k->second;
1452 I = i->second;
1453 J = j->second;
1454
1455 value_type c = value_type(1.0);
1456 for (ordinal_type dim=0; dim<d; dim++) {
1457 c *= Cijk_1d[dim]->getValue(terms_i[dim], terms_j[dim], terms_k[dim]);
1458 }
1459 if (std::abs(c) > sparse_tol) {
1460 Cijk->add_term(I,J,K,c);
1461 // Cijk->add_term(I,K,J,c);
1462 // Cijk->add_term(J,I,K,c);
1463 // Cijk->add_term(J,K,I,c);
1464 // Cijk->add_term(K,I,J,c);
1465 // Cijk->add_term(K,J,I,c);
1466 }
1467
1468 // Increment iterators to the next valid term
1469 // We keep i >= j >= k for each dimension
1470 ordinal_type cdim = 0;
1471 bool inc = true;
1472 while (inc && cdim < d) {
1473 bool more = Cijk_1d_iterators[cdim].increment();
1474 terms_i[cdim] = Cijk_1d_iterators[cdim].i;
1475 terms_j[cdim] = Cijk_1d_iterators[cdim].j;
1476 terms_k[cdim] = Cijk_1d_iterators[cdim].k;
1477 if (!more) {
1478 ++cdim;
1479 }
1480 else {
1481 valid_i = coeff_pred(terms_i);
1482 valid_j = coeff_pred(terms_j);
1483 valid_k = k_coeff_pred(terms_k);
1484
1485 if (valid_i && valid_j && valid_k)
1486 inc = false;
1487 }
1488 }
1489
1490 if (cdim == d)
1491 stop = true;
1492
1493 cnt++;
1494 }
1495
1496 Cijk->fillComplete();
1497
1498 return Cijk;
1499 }
1500 };
1501
1505 template <typename ordinal_type, typename value_type>
1507 public:
1508
1513 /*
1514 * Returns expansion total order.
1515 * This version is for an isotropic expansion of total order \c p in
1516 * \c d dimensions.
1517 */
1518 static ordinal_type
1519 compute_terms(ordinal_type p, ordinal_type d,
1520 ordinal_type& sz,
1521 Teuchos::Array< MultiIndex<ordinal_type> >& terms,
1522 Teuchos::Array<ordinal_type>& num_terms);
1523
1528 /*
1529 * Returns expansion total order.
1530 * This version allows for anisotropy in the maximum order in each
1531 * dimension.
1532 */
1533 static ordinal_type
1534 compute_terms(const Teuchos::Array<ordinal_type>& basis_orders,
1535 ordinal_type& sz,
1536 Teuchos::Array< MultiIndex<ordinal_type> >& terms,
1537 Teuchos::Array<ordinal_type>& num_terms);
1538
1543 static ordinal_type
1545 const Teuchos::Array< MultiIndex<ordinal_type> >& terms,
1546 const Teuchos::Array<ordinal_type>& num_terms,
1547 ordinal_type max_p);
1548
1549 };
1550
1551}
1552
1553template <typename ordinal_type, typename value_type>
1554ordinal_type
1556compute_terms(ordinal_type p, ordinal_type d,
1557 ordinal_type& sz,
1558 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> >& terms,
1559 Teuchos::Array<ordinal_type>& num_terms)
1560{
1561 Teuchos::Array<ordinal_type> basis_orders(d, p);
1562 return compute_terms(basis_orders, sz, terms, num_terms);
1563}
1564
1565template <typename ordinal_type, typename value_type>
1568compute_terms(const Teuchos::Array<ordinal_type>& basis_orders,
1569 ordinal_type& sz,
1570 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> >& terms,
1571 Teuchos::Array<ordinal_type>& num_terms)
1572{
1573 // The approach here for ordering the terms is inductive on the total
1574 // order p. We get the terms of total order p from the terms of total
1575 // order p-1 by incrementing the orders of the first dimension by 1.
1576 // We then increment the orders of the second dimension by 1 for all of the
1577 // terms whose first dimension order is 0. We then repeat for the third
1578 // dimension whose first and second dimension orders are 0, and so on.
1579 // How this is done is most easily illustrated by an example of dimension 3:
1580 //
1581 // Order terms cnt Order terms cnt
1582 // 0 0 0 0 4 4 0 0 15 5 1
1583 // 3 1 0
1584 // 1 1 0 0 3 2 1 3 0 1
1585 // 0 1 0 2 2 0
1586 // 0 0 1 2 1 1
1587 // 2 0 2
1588 // 2 2 0 0 6 3 1 1 3 0
1589 // 1 1 0 1 2 1
1590 // 1 0 1 1 1 2
1591 // 0 2 0 1 0 3
1592 // 0 1 1 0 4 0
1593 // 0 0 2 0 3 1
1594 // 0 2 2
1595 // 3 3 0 0 10 4 1 0 1 3
1596 // 2 1 0 0 0 4
1597 // 2 0 1
1598 // 1 2 0
1599 // 1 1 1
1600 // 1 0 2
1601 // 0 3 0
1602 // 0 2 1
1603 // 0 1 2
1604 // 0 0 3
1605
1606 // Compute total order
1607 ordinal_type d = basis_orders.size();
1608 ordinal_type p = 0;
1609 for (ordinal_type i=0; i<d; i++) {
1610 if (basis_orders[i] > p)
1611 p = basis_orders[i];
1612 }
1613
1614 // Temporary array of terms grouped in terms of same order
1615 Teuchos::Array< Teuchos::Array< MultiIndex<ordinal_type> > > terms_order(p+1);
1616
1617 // Store number of terms up to each order
1618 num_terms.resize(p+2, ordinal_type(0));
1619
1620 // Set order 0
1621 terms_order[0].resize(1);
1622 terms_order[0][0].resize(d, ordinal_type(0));
1623 num_terms[0] = 1;
1624
1625 // The array "cnt" stores the number of terms we need to increment for each
1626 // dimension.
1627 Teuchos::Array<ordinal_type> cnt(d), cnt_next(d);
1628 MultiIndex<ordinal_type> term(d);
1629 for (ordinal_type j=0; j<d; j++) {
1630 if (basis_orders[j] >= 1)
1631 cnt[j] = 1;
1632 else
1633 cnt[j] = 0;
1634 cnt_next[j] = 0;
1635 }
1636
1637 sz = 1;
1638 // Loop over orders
1639 for (ordinal_type k=1; k<=p; k++) {
1640
1641 num_terms[k] = num_terms[k-1];
1642
1643 // Stores the index of the term we copying
1644 ordinal_type prev = 0;
1645
1646 // Loop over dimensions
1647 for (ordinal_type j=0; j<d; j++) {
1648
1649 // Increment orders of cnt[j] terms for dimension j
1650 for (ordinal_type i=0; i<cnt[j]; i++) {
1651 if (terms_order[k-1][prev+i][j] < basis_orders[j]) {
1652 term = terms_order[k-1][prev+i];
1653 ++term[j];
1654 terms_order[k].push_back(term);
1655 ++sz;
1656 num_terms[k]++;
1657 for (ordinal_type l=0; l<=j; l++)
1658 ++cnt_next[l];
1659 }
1660 }
1661
1662 // Move forward to where all orders for dimension j are 0
1663 if (j < d-1)
1664 prev += cnt[j] - cnt[j+1];
1665
1666 }
1667
1668 // Update the number of terms we must increment for the new order
1669 for (ordinal_type j=0; j<d; j++) {
1670 cnt[j] = cnt_next[j];
1671 cnt_next[j] = 0;
1672 }
1673
1674 }
1675
1676 num_terms[p+1] = sz;
1677
1678 // Copy into final terms array
1679 terms.resize(sz);
1680 ordinal_type i = 0;
1681 for (ordinal_type k=0; k<=p; k++) {
1682 ordinal_type num_k = terms_order[k].size();
1683 for (ordinal_type j=0; j<num_k; j++)
1684 terms[i++] = terms_order[k][j];
1685 }
1686
1687 return p;
1688}
1689
1690template <typename ordinal_type, typename value_type>
1694 const Teuchos::Array< Stokhos::MultiIndex<ordinal_type> >& terms,
1695 const Teuchos::Array<ordinal_type>& num_terms,
1696 ordinal_type max_p)
1697{
1698 // The approach here for computing the index is to find the order block
1699 // corresponding to this term by adding up the component orders. We then
1700 // do a linear search through the terms_order array for this order
1701
1702 // First compute order of term
1703 ordinal_type d = term.dimension();
1704 ordinal_type ord = 0;
1705 for (ordinal_type i=0; i<d; i++)
1706 ord += term[i];
1707 TEUCHOS_TEST_FOR_EXCEPTION(ord < 0 || ord > max_p, std::logic_error,
1708 "Stokhos::CompletePolynomialBasis::compute_index(): " <<
1709 "Term has invalid order " << ord);
1710
1711 // Now search through terms of that order to find a match
1712 ordinal_type k;
1713 if (ord == 0)
1714 k = 0;
1715 else
1716 k = num_terms[ord-1];
1717 ordinal_type k_max=num_terms[ord];
1718 bool found = false;
1719 while (k < k_max && !found) {
1720 bool found_term = true;
1721 for (ordinal_type j=0; j<d; j++) {
1722 found_term = found_term && (term[j] == terms[k][j]);
1723 if (!found_term)
1724 break;
1725 }
1726 found = found_term;
1727 ++k;
1728 }
1729 TEUCHOS_TEST_FOR_EXCEPTION(k >= k_max && !found, std::logic_error,
1730 "Stokhos::CompletePolynomialBasis::compute_index(): " <<
1731 "Could not find specified term.");
1732
1733 return k-1;
1734}
1735
1736#endif
expr val()
Iterator class for iterating over elements of the index set.
multiindex_type index
Current value of iterator.
Iterator(ordinal_type max_order_, const multiindex_type &component_max_order_, const multiindex_type &index_)
Constructor.
multiindex_type component_max_order
Maximum order for each component.
Iterator & operator++()
Prefix increment, i.e., ++iterator.
std::iterator< std::input_iterator_tag, multiindex_type > base_type
bool operator==(const Iterator &it) const
Compare equality of iterators.
Teuchos::Array< ordinal_type > orders
Maximum orders for each term to determine how to increment.
bool operator!=(const Iterator &it) const
Compare inequality of iterators.
Iterator & operator++(int)
Postfix increment, i.e., iterator++.
An anisotropic total order index set.
multiindex_type max_orders() const
Return maximum order for each dimension.
const_iterator end() const
Return iterator for end of the index set.
ordinal_type lower_order
Lower order of index set.
multiindex_type lower
Component-wise lower bounds.
ordinal_type dimension() const
Return dimension.
AnisotropicTotalOrderIndexSet(ordinal_type upper_order_, const multiindex_type &lower_, const multiindex_type &upper_)
Constructor.
AnisotropicTotalOrderIndexSet(ordinal_type upper_order_, const multiindex_type &upper_)
Constructor.
multiindex_type upper
Component-wise upper bounds.
const_iterator begin() const
Return iterator for first element in the set.
ordinal_type upper_order
Upper order of index set.
Utilities for indexing a multi-variate complete polynomial basis.
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...
static ordinal_type compute_terms(const Teuchos::Array< ordinal_type > &basis_orders, 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...
static ordinal_type compute_index(const MultiIndex< ordinal_type > &term, const Teuchos::Array< MultiIndex< ordinal_type > > &terms, const Teuchos::Array< ordinal_type > &num_terms, ordinal_type max_p)
Compute basis index given the orders for each basis dimension.
A functor for comparing floating-point numbers to some tolerance.
FloatingPointLess(const value_type &tol_=1.0e-12)
Constructor.
bool operator()(const value_type &a, const value_type &b) const
Compare if a < b.
A comparison functor implementing a strict weak ordering based lexographic ordering.
bool operator()(const term_type &a, const term_type &b) const
Determine if a is less than b.
LexographicLess(const compare_type &cmp_=compare_type())
Constructor.
compare_type cmp
Element comparison functor.
bool equal(const element_type &a, const element_type &b) const
Determine if two elements a and b are equal.
term_type::ordinal_type ordinal_type
term_type::element_type element_type
A comparison functor implementing a strict weak ordering based Morton Z-ordering.
term_type::element_type element_type
term_type::ordinal_type ordinal_type
bool operator()(const term_type &a, const term_type &b) const
A multidimensional index.
void resize(ordinal_type d, ordinal_type v=ordinal_type(0))
Resize.
Teuchos::Array< ordinal_type > index
index terms
Teuchos::Array< element_type > & getTerm()
Term access.
ordinal_type dimension() const
Dimension.
const Teuchos::Array< element_type > & getTerm() const
Term access.
bool operator==(const MultiIndex &idx) const
Compare equality.
MultiIndex & termWiseMin(const MultiIndex &idx)
Replace multiindex with min of this and other multiindex.
bool operator!=(const MultiIndex &idx) const
Compare equality.
MultiIndex(ordinal_type dim, ordinal_type v=ordinal_type(0))
Constructor.
ordinal_type order() const
Compute total order of index.
std::ostream & print(std::ostream &os) const
Print multiindex.
MultiIndex & termWiseMin(const ordinal_type idx)
Replace multiindex with min of this and given value.
ordinal_type size() const
Size.
void init(ordinal_type v)
Initialize.
MultiIndex & termWiseMax(const MultiIndex &idx)
Replace multiindex with max of this and other multiindex.
bool termWiseLEQ(const MultiIndex &idx) const
Compare term-wise less-than or equal-to.
MultiIndex & termWiseMax(const ordinal_type idx)
Replace multiindex with max of this and given value.
const ordinal_type & operator[](ordinal_type i) const
Term access.
Abstract base class for 1-D orthogonal polynomials.
Utilities for indexing a multi-variate complete polynomial basis.
static void buildProductBasis(const index_set_type &index_set, basis_set_type &basis_set, basis_map_type &basis_map)
Generate a product basis from an index set.
static Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorNew(const Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > &bases, const basis_set_type &basis_set, const basis_map_type &basis_map, const coeff_predicate_type &coeff_pred, const k_coeff_predicate_type &k_coeff_pred, bool symmetric=false, const value_type sparse_tol=1.0e-12)
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.
static Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor(const Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > &bases, const basis_set_type &basis_set, const basis_map_type &basis_map, const coeff_predicate_type &coeff_pred, const k_coeff_predicate_type &k_coeff_pred, const value_type sparse_tol=1.0e-12)
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.
element_type order() const
Compute total order of tensor product element.
std::ostream & print(std::ostream &os) const
Print multiindex.
Teuchos::Array< element_type > & getTerm()
Term access.
ordinal_type dimension() const
Return dimension.
ordinal_type size() const
Return dimension.
const Teuchos::Array< element_type > & getTerm() const
Term access.
const element_type & operator[](ordinal_type i) const
Term access.
TensorProductElement(ordinal_type dim, const element_type &val=element_type(0))
Constructor.
Teuchos::Array< element_type > term
Array storing term elements.
Iterator class for iterating over elements of the index set.
multiindex_type index
Current value of iterator.
std::iterator< std::input_iterator_tag, multiindex_type > base_type
Iterator & operator++(int)
Postfix increment, i.e., iterator++.
Iterator(const multiindex_type &upper_, const multiindex_type &index_)
Constructor.
multiindex_type upper
Upper bound of iterator.
const_reference operator*() const
Dereference.
bool operator!=(const Iterator &it) const
Compare inequality of iterators.
Iterator & operator++()
Prefix increment, i.e., ++iterator.
bool operator==(const Iterator &it) const
Compare equality of iterators.
multiindex_type upper
Upper bound of index set.
TensorProductIndexSet(const multiindex_type &lower_, const multiindex_type &upper_)
Constructor.
ordinal_type dimension() const
Return dimension.
TensorProductIndexSet(ordinal_type dim_, ordinal_type upper_)
Constructor.
TensorProductIndexSet(ordinal_type dim_, ordinal_type lower_, ordinal_type upper_)
Constructor.
const_iterator begin() const
Return iterator for first element in the set.
TensorProductIndexSet(const multiindex_type &upper_)
Constructor.
const_iterator end() const
Return iterator for end of the index set.
MultiIndex< ordinal_type > multiindex_type
multiindex_type lower
Lower bound of index set.
multiindex_type max_orders() const
Return maximum order for each dimension.
Iterator class for iterating over elements of the index set.
bool operator==(const Iterator &it) const
Compare equality of iterators.
std::iterator< std::input_iterator_tag, multiindex_type > base_type
Iterator(ordinal_type max_order_, const multiindex_type &index_)
Constructor.
Iterator & operator++()
Prefix increment, i.e., ++iterator.
const_pointer operator->() const
Dereference.
multiindex_type index
Current value of iterator.
Teuchos::Array< ordinal_type > orders
Maximum orders for each term to determine how to increment.
bool operator!=(const Iterator &it) const
Compare inequality of iterators.
const_reference operator*() const
Dereference.
ordinal_type max_order
Maximum order of iterator.
Iterator & operator++(int)
Postfix increment, i.e., iterator++.
An isotropic total order index set.
const_iterator end() const
Return iterator for end of the index set.
const_iterator begin() const
Return iterator for first element in the set.
MultiIndex< ordinal_type > multiindex_type
ordinal_type upper
Upper order of index set.
multiindex_type max_orders() const
Return maximum order for each dimension.
ordinal_type lower
Lower order of index set.
ordinal_type dimension() const
Return dimension.
TotalOrderIndexSet(ordinal_type dim_, ordinal_type lower_, ordinal_type upper_)
Constructor.
TotalOrderIndexSet(ordinal_type dim_, ordinal_type upper_)
Constructor.
A comparison functor implementing a strict weak ordering based total-order ordering,...
term_type::element_type element_type
term_type::ordinal_type ordinal_type
TotalOrderLess(const compare_type &cmp_=compare_type())
Constructor.
bool operator()(const term_type &a, const term_type &b) const
bool equal(const element_type &a, const element_type &b) const
Determine if two elements a and b are equal.
compare_type cmp
Element comparison functor.
Top-level namespace for Stokhos classes and functions.
ordinal_type totalOrderMapping(const Stokhos::MultiIndex< ordinal_type > &index)
ordinal_type lexicographicMapping(const Stokhos::MultiIndex< ordinal_type > &index, ordinal_type max_order)
std::ostream & operator<<(std::ostream &os, const ProductContainer< coeff_type > &vec)
ordinal_type n_choose_k(const ordinal_type &n, const ordinal_type &k)
Compute bionomial coefficient (n ; k) = n!/( k! (n-k)! )
bool increment(ordinal_type &delta_i, ordinal_type &delta_j, ordinal_type &delta_k)
Cijk_1D_Iterator(ordinal_type p_i, ordinal_type p_j, ordinal_type p_k, bool sym)
Predicate functor for building sparse triple products.
bool operator()(const coeff_type &term) const
TensorProductPredicate(const coeff_type &orders_)
Predicate functor for building sparse triple products based on total order.
MultiIndex< ordinal_type > coeff_type
bool operator()(const coeff_type &term) const
TotalOrderPredicate(ordinal_type max_order_, const coeff_type &orders_)