Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
cijk_partition_zoltan_3d.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 3-D decomposition
94 // For this hypergraph model
95 // * the nnz vertices are the Cijk non-zeros
96 // * the 3*n hyperedges are the i, j, and k values
97 // each Cijk non-zero belongs to 3 hyperedges: i, j, and k
98 namespace HG_3D {
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->Cijk->num_entries();
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 nnz = td->Cijk->num_entries();
116  for (int i=0; i<nnz; ++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_pins,
126  int *format, int *ierr) {
127  TensorData *td = static_cast<TensorData*>(data);
128  *ierr = ZOLTAN_OK;
129 
130  //int n = td->basis->size();
131  int nnz = td->Cijk->num_entries();
132 
133  // Number of vertices
134  *num_lists = nnz;
135 
136  // Number of pins. Each nonzero belongs creates 1 pin in 3 hyperedges
137  // thus there are 3*nnz pins
138  *num_pins = 3*nnz;
139 
140  // hypergraph will be stored in compressed-vertex format
141  *format = ZOLTAN_COMPRESSED_VERTEX;
142  }
143 
144  // Compute hypergraph
145  void get_hypergraph(void *data, int sizeGID, int num_vtx, int num_pins,
146  int format, ZOLTAN_ID_PTR vtxGID, int *edgePtr,
147  ZOLTAN_ID_PTR edgeGID, int *ierr) {
148  TensorData *td = static_cast<TensorData*>(data);
149  *ierr = ZOLTAN_OK;
150 
151  int n = td->basis->size();
152 
153  TEUCHOS_ASSERT(sizeGID == 1);
154  TEUCHOS_ASSERT(num_vtx == td->Cijk->num_entries());
155  TEUCHOS_ASSERT(num_pins == 3*(td->Cijk->num_entries()));
156 
157  // Compute pins in each hyperedge stored in compressed-vertex format.
158  // For each vertex we store the GIDs of the 3 edges that it connects to.
159  // Edges are ordered as follows:
160  // [0,n) -- i edges
161  // [n,2*n) -- j edges
162  // [2*n,3*n) -- k edges
163  int vtx_idx = 0;
164  int pin_idx = 0;
165  TensorData::Cijk_type::k_iterator k_begin = td->Cijk->k_begin();
166  TensorData::Cijk_type::k_iterator k_end = td->Cijk->k_end();
167  for (TensorData::Cijk_type::k_iterator k_it=k_begin; k_it!=k_end;
168  ++k_it) {
169  int k = index(k_it);
170  TensorData::Cijk_type::kj_iterator j_begin = td->Cijk->j_begin(k_it);
171  TensorData::Cijk_type::kj_iterator j_end = td->Cijk->j_end(k_it);
172  for (TensorData::Cijk_type::kj_iterator j_it = j_begin; j_it != j_end;
173  ++j_it) {
174  int j = index(j_it);
175  TensorData::Cijk_type::kji_iterator i_begin = td->Cijk->i_begin(j_it);
176  TensorData::Cijk_type::kji_iterator i_end = td->Cijk->i_end(j_it);
177  for (TensorData::Cijk_type::kji_iterator i_it = i_begin; i_it != i_end;
178  ++i_it) {
179  int i = index(i_it);
180  vtxGID[vtx_idx] = vtx_idx;
181  edgePtr[vtx_idx++] = pin_idx;
182  edgeGID[pin_idx++] = i;
183  edgeGID[pin_idx++] = n + j;
184  edgeGID[pin_idx++] = 2*n + k;
185  }
186  }
187  }
188  }
189 }
190 
191 
192 int main(int argc, char **argv)
193 {
194  try {
195 
196  // Initialize Zoltan
197  float version;
198  int rc = Zoltan_Initialize(argc,argv,&version);
199  TEUCHOS_ASSERT(rc == 0);
200 
201  // Setup command line options
202  Teuchos::CommandLineProcessor CLP;
203  CLP.setDocString(
204  "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
205  int d = 5;
206  CLP.setOption("dimension", &d, "Stochastic dimension");
207  int p = 3;
208  CLP.setOption("order", &p, "Polynomial order");
209  double drop = 1.0e-12;
210  CLP.setOption("drop", &drop, "Drop tolerance");
211  bool symmetric = true;
212  CLP.setOption("symmetric", "asymmetric", &symmetric, "Use basis polynomials with symmetric PDF");
214  CLP.setOption("growth", &growth_type,
216  "Growth type");
217  ProductBasisType prod_basis_type = TOTAL;
218  CLP.setOption("product_basis", &prod_basis_type,
221  "Product basis type");
222  OrderingType ordering_type = LEXICOGRAPHIC_ORDERING;
223  CLP.setOption("ordering", &ordering_type,
226  "Product basis ordering");
227  PartitioningType partitioning_type = RCB;
228  CLP.setOption("partitioning", &partitioning_type,
231  "Partitioning Method");
232  double imbalance_tol = 1.2;
233  CLP.setOption("imbalance", &imbalance_tol, "Imbalance tolerance");
234  bool full = true;
235  CLP.setOption("full", "linear", &full, "Use full or linear expansion");
236  int tile_size = 32;
237  CLP.setOption("tile_size", &tile_size, "Tile size");
238  bool save_3tensor = false;
239  CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
240  "Save full 3tensor to file");
241  std::string file_3tensor = "Cijk.dat";
242  CLP.setOption("filename_3tensor", &file_3tensor,
243  "Filename to store full 3-tensor");
244 
245  // Parse arguments
246  CLP.parse( argc, argv );
247 
248  // Basis
249  Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
250  const double alpha = 1.0;
251  const double beta = symmetric ? 1.0 : 2.0 ;
252  for (int i=0; i<d; i++) {
253  bases[i] = rcp(new Stokhos::JacobiBasis<int,double>(
254  p, alpha, beta, true, growth_type));
255  }
256  RCP<const Stokhos::ProductBasis<int,double> > basis;
259  if (prod_basis_type == COMPLETE)
260  basis =
262  bases, drop));
263  else if (prod_basis_type == TENSOR) {
264  if (ordering_type == TOTAL_ORDERING)
265  basis =
267  bases, drop));
268  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
269  basis =
271  bases, drop));
272  }
273  else if (prod_basis_type == TOTAL) {
274  if (ordering_type == TOTAL_ORDERING)
275  basis =
277  bases, drop));
278  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
279  basis =
281  bases, drop));
282  }
283  else if (prod_basis_type == SMOLYAK) {
284  Stokhos::TotalOrderIndexSet<int> index_set(d, p);
285  if (ordering_type == TOTAL_ORDERING)
286  basis =
288  bases, index_set, drop));
289  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
290  basis =
292  bases, index_set, drop));
293  }
294 
295  // Triple product tensor
297  RCP<Cijk_type> Cijk;
298  if (full)
299  Cijk = basis->computeTripleProductTensor();
300  else
301  Cijk = basis->computeLinearTripleProductTensor();
302 
303  int basis_size = basis->size();
304  std::cout << "basis size = " << basis_size
305  << " num nonzero Cijk entries = " << Cijk->num_entries()
306  << std::endl;
307 
308  // File for saving sparse Cijk tensor and parts
309  std::ofstream cijk_file;
310  if (save_3tensor) {
311  cijk_file.open(file_3tensor.c_str());
312  cijk_file.precision(14);
313  cijk_file.setf(std::ios::scientific);
314  cijk_file << "i, j, k, part" << std::endl;
315  }
316 
317  // Create zoltan
318  Zoltan_Struct *zz = Zoltan_Create(MPI_COMM_WORLD);
319 
320  // Setup Zoltan parameters
321  Zoltan_Set_Param(zz, "DEBUG_LEVEL", "2");
322 
323  // partitioning method
324  Zoltan_Set_Param(zz, "LB_METHOD", "HYPERGRAPH");
325  Zoltan_Set_Param(zz, "HYPERGRAPH_PACKAGE", "PHG"); // version of method
326  Zoltan_Set_Param(zz, "NUM_GID_ENTRIES", "1");// global IDs are integers
327  Zoltan_Set_Param(zz, "NUM_LID_ENTRIES", "1");// local IDs are integers
328  //Zoltan_Set_Param(zz, "RETURN_LISTS", "ALL"); // export AND import lists
329  Zoltan_Set_Param(zz, "RETURN_LISTS", "PARTS");
330  Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", "0"); // use Zoltan default vertex weights
331  Zoltan_Set_Param(zz, "EDGE_WEIGHT_DIM", "0");// use Zoltan default hyperedge weights
332  int num_parts = basis_size / tile_size;
333  Zoltan_Set_Param(zz, "NUM_GLOBAL_PARTS", toString(num_parts).c_str());
334  Zoltan_Set_Param(zz, "NUM_LOCAL_PARTS", toString(num_parts).c_str());
335  Zoltan_Set_Param(zz, "IMBALANCE_TOL", toString(imbalance_tol).c_str());
336 
337  // Set query functions
338  TensorData td; td.basis = basis; td.Cijk = Cijk;
339  Zoltan_Set_Num_Obj_Fn(zz, HG_3D::get_number_of_vertices, &td);
340  Zoltan_Set_Obj_List_Fn(zz, HG_3D::get_vertex_list, &td);
341  Zoltan_Set_HG_Size_CS_Fn(zz, HG_3D::get_hypergraph_size, &td);
342  Zoltan_Set_HG_CS_Fn(zz, HG_3D::get_hypergraph, &td);
343 
344  // Partition
345  int changes, numGidEntries, numLidEntries, numImport, numExport;
346  ZOLTAN_ID_PTR importGlobalGids, importLocalGids, exportGlobalGids, exportLocalGids;
347  int *importProcs, *importToPart, *exportProcs, *exportToPart;
348  rc =
349  Zoltan_LB_Partition(
350  zz, // input (all remaining fields are output)
351  &changes, // 1 if partitioning was changed, 0 otherwise
352  &numGidEntries, // Number of integers used for a global ID
353  &numLidEntries, // Number of integers used for a local ID
354  &numImport, // Number of vertices to be sent to me
355  &importGlobalGids, // Global IDs of vertices to be sent to me
356  &importLocalGids, // Local IDs of vertices to be sent to me
357  &importProcs, // Process rank for source of each incoming vertex
358  &importToPart, // New partition for each incoming vertex
359  &numExport, // Number of vertices I must send to other processes*/
360  &exportGlobalGids, // Global IDs of the vertices I must send
361  &exportLocalGids, // Local IDs of the vertices I must send
362  &exportProcs, // Process to which I send each of the vertices
363  &exportToPart); // Partition to which each vertex will belong
364  TEUCHOS_ASSERT(rc == 0);
365 
366  std::cout << "num parts requested = " << num_parts
367  << " changes= " << changes
368  << " num import = " << numImport
369  << " num export = " << numExport << std::endl;
370 
371  // Build list of rows that belong to each part based on diagonal
372  Array< Array<int> > part_map(num_parts);
373  int idx = 0;
374  int num_diag = 0;
375  Cijk_type::k_iterator k_begin = Cijk->k_begin();
376  Cijk_type::k_iterator k_end = Cijk->k_end();
377  for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
378  int k = index(k_it);
379  Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
380  Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
381  for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
382  int j = index(j_it);
383  Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
384  Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
385  for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
386  int i = index(i_it);
387  if (i == j && j == k) {
388  part_map[ exportToPart[idx] ].push_back(i);
389  ++num_diag;
390  }
391  idx++;
392  }
393  }
394  }
395 
396  std::cout << "basis_size = " << basis_size << " num_diag = " << num_diag
397  << std::endl;
398 
399  // Build permuation array mapping reoredered to original
400  Array<int> perm_new_to_old;
401  for (int part=0; part<num_parts; ++part) {
402  int num_row = part_map[part].size();
403  for (int i=0; i<num_row; ++i)
404  perm_new_to_old.push_back(part_map[part][i]);
405  }
406  TEUCHOS_ASSERT(perm_new_to_old.size() == basis_size);
407 
408  // Build permuation array mapping original to reordered
409  Array<int> perm_old_to_new(basis_size);
410  for (int i=0; i<basis_size; ++i)
411  perm_old_to_new[ perm_new_to_old[i] ] = i;
412 
413  if (save_3tensor) {
414  idx = 0;
415  Cijk_type::k_iterator k_begin = Cijk->k_begin();
416  Cijk_type::k_iterator k_end = Cijk->k_end();
417  for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
418  int k = index(k_it);
419  Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
420  Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
421  for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
422  int j = index(j_it);
423  Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
424  Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
425  for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
426  int i = index(i_it);
427  cijk_file << perm_old_to_new[i] << ", "
428  << perm_old_to_new[j] << ", "
429  << perm_old_to_new[k] << ", "
430  << exportToPart[idx++] << std::endl;
431  // cijk_file << i << ", "
432  // << j << ", "
433  // << k << ", "
434  // << exportToPart[idx++] << std::endl;
435  }
436  }
437  }
438  cijk_file.close();
439  }
440 
441  // Clean-up
442  Zoltan_LB_Free_Part(&importGlobalGids, &importLocalGids,
443  &importProcs, &importToPart);
444  Zoltan_LB_Free_Part(&exportGlobalGids, &exportLocalGids,
445  &exportProcs, &exportToPart);
446  Zoltan_Destroy(&zz);
447 
448  //Teuchos::TimeMonitor::summarize(std::cout);
449 
450  }
451  catch (std::exception& e) {
452  std::cout << e.what() << std::endl;
453  }
454 
455  return 0;
456 }
void get_hypergraph_size(void *data, int *num_lists, int *num_pins, int *format, int *ierr)
PartitioningType
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 int num_partitioning_types
int main(int argc, char **argv)
const char * partitioning_type_names[]
int get_number_of_vertices(void *data, int *ierr)
const OrderingType ordering_type_values[]
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
const int num_growth_types
const Stokhos::GrowthPolicy growth_type_values[]
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)
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
void get_hypergraph(void *data, int sizeGID, int num_vtx, int num_pins, int format, ZOLTAN_ID_PTR vtxGID, int *edgePtr, ZOLTAN_ID_PTR edgeGID, int *ierr)
ProductBasisType
Jacobi polynomial basis.
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.
const char * growth_type_names[]
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials...
const ProductBasisType prod_basis_type_values[]
Stokhos::Sparse3Tensor< int, double > Cijk_type
const int num_prod_basis_types
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
const char * ordering_type_names[]
RCP< const Stokhos::Sparse3Tensor< int, double > > Cijk
const PartitioningType partitioning_type_values[]
const char * prod_basis_type_names[]