201 int rc = Zoltan_Initialize(argc,
argv,&version);
202 TEUCHOS_ASSERT(rc == 0);
205 Teuchos::CommandLineProcessor
CLP;
207 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
209 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
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,
221 CLP.setOption(
"product_basis", &prod_basis_type,
224 "Product basis type");
226 CLP.setOption(
"ordering", &ordering_type,
229 "Product basis ordering");
231 CLP.setOption(
"partitioning", &partitioning_type,
234 "Partitioning Method");
235 double imbalance_tol = 1.2;
236 CLP.setOption(
"imbalance", &imbalance_tol,
"Imbalance tolerance");
238 CLP.setOption(
"full",
"linear", &
full,
"Use full or linear expansion");
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");
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++) {
257 p, alpha, beta,
true, growth_type));
259 RCP<const Stokhos::ProductBasis<int,double> > basis;
266 else if (prod_basis_type ==
TENSOR) {
276 else if (prod_basis_type ==
TOTAL) {
286 else if (prod_basis_type ==
SMOLYAK) {
291 bases, index_set, drop));
295 bases, index_set, drop));
302 Cijk = basis->computeTripleProductTensor();
304 Cijk = basis->computeLinearTripleProductTensor();
306 int basis_size = basis->size();
307 std::cout <<
"basis size = " << basis_size
308 <<
" num nonzero Cijk entries = " << Cijk->num_entries()
312 std::ofstream cijk_file;
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;
321 Zoltan_Struct *zz = Zoltan_Create(MPI_COMM_WORLD);
324 Zoltan_Set_Param(zz,
"DEBUG_LEVEL",
"2");
327 Zoltan_Set_Param(zz,
"LB_METHOD",
"HYPERGRAPH");
328 Zoltan_Set_Param(zz,
"HYPERGRAPH_PACKAGE",
"PHG");
329 Zoltan_Set_Param(zz,
"NUM_GID_ENTRIES",
"1");
330 Zoltan_Set_Param(zz,
"NUM_LID_ENTRIES",
"1");
332 Zoltan_Set_Param(zz,
"RETURN_LISTS",
"PARTS");
333 Zoltan_Set_Param(zz,
"OBJ_WEIGHT_DIM",
"0");
334 Zoltan_Set_Param(zz,
"EDGE_WEIGHT_DIM",
"0");
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());
348 int changes, numGidEntries, numLidEntries, numImport, numExport;
349 ZOLTAN_ID_PTR importGlobalGids, importLocalGids, exportGlobalGids, exportLocalGids;
350 int *importProcs, *importToPart, *exportProcs, *exportToPart;
367 TEUCHOS_ASSERT(rc == 0);
369 std::cout <<
"num parts requested = " << num_parts
370 <<
" changes= " << changes
371 <<
" num import = " << numImport
372 <<
" num export = " << numExport << std::endl;
375 Array< Array<int> > part_map(num_parts);
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) {
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) {
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) {
390 if (i ==
j &&
j == k) {
391 part_map[ exportToPart[idx] ].push_back(i);
399 std::cout <<
"basis_size = " << basis_size <<
" num_diag = " << num_diag
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]);
409 TEUCHOS_ASSERT(perm_new_to_old.size() == basis_size);
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;
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) {
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) {
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) {
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;
445 Zoltan_LB_Free_Part(&importGlobalGids, &importLocalGids,
446 &importProcs, &importToPart);
447 Zoltan_LB_Free_Part(&exportGlobalGids, &exportLocalGids,
448 &exportProcs, &exportToPart);
454 catch (std::exception& e) {
455 std::cout << e.what() << std::endl;