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