Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
cijk_simple_tile.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"
43#include "Teuchos_CommandLineProcessor.hpp"
44#include "Teuchos_ParameterList.hpp"
45
46#include <algorithm>
47#include <fstream>
48#include <iostream>
49
50// sparsity_example
51//
52// usage:
53// sparsity_example [options]
54//
55// output:
56// prints the sparsity of the sparse 3 tensor specified by the basis,
57// dimension, order, given by summing over the third index, to a matrix
58// market file. This sparsity pattern yields the sparsity of the block
59// stochastic Galerkin matrix which can be visualized, e.g., by matlab.
60// The full/linear flag determines whether the third index ranges over
61// the full polynomial dimension, or only over the zeroth and first order
62// terms.
63
64// Growth policies
65const int num_growth_types = 2;
68const char *growth_type_names[] = { "slow", "moderate" };
69
70// Product Basis types
75const char *prod_basis_type_names[] = {
76 "complete", "tensor", "total", "smolyak" };
77
78// Ordering types
80const int num_ordering_types = 2;
83const char *ordering_type_names[] = {
84 "total", "lexicographic" };
85
86using Teuchos::rcp;
87using Teuchos::RCP;
88using Teuchos::ParameterList;
89using Teuchos::Array;
90
91struct Coord {
92 int i, j, k;
93 int gid;
94};
95
96template <typename coord_t>
97struct Tile {
99 Array<coord_t> parts;
100};
101
105
106int main(int argc, char **argv)
107{
108 try {
109
110 // Initialize MPI
111#ifdef HAVE_MPI
112 MPI_Init(&argc,&argv);
113#endif
114
115 // Setup command line options
116 Teuchos::CommandLineProcessor CLP;
117 CLP.setDocString(
118 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
119 int d = 3;
120 CLP.setOption("dimension", &d, "Stochastic dimension");
121 int p = 5;
122 CLP.setOption("order", &p, "Polynomial order");
123 double drop = 1.0e-12;
124 CLP.setOption("drop", &drop, "Drop tolerance");
125 std::string file = "A.mm";
126 CLP.setOption("filename", &file, "Matrix Market filename");
127 bool symmetric = true;
128 CLP.setOption("symmetric", "asymmetric", &symmetric,
129 "Use basis polynomials with symmetric PDF");
131 CLP.setOption("growth", &growth_type,
133 "Growth type");
134 ProductBasisType prod_basis_type = TOTAL;
135 CLP.setOption("product_basis", &prod_basis_type,
138 "Product basis type");
139 OrderingType ordering_type = LEXICOGRAPHIC_ORDERING;
140 CLP.setOption("ordering", &ordering_type,
143 "Product basis ordering");
144 int i_tile_size = 128;
145 CLP.setOption("tile_size", &i_tile_size, "Tile size");
146 bool save_3tensor = false;
147 CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
148 "Save full 3tensor to file");
149 std::string file_3tensor = "Cijk.dat";
150 CLP.setOption("filename_3tensor", &file_3tensor,
151 "Filename to store full 3-tensor");
152
153 // Parse arguments
154 CLP.parse( argc, argv );
155
156 // Basis
157 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
158 const double alpha = 1.0;
159 const double beta = symmetric ? 1.0 : 2.0 ;
160 for (int i=0; i<d; i++) {
161 bases[i] = rcp(new Stokhos::JacobiBasis<int,double>(
162 p, alpha, beta, true, growth_type));
163 }
164 RCP<const Stokhos::ProductBasis<int,double> > basis;
167 if (prod_basis_type == COMPLETE)
168 basis =
170 bases, drop));
171 else if (prod_basis_type == TENSOR) {
172 if (ordering_type == TOTAL_ORDERING)
173 basis =
175 bases, drop));
176 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
177 basis =
179 bases, drop));
180 }
181
182 else if (prod_basis_type == TOTAL) {
183 if (ordering_type == TOTAL_ORDERING)
184 basis =
186 bases, drop));
187 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
188 basis =
190 bases, drop));
191 }
192 else if (prod_basis_type == SMOLYAK) {
193 Stokhos::TotalOrderIndexSet<int> index_set(d, p);
194 if (ordering_type == TOTAL_ORDERING)
195 basis =
197 bases, index_set, drop));
198 else if (ordering_type == LEXICOGRAPHIC_ORDERING)
199 basis =
201 bases, index_set, drop));
202 }
203
204 // Triple product tensor
205 typedef Stokhos::Sparse3Tensor<int,double> Cijk_type;
206 RCP<Cijk_type> Cijk = basis->computeTripleProductTensor();
207
208 int basis_size = basis->size();
209 std::cout << "basis size = " << basis_size
210 << " num nonzero Cijk entries = " << Cijk->num_entries()
211 << std::endl;
212
213 // Build 2-way symmetric Cijk tensor
214 RCP<Cijk_type> Cijk_sym = rcp(new Cijk_type);
215 Cijk_type::i_iterator i_begin = Cijk->i_begin();
216 Cijk_type::i_iterator i_end = Cijk->i_end();
217 for (Cijk_type::i_iterator i_it=i_begin; i_it!=i_end; ++i_it) {
218 int i = index(i_it);
219 Cijk_type::ik_iterator k_begin = Cijk->k_begin(i_it);
220 Cijk_type::ik_iterator k_end = Cijk->k_end(i_it);
221 for (Cijk_type::ik_iterator k_it = k_begin; k_it != k_end; ++k_it) {
222 int k = index(k_it);
223 Cijk_type::ikj_iterator j_begin = Cijk->j_begin(k_it);
224 Cijk_type::ikj_iterator j_end = Cijk->j_end(k_it);
225 for (Cijk_type::ikj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
226 int j = index(j_it);
227 if (k <= j) {
228 double c = value(j_it);
229 Cijk_sym->add_term(i, j, k, c);
230 }
231 }
232 }
233 }
234 Cijk_sym->fillComplete();
235
236 // First partition based on i
237 int j_tile_size = i_tile_size / 2;
238 int num_i_parts = (basis_size + i_tile_size-1) / i_tile_size;
239 int its = basis_size / num_i_parts;
240 Array<ITile> i_tiles(num_i_parts);
241 for (int i=0; i<num_i_parts; ++i) {
242 i_tiles[i].lower = i*its;
243 i_tiles[i].upper = (i+1)*its;
244 i_tiles[i].parts.resize(1);
245 i_tiles[i].parts[0].lower = basis_size;
246 i_tiles[i].parts[0].upper = 0;
247 }
248
249 // Next partition j
250 for (Cijk_type::i_iterator i_it=Cijk_sym->i_begin();
251 i_it!=Cijk_sym->i_end(); ++i_it) {
252 int i = index(i_it);
253
254 // Find which part i belongs to
255 int idx = 0;
256 while (idx < num_i_parts && i >= i_tiles[idx].lower) ++idx;
257 --idx;
258 TEUCHOS_ASSERT(idx >= 0 && idx < num_i_parts);
259
260 Cijk_type::ik_iterator k_begin = Cijk_sym->k_begin(i_it);
261 Cijk_type::ik_iterator k_end = Cijk_sym->k_end(i_it);
262 for (Cijk_type::ik_iterator k_it = k_begin; k_it != k_end; ++k_it) {
263 int j = index(k_it); // using symmetry to interchange j and k
264
265 if (j < i_tiles[idx].parts[0].lower)
266 i_tiles[idx].parts[0].lower = j;
267 if (j > i_tiles[idx].parts[0].upper)
268 i_tiles[idx].parts[0].upper = j;
269 }
270 }
271 for (int idx=0; idx<num_i_parts; ++idx) {
272 int lower = i_tiles[idx].parts[0].lower;
273 int upper = i_tiles[idx].parts[0].upper;
274 int range = upper - lower + 1;
275 int num_j_parts = (range + j_tile_size-1) / j_tile_size;
276 int jts = range / num_j_parts;
277 Array<JTile> j_tiles(num_j_parts);
278 for (int j=0; j<num_j_parts; ++j) {
279 j_tiles[j].lower = lower + j*jts;
280 j_tiles[j].upper = lower + (j+1)*jts;
281 j_tiles[j].parts.resize(1);
282 j_tiles[j].parts[0].lower = basis_size;
283 j_tiles[j].parts[0].upper = 0;
284 }
285 i_tiles[idx].parts.swap(j_tiles);
286 }
287
288 // Now partition k
289 for (Cijk_type::i_iterator i_it=Cijk_sym->i_begin();
290 i_it!=Cijk_sym->i_end(); ++i_it) {
291 int i = index(i_it);
292
293 // Find which part i belongs to
294 int idx = 0;
295 while (idx < num_i_parts && i >= i_tiles[idx].lower) ++idx;
296 --idx;
297 TEUCHOS_ASSERT(idx >= 0 && idx < num_i_parts);
298
299 Cijk_type::ik_iterator k_begin = Cijk_sym->k_begin(i_it);
300 Cijk_type::ik_iterator k_end = Cijk_sym->k_end(i_it);
301 for (Cijk_type::ik_iterator k_it = k_begin; k_it != k_end; ++k_it) {
302 int j = index(k_it); // using symmetry to interchange j and k
303
304 // Find which part j belongs to
305 int num_j_parts = i_tiles[idx].parts.size();
306 int jdx = 0;
307 while (jdx < num_j_parts && j >= i_tiles[idx].parts[jdx].lower) ++jdx;
308 --jdx;
309 TEUCHOS_ASSERT(jdx >= 0 && jdx < num_j_parts);
310
311 Cijk_type::ikj_iterator j_begin = Cijk_sym->j_begin(k_it);
312 Cijk_type::ikj_iterator j_end = Cijk_sym->j_end(k_it);
313 for (Cijk_type::ikj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
314 int k = index(j_it); // using symmetry to interchange j and k
315
316 if (k >= j) {
317 Coord coord;
318 coord.i = i; coord.j = j; coord.k = k;
319 i_tiles[idx].parts[jdx].parts[0].parts.push_back(coord);
320 if (k < i_tiles[idx].parts[jdx].parts[0].lower)
321 i_tiles[idx].parts[jdx].parts[0].lower = k;
322 if (k > i_tiles[idx].parts[jdx].parts[0].upper)
323 i_tiles[idx].parts[jdx].parts[0].upper = k;
324 }
325 }
326 }
327 }
328
329 // Now need to divide up k-parts based on lower/upper bounds
330 int num_parts = 0;
331 int num_coord = 0;
332 for (int idx=0; idx<num_i_parts; ++idx) {
333 int num_j_parts = i_tiles[idx].parts.size();
334 for (int jdx=0; jdx<num_j_parts; ++jdx) {
335 int lower = i_tiles[idx].parts[jdx].parts[0].lower;
336 int upper = i_tiles[idx].parts[jdx].parts[0].upper;
337 int range = upper - lower + 1;
338 int num_k_parts = (range + j_tile_size-1) / j_tile_size;
339 int kts = range / num_k_parts;
340 Array<KTile> k_tiles(num_k_parts);
341 for (int k=0; k<num_k_parts; ++k) {
342 k_tiles[k].lower = lower + k*kts;
343 k_tiles[k].upper = lower + (k+1)*kts;
344 }
345 int num_k = i_tiles[idx].parts[jdx].parts[0].parts.size();
346 for (int l=0; l<num_k; ++l) {
347 int i = i_tiles[idx].parts[jdx].parts[0].parts[l].i;
348 int j = i_tiles[idx].parts[jdx].parts[0].parts[l].j;
349 int k = i_tiles[idx].parts[jdx].parts[0].parts[l].k;
350
351 // Find which part k belongs to
352 int kdx = 0;
353 while (kdx < num_k_parts && k >= k_tiles[kdx].lower) ++kdx;
354 --kdx;
355 TEUCHOS_ASSERT(kdx >= 0 && kdx < num_k_parts);
356
357 Coord coord;
358 coord.i = i; coord.j = j; coord.k = k;
359 k_tiles[kdx].parts.push_back(coord);
360 ++num_coord;
361 if (j != k) ++num_coord;
362 }
363
364 // Eliminate parts with zero size
365 Array<KTile> k_tiles2;
366 for (int k=0; k<num_k_parts; ++k) {
367 if (k_tiles[k].parts.size() > 0)
368 k_tiles2.push_back(k_tiles[k]);
369 }
370 num_parts += k_tiles2.size();
371 i_tiles[idx].parts[jdx].parts.swap(k_tiles2);
372 }
373 }
374 TEUCHOS_ASSERT(num_coord == Cijk->num_entries());
375
376 std::cout << "num parts requested = " << num_parts << std::endl;
377
378 // Form list of part IDs
379 Teuchos::Array<int> part_IDs(num_parts);
380 for (int i=0; i<num_parts; ++i)
381 part_IDs[i] = i;
382 std::random_shuffle(part_IDs.begin(), part_IDs.end());
383
384 // Assign part IDs
385 int pp = 0;
386 for (int idx=0; idx<num_i_parts; ++idx) {
387 int num_j_parts = i_tiles[idx].parts.size();
388 for (int jdx=0; jdx<num_j_parts; ++jdx) {
389 int num_k_parts = i_tiles[idx].parts[jdx].parts.size();
390 for (int kdx=0; kdx<num_k_parts; ++kdx) {
391 int num_k = i_tiles[idx].parts[jdx].parts[kdx].parts.size();
392 for (int l=0; l<num_k; ++l) {
393 i_tiles[idx].parts[jdx].parts[kdx].parts[l].gid = part_IDs[pp];
394 }
395 ++pp;
396 }
397 }
398 }
399
400 int part = 0;
401 for (int idx=0; idx<num_i_parts; ++idx) {
402 int num_j_parts = i_tiles[idx].parts.size();
403 for (int jdx=0; jdx<num_j_parts; ++jdx) {
404 int num_k_parts = i_tiles[idx].parts[jdx].parts.size();
405 for (int kdx=0; kdx<num_k_parts; ++kdx) {
406 std::cout << part++ << " : ["
407 << i_tiles[idx].lower << ","
408 << i_tiles[idx].upper << ") x ["
409 << i_tiles[idx].parts[jdx].lower << ","
410 << i_tiles[idx].parts[jdx].upper << ") x ["
411 << i_tiles[idx].parts[jdx].parts[kdx].lower << ","
412 << i_tiles[idx].parts[jdx].parts[kdx].upper << ") : "
413 << i_tiles[idx].parts[jdx].parts[kdx].parts.size()
414 << std::endl;
415 }
416 }
417 }
418
419 // Print full 3-tensor to file
420 std::ofstream cijk_file;
421 if (save_3tensor) {
422 cijk_file.open(file_3tensor.c_str());
423 cijk_file.precision(14);
424 cijk_file.setf(std::ios::scientific);
425 cijk_file << "i, j, k, part" << std::endl;
426 for (int idx=0; idx<num_i_parts; ++idx) {
427 int num_j_parts = i_tiles[idx].parts.size();
428 for (int jdx=0; jdx<num_j_parts; ++jdx) {
429 int num_k_parts = i_tiles[idx].parts[jdx].parts.size();
430 for (int kdx=0; kdx<num_k_parts; ++kdx) {
431 int num_k = i_tiles[idx].parts[jdx].parts[kdx].parts.size();
432 for (int l=0; l<num_k; ++l) {
433 Coord c = i_tiles[idx].parts[jdx].parts[kdx].parts[l];
434 cijk_file << c.i << ", " << c.j << ", " << c.k << ", " << c.gid
435 << std::endl;
436 if (c.j != c.k)
437 cijk_file << c.i << ", " << c.k << ", " << c.j << ", " << c.gid
438 << std::endl;
439 }
440 }
441 }
442 }
443 cijk_file.close();
444 }
445
446 }
447
448 catch (std::exception& e) {
449 std::cout << e.what() << std::endl;
450 }
451
452 return 0;
453}
ProductBasisType
OrderingType
int main(int argc, char **argv)
@ LEXICOGRAPHIC_ORDERING
@ TOTAL_ORDERING
const int num_ordering_types
Tile< JTile > ITile
const char * ordering_type_names[]
const OrderingType ordering_type_values[]
Tile< Coord > KTile
Tile< KTile > JTile
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
@ SMOLYAK
@ TOTAL
@ TENSOR
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,...
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
Array< coord_t > parts