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