Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
cijk_partition_zoltan.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_Epetra.hpp"
43 #include "Teuchos_CommandLineProcessor.hpp"
44 #include "Teuchos_ParameterList.hpp"
45 #include "Teuchos_toString.hpp"
46 
47 extern "C" {
48 #include "zoltan.h"
49 }
50 
51 // Growth policies
52 const int num_growth_types = 2;
55 const char *growth_type_names[] = { "slow", "moderate" };
56 
57 // Product Basis types
59 const int num_prod_basis_types = 4;
62 const char *prod_basis_type_names[] = {
63  "complete", "tensor", "total", "smolyak" };
64 
65 // Ordering types
67 const int num_ordering_types = 2;
70 const char *ordering_type_names[] = {
71  "total", "lexicographic" };
72 
73 // Partitioning types
75 const int num_partitioning_types = 2;
77  RCB, HG_FLAT_J };
78 const char *partitioning_type_names[] = {
79  "rcb", "hg_flat_j" };
80 
81 using Teuchos::rcp;
82 using Teuchos::RCP;
83 using Teuchos::ParameterList;
84 using Teuchos::Array;
85 using Teuchos::toString;
86 
87 struct TensorData {
89  RCP<const Stokhos::ProductBasis<int,double> > basis;
90  RCP<const Stokhos::Sparse3Tensor<int,double> > Cijk;
91 };
92 
93 // Functions implementing hypergraph for 1-D i-wise decomposition
94 // with flattened j. For this hypergraph model
95 // * the n vertices are the i-indices (n = basis size)
96 // * the n_k hyperedges are the flattened j-k planes:
97 // hyperedge k contains vertex i if C_ijk \neq 0 for any j
98 namespace HG_1D_Flat_J {
99 
100  // Return number of vertices
101  int get_number_of_vertices(void *data, int *ierr) {
102  TensorData *td = static_cast<TensorData*>(data);
103  *ierr = ZOLTAN_OK;
104 
105  return td->basis->size();
106  }
107 
108  // Compute IDs and weights of each vertex
109  void get_vertex_list(void *data, int sizeGID, int sizeLID,
110  ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID,
111  int wgt_dim, float *obj_wgts, int *ierr) {
112  TensorData *td = static_cast<TensorData*>(data);
113  *ierr = ZOLTAN_OK;
114 
115  int n = td->basis->size();
116  for (int i=0; i<n; ++i) {
117  globalID[i] = i;
118  localID[i] = i;
119  }
120 
121  // Do not set weights so Zoltan assumes equally weighted vertices
122  }
123 
124  // Compute number of hyperedges and pins
125  void get_hypergraph_size(void *data, int *num_lists, int *num_nonzeroes,
126  int *format, int *ierr) {
127  TensorData *td = static_cast<TensorData*>(data);
128  *ierr = ZOLTAN_OK;
129 
130  // Number of hyperedges
131  *num_lists = td->Cijk->num_k();
132 
133  // Number of pins == number of i's for all k's computing using
134  // the i-j symmetry
135  int num_pins = 0;
136  TensorData::Cijk_type::k_iterator k_begin = td->Cijk->k_begin();
137  TensorData::Cijk_type::k_iterator k_end = td->Cijk->k_end();
138  for (TensorData::Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it)
139  num_pins += td->Cijk->num_j(k_it);
140  *num_nonzeroes = num_pins;
141 
142  // hypergraph will be stored in compressed-edge format
143  *format = ZOLTAN_COMPRESSED_EDGE;
144  }
145 
146  // Compute hypergraph
147  void get_hypergraph(void *data, int sizeGID, int num_edges, int num_nonzeroes,
148  int format, ZOLTAN_ID_PTR edgeGID, int *vtxPtr,
149  ZOLTAN_ID_PTR vtxGID, int *ierr) {
150  TensorData *td = static_cast<TensorData*>(data);
151  *ierr = ZOLTAN_OK;
152 
153  // Compute pins in each hyperedge. For each hyperedge (k), these are
154  // all of the vertices (i) such that Cijk \neq 0 for any j. Due to i-j
155  // symmetry, this is all of the j's for each k such that Cijk \neq 0 for
156  // any i.
157  int kdx = 0, jdx = 0;
158  int num_pins = 0;
159  TensorData::Cijk_type::k_iterator k_begin = td->Cijk->k_begin();
160  TensorData::Cijk_type::k_iterator k_end = td->Cijk->k_end();
161  for (TensorData::Cijk_type::k_iterator k_it=k_begin; k_it!=k_end;
162  ++k_it, ++kdx) {
163  int k = index(k_it);
164  edgeGID[kdx] = k;
165  vtxPtr[kdx] = num_pins;
166  num_pins += td->Cijk->num_j(k_it);
167  TensorData::Cijk_type::kj_iterator j_begin = td->Cijk->j_begin(k_it);
168  TensorData::Cijk_type::kj_iterator j_end = td->Cijk->j_end(k_it);
169  for (TensorData::Cijk_type::kj_iterator j_it = j_begin; j_it != j_end;
170  ++j_it) {
171  int j = index(j_it);
172  vtxGID[jdx++] = j;
173  }
174  }
175  }
176 }
177 
178 
179 int main(int argc, char **argv)
180 {
181  try {
182 
183  // Initialize Zoltan
184  float version;
185  int rc = Zoltan_Initialize(argc,argv,&version);
186  TEUCHOS_ASSERT(rc == 0);
187 
188  // Setup command line options
189  Teuchos::CommandLineProcessor CLP;
190  CLP.setDocString(
191  "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
192  int d = 5;
193  CLP.setOption("dimension", &d, "Stochastic dimension");
194  int p = 3;
195  CLP.setOption("order", &p, "Polynomial order");
196  double drop = 1.0e-12;
197  CLP.setOption("drop", &drop, "Drop tolerance");
198  bool symmetric = true;
199  CLP.setOption("symmetric", "asymmetric", &symmetric, "Use basis polynomials with symmetric PDF");
201  CLP.setOption("growth", &growth_type,
203  "Growth type");
204  ProductBasisType prod_basis_type = TOTAL;
205  CLP.setOption("product_basis", &prod_basis_type,
208  "Product basis type");
209  OrderingType ordering_type = LEXICOGRAPHIC_ORDERING;
210  CLP.setOption("ordering", &ordering_type,
213  "Product basis ordering");
214  PartitioningType partitioning_type = RCB;
215  CLP.setOption("partitioning", &partitioning_type,
218  "Partitioning Method");
219  double imbalance_tol = 1.2;
220  CLP.setOption("imbalance", &imbalance_tol, "Imbalance tolerance");
221  bool full = true;
222  CLP.setOption("full", "linear", &full, "Use full or linear expansion");
223  int tile_size = 32;
224  CLP.setOption("tile_size", &tile_size, "Tile size");
225  bool save_3tensor = false;
226  CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
227  "Save full 3tensor to file");
228  std::string file_3tensor = "Cijk.dat";
229  CLP.setOption("filename_3tensor", &file_3tensor,
230  "Filename to store full 3-tensor");
231 
232  // Parse arguments
233  CLP.parse( argc, argv );
234 
235  // Basis
236  Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
237  const double alpha = 1.0;
238  const double beta = symmetric ? 1.0 : 2.0 ;
239  for (int i=0; i<d; i++) {
240  bases[i] = rcp(new Stokhos::JacobiBasis<int,double>(
241  p, alpha, beta, true, growth_type));
242  }
243  RCP<const Stokhos::ProductBasis<int,double> > basis;
246  if (prod_basis_type == COMPLETE)
247  basis =
249  bases, drop));
250  else if (prod_basis_type == TENSOR) {
251  if (ordering_type == TOTAL_ORDERING)
252  basis =
254  bases, drop));
255  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
256  basis =
258  bases, drop));
259  }
260  else if (prod_basis_type == TOTAL) {
261  if (ordering_type == TOTAL_ORDERING)
262  basis =
264  bases, drop));
265  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
266  basis =
268  bases, drop));
269  }
270  else if (prod_basis_type == SMOLYAK) {
271  Stokhos::TotalOrderIndexSet<int> index_set(d, p);
272  if (ordering_type == TOTAL_ORDERING)
273  basis =
275  bases, index_set, drop));
276  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
277  basis =
279  bases, index_set, drop));
280  }
281 
282  // Triple product tensor
284  RCP<Cijk_type> Cijk;
285  if (full)
286  Cijk = basis->computeTripleProductTensor();
287  else
288  Cijk = basis->computeLinearTripleProductTensor();
289 
290  int basis_size = basis->size();
291  std::cout << "basis size = " << basis_size
292  << " num nonzero Cijk entries = " << Cijk->num_entries()
293  << std::endl;
294 
295  // File for saving sparse Cijk tensor and parts
296  std::ofstream cijk_file;
297  if (save_3tensor) {
298  cijk_file.open(file_3tensor.c_str());
299  cijk_file.precision(14);
300  cijk_file.setf(std::ios::scientific);
301  cijk_file << "i, j, k, part" << std::endl;
302  }
303 
304  // Create zoltan
305  Zoltan_Struct *zz = Zoltan_Create(MPI_COMM_WORLD);
306 
307  // Setup Zoltan parameters
308  Zoltan_Set_Param(zz, "DEBUG_LEVEL", "2");
309 
310  // partitioning method
311  Zoltan_Set_Param(zz, "LB_METHOD", "HYPERGRAPH");
312  Zoltan_Set_Param(zz, "HYPERGRAPH_PACKAGE", "PHG"); // version of method
313  Zoltan_Set_Param(zz, "NUM_GID_ENTRIES", "1");// global IDs are integers
314  Zoltan_Set_Param(zz, "NUM_LID_ENTRIES", "1");// local IDs are integers
315  //Zoltan_Set_Param(zz, "RETURN_LISTS", "ALL"); // export AND import lists
316  Zoltan_Set_Param(zz, "RETURN_LISTS", "PARTS");
317  Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", "0"); // use Zoltan default vertex weights
318  Zoltan_Set_Param(zz, "EDGE_WEIGHT_DIM", "0");// use Zoltan default hyperedge weights
319  int num_parts = basis_size / tile_size;
320  Zoltan_Set_Param(zz, "NUM_GLOBAL_PARTS", toString(num_parts).c_str());
321  Zoltan_Set_Param(zz, "NUM_LOCAL_PARTS", toString(num_parts).c_str());
322  Zoltan_Set_Param(zz, "IMBALANCE_TOL", toString(imbalance_tol).c_str());
323 
324  // Set query functions
325  TensorData td; td.basis = basis; td.Cijk = Cijk;
326  Zoltan_Set_Num_Obj_Fn(zz, HG_1D_Flat_J::get_number_of_vertices, &td);
327  Zoltan_Set_Obj_List_Fn(zz, HG_1D_Flat_J::get_vertex_list, &td);
328  Zoltan_Set_HG_Size_CS_Fn(zz, HG_1D_Flat_J::get_hypergraph_size, &td);
329  Zoltan_Set_HG_CS_Fn(zz, HG_1D_Flat_J::get_hypergraph, &td);
330 
331  // Partition
332  int changes, numGidEntries, numLidEntries, numImport, numExport;
333  ZOLTAN_ID_PTR importGlobalGids, importLocalGids, exportGlobalGids, exportLocalGids;
334  int *importProcs, *importToPart, *exportProcs, *exportToPart;
335  rc =
336  Zoltan_LB_Partition(
337  zz, // input (all remaining fields are output)
338  &changes, // 1 if partitioning was changed, 0 otherwise
339  &numGidEntries, // Number of integers used for a global ID
340  &numLidEntries, // Number of integers used for a local ID
341  &numImport, // Number of vertices to be sent to me
342  &importGlobalGids, // Global IDs of vertices to be sent to me
343  &importLocalGids, // Local IDs of vertices to be sent to me
344  &importProcs, // Process rank for source of each incoming vertex
345  &importToPart, // New partition for each incoming vertex
346  &numExport, // Number of vertices I must send to other processes*/
347  &exportGlobalGids, // Global IDs of the vertices I must send
348  &exportLocalGids, // Local IDs of the vertices I must send
349  &exportProcs, // Process to which I send each of the vertices
350  &exportToPart); // Partition to which each vertex will belong
351  TEUCHOS_ASSERT(rc == 0);
352 
353  std::cout << "num parts requested = " << num_parts
354  << " changes= " << changes
355  << " num import = " << numImport
356  << " num export = " << numExport << std::endl;
357 
358  // for (int i=0; i<numExport; ++i)
359  // std::cout << exportGlobalGids[i] << " " << exportToPart[i] << std::endl;
360 
361  // Build list of rows that belong to each part
362  Array< Array<int> > part_map(num_parts);
363  for (int i=0; i<numExport; ++i) {
364  part_map[ exportToPart[i] ].push_back( exportGlobalGids[i] );
365  }
366 
367  // Build permuation array mapping reoredered to original
368  Array<int> perm_new_to_old;
369  for (int part=0; part<num_parts; ++part) {
370  int num_vtx = part_map[part].size();
371  for (int i=0; i<num_vtx; ++i)
372  perm_new_to_old.push_back(part_map[part][i]);
373  }
374  TEUCHOS_ASSERT(perm_new_to_old.size() == basis_size);
375 
376  // Build permuation array mapping original to reordered
377  Array<int> perm_old_to_new(basis_size);
378  for (int i=0; i<basis_size; ++i)
379  perm_old_to_new[ perm_new_to_old[i] ] = i;
380 
381  if (save_3tensor) {
382  Cijk_type::k_iterator k_begin = Cijk->k_begin();
383  Cijk_type::k_iterator k_end = Cijk->k_end();
384  for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
385  int k = index(k_it);
386  Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
387  Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
388  for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
389  int j = index(j_it);
390  Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
391  Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
392  for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
393  int i = index(i_it);
394  cijk_file << perm_old_to_new[i] << ", "
395  << perm_old_to_new[j] << ", "
396  << perm_old_to_new[k] << ", "
397  << exportToPart[i] << std::endl;
398  }
399  }
400  }
401  cijk_file.close();
402  }
403 
404  // Clean-up
405  Zoltan_LB_Free_Part(&importGlobalGids, &importLocalGids,
406  &importProcs, &importToPart);
407  Zoltan_LB_Free_Part(&exportGlobalGids, &exportLocalGids,
408  &exportProcs, &exportToPart);
409  Zoltan_Destroy(&zz);
410 
411  //Teuchos::TimeMonitor::summarize(std::cout);
412 
413  }
414  catch (std::exception& e) {
415  std::cout << e.what() << std::endl;
416  }
417 
418  return 0;
419 }
int main(int argc, char **argv)
const ProductBasisType prod_basis_type_values[]
PartitioningType
const int num_partitioning_types
const char * ordering_type_names[]
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...
void get_vertex_list(void *data, int sizeGID, int sizeLID, ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID, int wgt_dim, float *obj_wgts, int *ierr)
const int num_growth_types
const Stokhos::GrowthPolicy growth_type_values[]
const char * prod_basis_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.
Bi-directional iterator for traversing a sparse array.
RCP< const Stokhos::ProductBasis< int, double > > basis
OrderingType
const char * partitioning_type_names[]
ProductBasisType
const OrderingType ordering_type_values[]
const char * growth_type_names[]
Jacobi polynomial basis.
int get_number_of_vertices(void *data, int *ierr)
ordinal_type num_entries() const
Return number of non-zero entries.
expr expr expr expr j
const PartitioningType partitioning_type_values[]
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.
A comparison functor implementing a strict weak ordering based lexographic ordering.
const int num_ordering_types
Stokhos::Sparse3Tensor< int, double > Cijk_type
void get_hypergraph_size(void *data, int *num_lists, int *num_nonzeroes, int *format, int *ierr)
RCP< const Stokhos::Sparse3Tensor< int, double > > Cijk
void get_hypergraph(void *data, int sizeGID, int num_edges, int num_nonzeroes, int format, ZOLTAN_ID_PTR edgeGID, int *vtxPtr, ZOLTAN_ID_PTR vtxGID, int *ierr)
const int num_prod_basis_types