45 #include "Teuchos_CommandLineProcessor.hpp" 46 #include "Teuchos_Array.hpp" 47 #include "Teuchos_RCP.hpp" 58 Teuchos::CommandLineProcessor
CLP;
60 "This example generates partitions the Cijk tensor for a lexicographic tree basis.\n");
62 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
64 CLP.setOption(
"order", &p,
"Polynomial order");
65 double drop = 1.0e-12;
66 CLP.setOption(
"drop", &drop,
"Drop tolerance");
67 bool symmetric =
true;
68 CLP.setOption(
"symmetric",
"asymmetric", &symmetric,
"Use basis polynomials with symmetric PDF");
70 CLP.setOption(
"level", &level,
"Level to partition");
71 bool save_3tensor =
false;
72 CLP.setOption(
"save_3tensor",
"no-save_3tensor", &save_3tensor,
73 "Save full 3tensor to file");
74 std::string file_3tensor =
"Cijk.dat";
75 CLP.setOption(
"filename_3tensor", &file_3tensor,
76 "Filename to store full 3-tensor");
82 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
83 const double alpha = 1.0;
84 const double beta = symmetric ? 1.0 : 2.0 ;
85 for (
int i=0; i<d; i++) {
87 p, alpha, beta,
true));
91 RCP<const basis_type> basis = Teuchos::rcp(
new basis_type(bases, drop));
95 typedef Cijk_LTB_type::CijkNode
node_type;
96 Teuchos::RCP<Cijk_LTB_type> Cijk =
99 int sz = basis->size();
100 std::cout <<
"basis size = " << sz
101 <<
" num nonzero Cijk entries = " << Cijk->num_entries()
105 Teuchos::Array< Teuchos::RCP<const node_type> > node_stack;
106 Teuchos::Array< int > index_stack;
107 node_stack.push_back(Cijk->getHeadNode());
108 index_stack.push_back(0);
109 Teuchos::RCP<const node_type> node;
111 Teuchos::Array< Teuchos::RCP<const node_type> > partition_stack;
113 while (node_stack.size() > 0) {
114 node = node_stack.back();
115 child_index = index_stack.back();
119 partition_stack.push_back(node);
120 node_stack.pop_back();
121 index_stack.pop_back();
126 else if (my_level == level) {
127 partition_stack.push_back(node);
128 node_stack.pop_back();
129 index_stack.pop_back();
134 else if (child_index < node->children.size()) {
135 ++index_stack.back();
136 node = node->children[child_index];
137 node_stack.push_back(node);
138 index_stack.push_back(0);
144 node_stack.pop_back();
145 index_stack.pop_back();
152 int max_i_size = 0, max_j_size = 0, max_k_size = 0;
153 for (
int part=0; part<partition_stack.size(); ++part) {
154 node = partition_stack[part];
155 if (node->i_size > max_i_size) max_i_size = node->i_size;
156 if (node->j_size > max_j_size) max_j_size = node->j_size;
157 if (node->k_size > max_k_size) max_k_size = node->k_size;
159 std::cout <<
"num partitions = " << partition_stack.size() << std::endl
160 <<
"max i size = " << max_i_size << std::endl
161 <<
"max j size = " << max_j_size << std::endl
162 <<
"max k size = " << max_k_size << std::endl;
166 Teuchos::Array< Teuchos::Array<int> > tuples;
167 for (
int part=0; part<partition_stack.size(); ++part) {
168 node = partition_stack[part];
169 node_stack.push_back(node);
170 index_stack.push_back(0);
171 while (node_stack.size() > 0) {
172 node = node_stack.back();
173 child_index = index_stack.back();
177 Cijk_Iterator cijk_iterator(node->p_i,
183 Teuchos::Array<int> t(4);
184 int I = node->i_begin + cijk_iterator.i;
185 int J = node->j_begin + cijk_iterator.j;
186 int K = node->k_begin + cijk_iterator.k;
192 more = cijk_iterator.increment();
194 node_stack.pop_back();
195 index_stack.pop_back();
199 else if (child_index < node->children.size()) {
200 ++index_stack.back();
201 node = node->children[child_index];
202 node_stack.push_back(node);
203 index_stack.push_back(0);
208 node_stack.pop_back();
209 index_stack.pop_back();
217 std::ofstream cijk_file(file_3tensor.c_str());
218 cijk_file.precision(14);
219 cijk_file.setf(std::ios::scientific);
220 cijk_file <<
"i, j, k, part" << std::endl;
221 for (
int i=0; i<tuples.size(); ++i) {
222 cijk_file << tuples[i][0] <<
", " 223 << tuples[i][1] <<
", " 224 << tuples[i][2] <<
", " 225 << tuples[i][3] << std::endl;
231 catch (std::exception& e) {
232 std::cout << e.what() << std::endl;
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTB(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...
Stokhos::LegendreBasis< int, double > basis_type
Data structure storing a sparse 3-tensor C(i,j,k) in a a tree-based format for lexicographically orde...
int main(int argc, char **argv)
A comparison functor implementing a strict weak ordering based lexographic ordering.