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