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