Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
cijk_partition_zoltan_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_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 using Teuchos::rcp;
74 using Teuchos::RCP;
75 using Teuchos::ParameterList;
76 using Teuchos::Array;
77 using Teuchos::toString;
78 
79 struct TensorData {
81  RCP<const Stokhos::ProductBasis<int,double> > basis;
82  RCP<const Stokhos::Sparse3Tensor<int,double> > Cijk;
83 };
84 
85 struct CijkData {
86  int gid;
87  int i, j, k;
88 };
89 
90 // Functions implementing Cijk tensor for geometric partitioning
91 namespace CoordPart {
92 
93  // Return number of vertices
94  int get_number_of_objects(void *data, int *ierr) {
95  TensorData *td = static_cast<TensorData*>(data);
96  *ierr = ZOLTAN_OK;
97 
98  return td->Cijk->num_entries();
99  }
100 
101  // Compute IDs and weights of each vertex
102  void get_object_list(void *data, int sizeGID, int sizeLID,
103  ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID,
104  int wgt_dim, float *obj_wgts, int *ierr) {
105  TensorData *td = static_cast<TensorData*>(data);
106  *ierr = ZOLTAN_OK;
107 
108  int n = td->Cijk->num_entries();
109  for (int i=0; i<n; ++i) {
110  globalID[i] = i;
111  localID[i] = i;
112  }
113 
114  // Do not set weights so Zoltan assumes equally weighted vertices
115  }
116 
117  // Geometry dimension -- here 3 for (i,j,k)
118  int get_num_geometry(void *data, int *ierr)
119  {
120  *ierr = ZOLTAN_OK;
121  return 3;
122  }
123 
124  // Get list of (i,j,k) tuples
125  void get_geometry_list(void *data, int sizeGID, int sizeLID,
126  int num_obj,
127  ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID,
128  int num_dim, double *geom_vec, int *ierr)
129  {
130  TensorData *td = static_cast<TensorData*>(data);
131  *ierr = ZOLTAN_OK;
132 
133  int idx = 0;
134  TensorData::Cijk_type::k_iterator k_begin = td->Cijk->k_begin();
135  TensorData::Cijk_type::k_iterator k_end = td->Cijk->k_end();
136  for (TensorData::Cijk_type::k_iterator k_it=k_begin; k_it!=k_end;
137  ++k_it) {
138  int k = index(k_it);
139  TensorData::Cijk_type::kj_iterator j_begin = td->Cijk->j_begin(k_it);
140  TensorData::Cijk_type::kj_iterator j_end = td->Cijk->j_end(k_it);
141  for (TensorData::Cijk_type::kj_iterator j_it = j_begin; j_it != j_end;
142  ++j_it) {
143  int j = index(j_it);
144  TensorData::Cijk_type::kji_iterator i_begin = td->Cijk->i_begin(j_it);
145  TensorData::Cijk_type::kji_iterator i_end = td->Cijk->i_end(j_it);
146  for (TensorData::Cijk_type::kji_iterator i_it = i_begin; i_it != i_end;
147  ++i_it) {
148  int i = index(i_it);
149  geom_vec[3*idx ] = i;
150  geom_vec[3*idx+1] = j;
151  geom_vec[3*idx+2] = k;
152  ++idx;
153  }
154  }
155  }
156  }
157 
158 }
159 
160 
161 int main(int argc, char **argv)
162 {
163  try {
164 
165  // Initialize Zoltan
166  float version;
167  int rc = Zoltan_Initialize(argc,argv,&version);
168  TEUCHOS_ASSERT(rc == 0);
169 
170  // Setup command line options
171  Teuchos::CommandLineProcessor CLP;
172  CLP.setDocString(
173  "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
174  int d = 5;
175  CLP.setOption("dimension", &d, "Stochastic dimension");
176  int p = 3;
177  CLP.setOption("order", &p, "Polynomial order");
178  double drop = 1.0e-12;
179  CLP.setOption("drop", &drop, "Drop tolerance");
180  bool symmetric = true;
181  CLP.setOption("symmetric", "asymmetric", &symmetric, "Use basis polynomials with symmetric PDF");
183  CLP.setOption("growth", &growth_type,
185  "Growth type");
186  ProductBasisType prod_basis_type = TOTAL;
187  CLP.setOption("product_basis", &prod_basis_type,
190  "Product basis type");
191  OrderingType ordering_type = LEXICOGRAPHIC_ORDERING;
192  CLP.setOption("ordering", &ordering_type,
195  "Product basis ordering");
196  double imbalance_tol = 1.2;
197  CLP.setOption("imbalance", &imbalance_tol, "Imbalance tolerance");
198  bool full = true;
199  CLP.setOption("full", "linear", &full, "Use full or linear expansion");
200  int tile_size = 32;
201  CLP.setOption("tile_size", &tile_size, "Tile size");
202  bool save_3tensor = false;
203  CLP.setOption("save_3tensor", "no-save_3tensor", &save_3tensor,
204  "Save full 3tensor to file");
205  std::string file_3tensor = "Cijk.dat";
206  CLP.setOption("filename_3tensor", &file_3tensor,
207  "Filename to store full 3-tensor");
208 
209  // Parse arguments
210  CLP.parse( argc, argv );
211 
212  // Basis
213  Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
214  const double alpha = 1.0;
215  const double beta = symmetric ? 1.0 : 2.0 ;
216  for (int i=0; i<d; i++) {
217  bases[i] = rcp(new Stokhos::JacobiBasis<int,double>(
218  p, alpha, beta, true, growth_type));
219  }
220  RCP<const Stokhos::ProductBasis<int,double> > basis;
223  if (prod_basis_type == COMPLETE)
224  basis =
226  bases, drop));
227  else if (prod_basis_type == TENSOR) {
228  if (ordering_type == TOTAL_ORDERING)
229  basis =
231  bases, drop));
232  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
233  basis =
235  bases, drop));
236  }
237  else if (prod_basis_type == TOTAL) {
238  if (ordering_type == TOTAL_ORDERING)
239  basis =
241  bases, drop));
242  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
243  basis =
245  bases, drop));
246  }
247  else if (prod_basis_type == SMOLYAK) {
248  Stokhos::TotalOrderIndexSet<int> index_set(d, p);
249  if (ordering_type == TOTAL_ORDERING)
250  basis =
252  bases, index_set, drop));
253  else if (ordering_type == LEXICOGRAPHIC_ORDERING)
254  basis =
256  bases, index_set, drop));
257  }
258 
259  // Triple product tensor
261  RCP<Cijk_type> Cijk;
262  if (full)
263  Cijk = basis->computeTripleProductTensor();
264  else
265  Cijk = basis->computeLinearTripleProductTensor();
266 
267  int basis_size = basis->size();
268  std::cout << "basis size = " << basis_size
269  << " num nonzero Cijk entries = " << Cijk->num_entries()
270  << std::endl;
271 
272  // File for saving sparse Cijk tensor and parts
273  std::ofstream cijk_file;
274  if (save_3tensor) {
275  cijk_file.open(file_3tensor.c_str());
276  cijk_file.precision(14);
277  cijk_file.setf(std::ios::scientific);
278  cijk_file << "i, j, k, part" << std::endl;
279  }
280 
281  // Create zoltan
282  Zoltan_Struct *zz = Zoltan_Create(MPI_COMM_WORLD);
283 
284  // Setup Zoltan parameters
285  Zoltan_Set_Param(zz, "DEBUG_LEVEL", "2");
286 
287  // partitioning method
288  Zoltan_Set_Param(zz, "LB_METHOD", "RCB");
289  //Zoltan_Set_Param(zz, "LB_METHOD", "HSFC");
290  Zoltan_Set_Param(zz, "NUM_GID_ENTRIES", "1");// global IDs are integers
291  Zoltan_Set_Param(zz, "NUM_LID_ENTRIES", "1");// local IDs are integers
292  //Zoltan_Set_Param(zz, "RETURN_LISTS", "ALL"); // export AND import lists
293  Zoltan_Set_Param(zz, "RETURN_LISTS", "PARTS");
294  Zoltan_Set_Param(zz, "OBJ_WEIGHT_DIM", "0"); // use Zoltan default vertex weights
295  int num_parts = basis_size / tile_size;
296  if (num_parts * tile_size < basis_size) ++num_parts;
297  Zoltan_Set_Param(zz, "NUM_GLOBAL_PARTS", toString(num_parts).c_str());
298  Zoltan_Set_Param(zz, "NUM_LOCAL_PARTS", toString(num_parts).c_str());
299  Zoltan_Set_Param(zz, "IMBALANCE_TOL", toString(imbalance_tol).c_str());
300 
301  Zoltan_Set_Param(zz, "KEEP_CUTS", "1");
302  Zoltan_Set_Param(zz, "RCB_RECTILINEAR_BLOCKS", "1");
303  Zoltan_Set_Param(zz, "RCB_RECOMPUTE_BOX", "1");
304  Zoltan_Set_Param(zz, "RCB_MAX_ASPECT_RATIO", "3");
305 
306  // Set query functions
307  TensorData td; td.basis = basis; td.Cijk = Cijk;
308 
309  Zoltan_Set_Num_Obj_Fn(zz, CoordPart::get_number_of_objects, &td);
310  Zoltan_Set_Obj_List_Fn(zz, CoordPart::get_object_list, &td);
311  Zoltan_Set_Num_Geom_Fn(zz, CoordPart::get_num_geometry, &td);
312  Zoltan_Set_Geom_Multi_Fn(zz, CoordPart::get_geometry_list, &td);
313 
314  // Partition
315  int changes, numGidEntries, numLidEntries, numImport, numExport;
316  ZOLTAN_ID_PTR importGlobalGids, importLocalGids, exportGlobalGids, exportLocalGids;
317  int *importProcs, *importToPart, *exportProcs, *exportToPart;
318  rc =
319  Zoltan_LB_Partition(
320  zz, // input (all remaining fields are output)
321  &changes, // 1 if partitioning was changed, 0 otherwise
322  &numGidEntries, // Number of integers used for a global ID
323  &numLidEntries, // Number of integers used for a local ID
324  &numImport, // Number of vertices to be sent to me
325  &importGlobalGids, // Global IDs of vertices to be sent to me
326  &importLocalGids, // Local IDs of vertices to be sent to me
327  &importProcs, // Process rank for source of each incoming vertex
328  &importToPart, // New partition for each incoming vertex
329  &numExport, // Number of vertices I must send to other processes*/
330  &exportGlobalGids, // Global IDs of the vertices I must send
331  &exportLocalGids, // Local IDs of the vertices I must send
332  &exportProcs, // Process to which I send each of the vertices
333  &exportToPart); // Partition to which each vertex will belong
334  TEUCHOS_ASSERT(rc == 0);
335 
336  std::cout << "num parts requested = " << num_parts
337  << " changes= " << changes
338  << " num import = " << numImport
339  << " num export = " << numExport << std::endl;
340 
341  // Compute max/min bounding box dimensions
342  // double dim_max = 0.0, dim_min = 0.0;
343  // for (int part=0; part<num_parts; ++part) {
344  // double xmin, ymin, zmin, xmax, ymax, zmax;
345  // int ndim;
346  // rc = Zoltan_RCB_Box(zz, part, &ndim, &xmin, &ymin, &zmin,
347  // &xmax, &ymax, &zmax);
348  // TEUCHOS_ASSERT(rc == 0);
349  // double delta_x = xmax - xmin;
350  // double delta_y = ymax - ymin;
351  // double delta_z = zmax - zmin;
352  // std::cout << delta_x << " " << delta_y << " " << delta_z << std::endl;
353  // if (delta_x > dim_max) dim_max = delta_x;
354  // if (delta_y > dim_max) dim_max = delta_y;
355  // if (delta_z > dim_max) dim_max = delta_z;
356  // if (delta_x < dim_min) dim_min = delta_x;
357  // if (delta_y < dim_min) dim_min = delta_y;
358  // if (delta_z < dim_min) dim_min = delta_z;
359  // }
360  // std::cout << "max part dimension = " << dim_max
361  // << " min part dimension = " << dim_min << std::endl;
362 
363  // Build part list
364  Teuchos::Array< Teuchos::Array<CijkData> > part_list(num_parts);
365  int idx = 0;
366  Cijk_type::k_iterator k_begin = Cijk->k_begin();
367  Cijk_type::k_iterator k_end = Cijk->k_end();
368  for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
369  int k = index(k_it);
370  Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
371  Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
372  for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
373  int j = index(j_it);
374  Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
375  Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
376  for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
377  int i = index(i_it);
378  int part = exportToPart[idx];
379  CijkData c;
380  c.i = i; c.j = j; c.k = k; c.gid = idx;
381  part_list[part].push_back(c);
382  ++idx;
383  }
384  }
385  }
386 
387  // Compute max/min bounding box dimensions
388  int dim_max = -1, dim_min = basis_size;
389  for (int part=0; part<num_parts; ++part) {
390  int imin = basis_size, jmin = basis_size, kmin = basis_size;
391  int imax = -1, jmax = -1, kmax = -1;
392  for (int idx=0; idx<part_list[part].size(); ++idx) {
393  if (part_list[part][idx].i < imin) imin = part_list[part][idx].i;
394  if (part_list[part][idx].j < jmin) jmin = part_list[part][idx].j;
395  if (part_list[part][idx].k < kmin) kmin = part_list[part][idx].k;
396  if (part_list[part][idx].i > imax) imax = part_list[part][idx].i;
397  if (part_list[part][idx].j > jmax) jmax = part_list[part][idx].j;
398  if (part_list[part][idx].k > kmax) kmax = part_list[part][idx].k;
399  }
400  int delta_i = imax - imin;
401  int delta_j = jmax - jmin;
402  int delta_k = kmax - kmin;
403  if (delta_i < dim_min) dim_min = delta_i;
404  if (delta_j < dim_min) dim_min = delta_j;
405  if (delta_k < dim_min) dim_min = delta_k;
406  if (delta_i > dim_max) dim_max = delta_i;
407  if (delta_j > dim_max) dim_max = delta_j;
408  if (delta_k > dim_max) dim_max = delta_k;
409  }
410  std::cout << "max part dimension = " << dim_max
411  << " min part dimension = " << dim_min << std::endl;
412 
413  if (save_3tensor) {
414  int 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 << i << ", " << j << ", " << k << ", "
428  << exportToPart[idx++] << std::endl;
429  }
430  }
431  }
432  cijk_file.close();
433  }
434 
435  // Clean-up
436  Zoltan_LB_Free_Part(&importGlobalGids, &importLocalGids,
437  &importProcs, &importToPart);
438  Zoltan_LB_Free_Part(&exportGlobalGids, &exportLocalGids,
439  &exportProcs, &exportToPart);
440  Zoltan_Destroy(&zz);
441 
442  //Teuchos::TimeMonitor::summarize(std::cout);
443 
444  }
445  catch (std::exception& e) {
446  std::cout << e.what() << std::endl;
447  }
448 
449  return 0;
450 }
const Stokhos::GrowthPolicy growth_type_values[]
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 char * prod_basis_type_names[]
void get_object_list(void *data, int sizeGID, int sizeLID, ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID, int wgt_dim, float *obj_wgts, int *ierr)
void get_geometry_list(void *data, int sizeGID, int sizeLID, int num_obj, ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID, int num_dim, double *geom_vec, int *ierr)
int main(int argc, char **argv)
int get_number_of_objects(void *data, int *ierr)
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
const OrderingType ordering_type_values[]
Bi-directional iterator for traversing a sparse array.
RCP< const Stokhos::ProductBasis< int, double > > basis
OrderingType
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.
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials...
Stokhos::Sparse3Tensor< int, double > Cijk_type
An isotropic total order index set.
const int num_ordering_types
const int num_growth_types
const int num_prod_basis_types
A comparison functor implementing a strict weak ordering based lexographic ordering.
const char * growth_type_names[]
int get_num_geometry(void *data, int *ierr)
Stokhos::Sparse3Tensor< int, double > Cijk_type
RCP< const Stokhos::Sparse3Tensor< int, double > > Cijk
const char * ordering_type_names[]
const ProductBasisType prod_basis_type_values[]