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(
"a_size", &a_size,
"Size of a (matrix) partition");
72 CLP.setOption(
"x_size", &x_size,
"Size of x (input vector) partition");
74 CLP.setOption(
"y_size", &y_size,
"Size of y (output vector) partition");
75 bool save_3tensor =
false;
76 CLP.setOption(
"save_3tensor",
"no-save_3tensor", &save_3tensor,
77 "Save full 3tensor to file");
78 std::string file_3tensor =
"Cijk.dat";
79 CLP.setOption(
"filename_3tensor", &file_3tensor,
80 "Filename to store full 3-tensor");
86 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
87 const double alpha = 1.0;
88 const double beta = symmetric ? 1.0 : 2.0 ;
89 for (
int i=0; i<d; i++) {
91 p, alpha, beta,
true));
95 RCP<const basis_type> basis = Teuchos::rcp(
new basis_type(bases, drop));
99 typedef Cijk_LTB_type::CijkNode
node_type;
100 Teuchos::RCP<Cijk_LTB_type> Cijk =
103 int sz = basis->size();
104 std::cout <<
"basis size = " << sz
105 <<
" num nonzero Cijk entries = " << Cijk->num_entries()
109 if (a_size > sz) a_size = sz;
110 if (x_size > sz) x_size = sz;
111 if (y_size > sz) y_size = sz;
113 Teuchos::Array< Teuchos::RCP<const node_type> > node_stack;
114 Teuchos::Array< int > index_stack;
115 node_stack.push_back(Cijk->getHeadNode());
116 index_stack.push_back(0);
117 Teuchos::RCP<const node_type> node;
119 Teuchos::Array< Teuchos::RCP<const node_type> > partition_stack;
120 while (node_stack.size() > 0) {
121 node = node_stack.back();
122 child_index = index_stack.back();
126 partition_stack.push_back(node);
127 node_stack.pop_back();
128 index_stack.pop_back();
132 else if (node->i_size <= y_size &&
133 node->j_size <= a_size &&
134 node->k_size <= x_size) {
135 partition_stack.push_back(node);
136 node_stack.pop_back();
137 index_stack.pop_back();
141 else if (child_index < node->children.size()) {
142 ++index_stack.back();
143 node = node->children[child_index];
144 node_stack.push_back(node);
145 index_stack.push_back(0);
150 node_stack.pop_back();
151 index_stack.pop_back();
156 std::cout <<
"num partitions = " << partition_stack.size() << std::endl;
173 Teuchos::Array< Teuchos::Array<int> > tuples;
174 for (
int part=0; part<partition_stack.size(); ++part) {
175 node = partition_stack[part];
176 node_stack.push_back(node);
177 index_stack.push_back(0);
178 while (node_stack.size() > 0) {
179 node = node_stack.back();
180 child_index = index_stack.back();
184 Cijk_Iterator cijk_iterator(node->p_i,
190 Teuchos::Array<int> t(4);
191 int I = node->i_begin + cijk_iterator.i;
192 int J = node->j_begin + cijk_iterator.j;
193 int K = node->k_begin + cijk_iterator.k;
199 more = cijk_iterator.increment();
201 node_stack.pop_back();
202 index_stack.pop_back();
206 else if (child_index < node->children.size()) {
207 ++index_stack.back();
208 node = node->children[child_index];
209 node_stack.push_back(node);
210 index_stack.push_back(0);
215 node_stack.pop_back();
216 index_stack.pop_back();
224 std::ofstream cijk_file(file_3tensor.c_str());
225 cijk_file.precision(14);
226 cijk_file.setf(std::ios::scientific);
227 cijk_file <<
"i, j, k, part" << std::endl;
228 for (
int i=0; i<tuples.size(); ++i) {
229 cijk_file << tuples[i][0] <<
", " 230 << tuples[i][1] <<
", " 231 << tuples[i][2] <<
", " 232 << tuples[i][3] << std::endl;
238 catch (std::exception& e) {
239 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.