Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
cijk_ltb_partition.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 a_size = 100;
70  CLP.setOption("a_size", &a_size, "Size of a (matrix) partition");
71  int x_size = 100;
72  CLP.setOption("x_size", &x_size, "Size of x (input vector) partition");
73  int y_size = 100;
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");
81 
82  // Parse arguments
83  CLP.parse( argc, argv );
84 
85  // Basis
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++) {
90  bases[i] = Teuchos::rcp(new Stokhos::JacobiBasis<int,double>(
91  p, alpha, beta, true));
92  }
95  RCP<const basis_type> basis = Teuchos::rcp(new basis_type(bases, drop));
96 
97  // Build LTB Cijk
98  typedef Stokhos::LTBSparse3Tensor<int,double> Cijk_LTB_type;
99  typedef Cijk_LTB_type::CijkNode node_type;
100  Teuchos::RCP<Cijk_LTB_type> Cijk =
101  computeTripleProductTensorLTB(*basis, symmetric);
102 
103  int sz = basis->size();
104  std::cout << "basis size = " << sz
105  << " num nonzero Cijk entries = " << Cijk->num_entries()
106  << std::endl;
107 
108  // Setup partitions
109  if (a_size > sz) a_size = sz;
110  if (x_size > sz) x_size = sz;
111  if (y_size > sz) y_size = sz;
112 
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;
118  int child_index;
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();
123 
124  // Leaf -- If we got here, just push this node into the partitions
125  if (node->is_leaf) {
126  partition_stack.push_back(node);
127  node_stack.pop_back();
128  index_stack.pop_back();
129  }
130 
131  // Once sizes are small enough, push node onto partition stack
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();
138  }
139 
140  // More children to process -- process them first
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);
146  }
147 
148  // No more children
149  else {
150  node_stack.pop_back();
151  index_stack.pop_back();
152  }
153 
154  }
155 
156  std::cout << "num partitions = " << partition_stack.size() << std::endl;
157  /*
158  for (int part=0; part<partition_stack.size(); ++part) {
159  node = partition_stack[part];
160  std::cout << "partition " << part << ":" << std::endl
161  << "\ti-range: [" << node->i_begin << ","
162  << node->i_begin+node->i_size << ")" << std::endl
163  << "\tj-range: [" << node->j_begin << ","
164  << node->j_begin+node->j_size << ")" << std::endl
165  << "\tk-range: [" << node->k_begin << ","
166  << node->k_begin+node->k_size << ")" << std::endl
167  << "\tnum_non_zeros = " << node->total_num_entries << std::endl;
168  }
169  */
170 
171  // Build flat list of (i,j,k,part) tuples
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();
181 
182  // Leaf -- store (i,j,k,part) tuples
183  if (node->is_leaf) {
184  Cijk_Iterator cijk_iterator(node->p_i,
185  node->p_j,
186  node->p_k,
187  symmetric);
188  bool more = true;
189  while (more) {
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;
194  t[0] = I;
195  t[1] = J;
196  t[2] = K;
197  t[3] = part;
198  tuples.push_back(t);
199  more = cijk_iterator.increment();
200  }
201  node_stack.pop_back();
202  index_stack.pop_back();
203  }
204 
205  // More children to process -- process them first
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);
211  }
212 
213  // No more children
214  else {
215  node_stack.pop_back();
216  index_stack.pop_back();
217  }
218 
219  }
220  }
221 
222  // Print full 3-tensor to file
223  if (save_3tensor) {
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;
233  }
234  cijk_file.close();
235  }
236 
237  }
238  catch (std::exception& e) {
239  std::cout << e.what() << std::endl;
240  }
241 
242  return 0;
243 }
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...
Jacobi polynomial basis.
int main(int argc, char **argv)
A comparison functor implementing a strict weak ordering based lexographic ordering.