Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
cijk_ltb_partition_level.cpp
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 #include "Stokhos.hpp"
44 
45 #include "Teuchos_CommandLineProcessor.hpp"
46 #include "Teuchos_Array.hpp"
47 #include "Teuchos_RCP.hpp"
48 
49 using Teuchos::Array;
50 using Teuchos::RCP;
51 using Teuchos::rcp;
52 
53 int main(int argc, char **argv)
54 {
55  try {
56 
57  // Setup command line options
58  Teuchos::CommandLineProcessor CLP;
59  CLP.setDocString(
60  "This example generates partitions the Cijk tensor for a lexicographic tree basis.\n");
61  int d = 3;
62  CLP.setOption("dimension", &d, "Stochastic dimension");
63  int p = 5;
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");
69  int level = 1;
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");
77 
78  // Parse arguments
79  CLP.parse( argc, argv );
80 
81  // Basis
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++) {
86  bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<int,double>(
87  p, alpha, beta, true));
88  }
91  RCP<const basis_type> basis = Teuchos::rcp(new basis_type(bases, drop));
92 
93  // Build LTB Cijk
94  typedef Stokhos::LTBSparse3Tensor<int,double> Cijk_LTB_type;
95  typedef Cijk_LTB_type::CijkNode node_type;
96  Teuchos::RCP<Cijk_LTB_type> Cijk =
97  computeTripleProductTensorLTB(*basis, symmetric);
98 
99  int sz = basis->size();
100  std::cout << "basis size = " << sz
101  << " num nonzero Cijk entries = " << Cijk->num_entries()
102  << std::endl;
103 
104  // Setup partitions
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;
110  int child_index;
111  Teuchos::Array< Teuchos::RCP<const node_type> > partition_stack;
112  int my_level = 0;
113  while (node_stack.size() > 0) {
114  node = node_stack.back();
115  child_index = index_stack.back();
116 
117  // Leaf -- If we got here, just push this node into the partitions
118  if (node->is_leaf) {
119  partition_stack.push_back(node);
120  node_stack.pop_back();
121  index_stack.pop_back();
122  --my_level;
123  }
124 
125  // Put nodes into partition if level matches
126  else if (my_level == level) {
127  partition_stack.push_back(node);
128  node_stack.pop_back();
129  index_stack.pop_back();
130  --my_level;
131  }
132 
133  // More children to process -- process them first
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);
139  ++my_level;
140  }
141 
142  // No more children
143  else {
144  node_stack.pop_back();
145  index_stack.pop_back();
146  --my_level;
147  }
148 
149  }
150 
151  // Print statistics
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;
158  }
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;
163 
164  // Build flat list of (i,j,k,part) tuples
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();
174 
175  // Leaf -- store (i,j,k,part) tuples
176  if (node->is_leaf) {
177  Cijk_Iterator cijk_iterator(node->p_i,
178  node->p_j,
179  node->p_k,
180  symmetric);
181  bool more = true;
182  while (more) {
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;
187  t[0] = I;
188  t[1] = J;
189  t[2] = K;
190  t[3] = part;
191  tuples.push_back(t);
192  more = cijk_iterator.increment();
193  }
194  node_stack.pop_back();
195  index_stack.pop_back();
196  }
197 
198  // More children to process -- process them first
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);
204  }
205 
206  // No more children
207  else {
208  node_stack.pop_back();
209  index_stack.pop_back();
210  }
211 
212  }
213  }
214 
215  // Print full 3-tensor to file
216  if (save_3tensor) {
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;
226  }
227  cijk_file.close();
228  }
229 
230  }
231  catch (std::exception& e) {
232  std::cout << e.what() << std::endl;
233  }
234 
235  return 0;
236 }
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
Kokkos::Serial node_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)
Jacobi polynomial basis.
A comparison functor implementing a strict weak ordering based lexographic ordering.