Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
cijk_partition_rcb.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 
47 // Growth policies
48 const int num_growth_types = 2;
51 const char *growth_type_names[] = { "slow", "moderate" };
52 
53 // Product Basis types
55 const int num_prod_basis_types = 4;
58 const char *prod_basis_type_names[] = {
59  "complete", "tensor", "total", "smolyak" };
60 
61 // Ordering types
63 const int num_ordering_types = 2;
66 const char *ordering_type_names[] = {
67  "total", "lexicographic" };
68 
69 // Symmetry types
70 const int num_symmetry_types = 3;
75 };
76 const char *symmetry_type_names[] = {
77  "none", "2-way", "6-way" };
78 
79 using Teuchos::rcp;
80 using Teuchos::RCP;
81 using Teuchos::Array;
82 using Teuchos::ArrayView;
83 using Teuchos::ArrayRCP;
84 
85 int main(int argc, char **argv)
86 {
87  try {
88 
89  // Setup command line options
90  Teuchos::CommandLineProcessor CLP;
91  CLP.setDocString(
92  "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
93  int d = 5;
94  CLP.setOption("dimension", &d, "Stochastic dimension");
95  int p = 3;
96  CLP.setOption("order", &p, "Polynomial order");
97  double drop = 1.0e-12;
98  CLP.setOption("drop", &drop, "Drop tolerance");
99  bool symmetric = true;
100  CLP.setOption("symmetric", "asymmetric", &symmetric, "Use basis polynomials with symmetric PDF");
102  CLP.setOption("growth", &growth_type,
104  "Growth type");
105  ProductBasisType prod_basis_type = TOTAL;
106  CLP.setOption("product_basis", &prod_basis_type,
109  "Product basis type");
110  OrderingType ordering_type = LEXICOGRAPHIC_ORDERING;
111  CLP.setOption("ordering", &ordering_type,
114  "Product basis ordering");
116  CLP.setOption("symmetry", &symmetry_type,
119  "Cijk symmetry type");
120  bool full = true;
121  CLP.setOption("full", "linear", &full, "Use full or linear expansion");
122  int tile_size = 32;
123  CLP.setOption("tile_size", &tile_size, "Tile size");
124  bool save_3tensor = false;
125  CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
126  "Save full 3tensor to file");
127  std::string file_3tensor = "Cijk.dat";
128  CLP.setOption("filename_3tensor", &file_3tensor,
129  "Filename to store full 3-tensor");
130 
131  // Parse arguments
132  CLP.parse( argc, argv );
133 
134  // Basis
135  Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
136  const double alpha = 1.0;
137  const double beta = symmetric ? 1.0 : 2.0 ;
138  for (int i=0; i<d; i++) {
139  bases[i] = rcp(new Stokhos::JacobiBasis<int,double>(
140  p, alpha, beta, true, growth_type));
141  }
142  RCP<const Stokhos::ProductBasis<int,double> > basis;
145  if (prod_basis_type == COMPLETE)
146  basis =
148  bases, drop));
149  else if (prod_basis_type == TENSOR) {
150  if (ordering_type == TOTAL_ORDERING)
151  basis =
153  bases, drop));
154  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
155  basis =
157  bases, drop));
158  }
159  else if (prod_basis_type == TOTAL) {
160  if (ordering_type == TOTAL_ORDERING)
161  basis =
163  bases, drop));
164  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
165  basis =
167  bases, drop));
168  }
169  else if (prod_basis_type == SMOLYAK) {
170  Stokhos::TotalOrderIndexSet<int> index_set(d, p);
171  if (ordering_type == TOTAL_ORDERING)
172  basis =
174  bases, index_set, drop));
175  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
176  basis =
178  bases, index_set, drop));
179  }
180 
181  // Triple product tensor
183  RCP<Cijk_type> Cijk;
184  if (full)
185  Cijk = basis->computeTripleProductTensor();
186  else
187  Cijk = basis->computeLinearTripleProductTensor();
188 
189  int basis_size = basis->size();
190  std::cout << "basis size = " << basis_size
191  << " num nonzero Cijk entries = " << Cijk->num_entries()
192  << std::endl;
193 
194  // File for saving sparse Cijk tensor and parts
195  std::ofstream cijk_file;
196  if (save_3tensor) {
197  cijk_file.open(file_3tensor.c_str());
198  cijk_file.precision(14);
199  cijk_file.setf(std::ios::scientific);
200  cijk_file << "i, j, k, part" << std::endl;
201  }
202 
203  // Build tensor data list
204  Teuchos::ArrayRCP< Stokhos::CijkData<int,double> > coordinate_list =
205  Stokhos::build_cijk_coordinate_list(*Cijk, symmetry_type);
206 
207  // Partition via RCB
209  rcb_type rcb(tile_size, 10000, coordinate_list());
210  int num_parts = rcb.get_num_parts();
211  std::cout << "num parts = " << num_parts << std::endl;
212 
213  // Print part sizes
214  RCP< Array< RCP<rcb_type::Box> > > parts = rcb.get_parts();
215  for (int i=0; i<num_parts; ++i) {
216  RCP<rcb_type::Box> box = (*parts)[i];
217  std::cout << "part " << i << " bounding box ="
218  << " [ " << box->delta_x << ", " << box->delta_y << ", "
219  << box->delta_z << " ]" << " nnz = "
220  << box->coords.size() << std::endl;
221  }
222 
223  if (save_3tensor) {
224  ArrayRCP<int> part_ids = rcb.get_part_IDs();
225  int idx = 0;
226  Cijk_type::k_iterator k_begin = Cijk->k_begin();
227  Cijk_type::k_iterator k_end = Cijk->k_end();
228  for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
229  int k = index(k_it);
230  Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
231  Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
232  for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
233  int j = index(j_it);
234  Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
235  Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
236  for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
237  int i = index(i_it);
238  if ((symmetry_type == Stokhos::CIJK_NO_SYMMETRY) ||
239  (symmetry_type == Stokhos::CIJK_TWO_WAY_SYMMETRY && j >= k) ||
240  (symmetry_type == Stokhos::CIJK_SIX_WAY_SYMMETRY && i >= j && j >= k)) {
241  cijk_file << i << ", " << j << ", " << k << ", "
242  << part_ids[idx++] << std::endl;
243  }
244  }
245  }
246  }
247  cijk_file.close();
248  }
249 
250  //Teuchos::TimeMonitor::summarize(std::cout);
251 
252  }
253  catch (std::exception& e) {
254  std::cout << e.what() << std::endl;
255  }
256 
257  return 0;
258 }
SparseArrayIterator< index_iterator, value_iterator >::value_type index(const SparseArrayIterator< index_iterator, value_iterator > &it)
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
const char * ordering_type_names[]
const char * symmetry_type_names[]
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
const OrderingType ordering_type_values[]
const ProductBasisType prod_basis_type_values[]
OrderingType
const Stokhos::CijkSymmetryType symmetry_type_values[]
ProductBasisType
const int num_symmetry_types
const Stokhos::GrowthPolicy growth_type_values[]
Jacobi polynomial basis.
const int num_ordering_types
ordinal_type num_entries() const
Return number of non-zero entries.
expr expr expr expr j
const int num_growth_types
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Multivariate orthogonal polynomial basis generated from a Smolyak sparse grid.
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials...
Stokhos::Sparse3Tensor< int, double > Cijk_type
An isotropic total order index set.
const char * prod_basis_type_names[]
int main(int argc, char **argv)
A comparison functor implementing a strict weak ordering based lexographic ordering.
const char * growth_type_names[]
Teuchos::ArrayRCP< CijkData< ordinal_type, scalar_type > > build_cijk_coordinate_list(const Sparse3Tensor< ordinal_type, scalar_type > &Cijk, CijkSymmetryType symmetry_type)
const int num_prod_basis_types