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