45#include "Teuchos_CommandLineProcessor.hpp"
46#include "Teuchos_Array.hpp"
47#include "Teuchos_RCP.hpp"
61 Teuchos::CommandLineProcessor
CLP;
63 "This example generates partitions the Cijk tensor for a lexicographic tree basis.\n");
65 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
67 CLP.setOption(
"order", &p,
"Polynomial order");
68 double drop = 1.0e-12;
69 CLP.setOption(
"drop", &drop,
"Drop tolerance");
70 bool symmetric =
true;
71 CLP.setOption(
"symmetric",
"asymmetric", &symmetric,
"Use basis polynomials with symmetric PDF");
73 CLP.setOption(
"level", &level,
"Level to partition");
74 bool save_3tensor =
false;
75 CLP.setOption(
"save_3tensor",
"no-save_3tensor", &save_3tensor,
76 "Save full 3tensor to file");
77 std::string file_3tensor =
"Cijk.dat";
78 CLP.setOption(
"filename_3tensor", &file_3tensor,
79 "Filename to store full 3-tensor");
85 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
86 const double alpha = 1.0;
87 const double beta = symmetric ? 1.0 : 2.0 ;
88 for (
int i=0; i<d; i++) {
90 p, alpha, beta,
true));
94 RCP<const basis_type> basis = Teuchos::rcp(
new basis_type(bases, drop));
98 typedef Cijk_LTB_type::CijkNode
node_type;
99 Teuchos::RCP<Cijk_LTB_type> Cijk =
100 computeTripleProductTensorLTB(*basis, symmetric);
102 int sz = basis->size();
103 std::cout <<
"basis size = " << sz
104 <<
" num nonzero Cijk entries = " << Cijk->num_entries()
108 Teuchos::Array< Teuchos::RCP<const node_type> > node_stack;
109 Teuchos::Array< int > index_stack;
110 node_stack.push_back(Cijk->getHeadNode());
111 index_stack.push_back(0);
112 Teuchos::RCP<const node_type> node;
114 Teuchos::Array< Teuchos::RCP<const node_type> > partition_stack;
116 while (node_stack.size() > 0) {
117 node = node_stack.back();
118 child_index = index_stack.back();
122 partition_stack.push_back(node);
123 node_stack.pop_back();
124 index_stack.pop_back();
129 else if (my_level == level) {
130 partition_stack.push_back(node);
131 node_stack.pop_back();
132 index_stack.pop_back();
137 else if (child_index < node->children.size()) {
138 ++index_stack.back();
139 node = node->children[child_index];
140 node_stack.push_back(node);
141 index_stack.push_back(0);
147 node_stack.pop_back();
148 index_stack.pop_back();
155 int max_i_size = 0, max_j_size = 0, max_k_size = 0;
156 for (
int part=0; part<partition_stack.size(); ++part) {
157 node = partition_stack[part];
158 if (node->i_size > max_i_size) max_i_size = node->i_size;
159 if (node->j_size > max_j_size) max_j_size = node->j_size;
160 if (node->k_size > max_k_size) max_k_size = node->k_size;
162 std::cout <<
"num partitions = " << partition_stack.size() << std::endl
163 <<
"max i size = " << max_i_size << std::endl
164 <<
"max j size = " << max_j_size << std::endl
165 <<
"max k size = " << max_k_size << std::endl;
169 Teuchos::Array< Teuchos::Array<int> > tuples;
170 for (
int part=0; part<partition_stack.size(); ++part) {
171 node = partition_stack[part];
172 node_stack.push_back(node);
173 index_stack.push_back(0);
174 while (node_stack.size() > 0) {
175 node = node_stack.back();
176 child_index = index_stack.back();
180 Cijk_Iterator cijk_iterator(node->p_i,
186 Teuchos::Array<int> t(4);
187 int I = node->i_begin + cijk_iterator.i;
188 int J = node->j_begin + cijk_iterator.j;
189 int K = node->k_begin + cijk_iterator.k;
195 more = cijk_iterator.increment();
197 node_stack.pop_back();
198 index_stack.pop_back();
202 else if (child_index < node->children.size()) {
203 ++index_stack.back();
204 node = node->children[child_index];
205 node_stack.push_back(node);
206 index_stack.push_back(0);
211 node_stack.pop_back();
212 index_stack.pop_back();
220 std::ofstream cijk_file(file_3tensor.c_str());
221 cijk_file.precision(14);
222 cijk_file.setf(std::ios::scientific);
223 cijk_file <<
"i, j, k, part" << std::endl;
224 for (
int i=0; i<tuples.size(); ++i) {
225 cijk_file << tuples[i][0] <<
", "
226 << tuples[i][1] <<
", "
227 << tuples[i][2] <<
", "
228 << tuples[i][3] << std::endl;
234 catch (std::exception& e) {
235 std::cout << e.what() << std::endl;
int main(int argc, char **argv)
Data structure storing a sparse 3-tensor C(i,j,k) in a a tree-based format for lexicographically orde...
A comparison functor implementing a strict weak ordering based lexographic ordering.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Stokhos::LegendreBasis< int, double > basis_type