170 int rc = Zoltan_Initialize(argc,
argv,&version);
171 TEUCHOS_ASSERT(rc == 0);
174 Teuchos::CommandLineProcessor
CLP;
176 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
178 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
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,
190 CLP.setOption(
"product_basis", &prod_basis_type,
193 "Product basis type");
195 CLP.setOption(
"ordering", &ordering_type,
198 "Product basis ordering");
199 double imbalance_tol = 1.2;
200 CLP.setOption(
"imbalance", &imbalance_tol,
"Imbalance tolerance");
202 CLP.setOption(
"full",
"linear", &
full,
"Use full or linear expansion");
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");
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++) {
221 p, alpha, beta,
true, growth_type));
223 RCP<const Stokhos::ProductBasis<int,double> > basis;
230 else if (prod_basis_type ==
TENSOR) {
240 else if (prod_basis_type ==
TOTAL) {
250 else if (prod_basis_type ==
SMOLYAK) {
255 bases, index_set, drop));
259 bases, index_set, drop));
266 Cijk = basis->computeTripleProductTensor();
268 Cijk = basis->computeLinearTripleProductTensor();
270 int basis_size = basis->size();
271 std::cout <<
"basis size = " << basis_size
272 <<
" num nonzero Cijk entries = " << Cijk->num_entries()
276 std::ofstream cijk_file;
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;
285 Zoltan_Struct *zz = Zoltan_Create(MPI_COMM_WORLD);
288 Zoltan_Set_Param(zz,
"DEBUG_LEVEL",
"2");
291 Zoltan_Set_Param(zz,
"LB_METHOD",
"RCB");
293 Zoltan_Set_Param(zz,
"NUM_GID_ENTRIES",
"1");
294 Zoltan_Set_Param(zz,
"NUM_LID_ENTRIES",
"1");
296 Zoltan_Set_Param(zz,
"RETURN_LISTS",
"PARTS");
297 Zoltan_Set_Param(zz,
"OBJ_WEIGHT_DIM",
"0");
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());
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");
318 int changes, numGidEntries, numLidEntries, numImport, numExport;
319 ZOLTAN_ID_PTR importGlobalGids, importLocalGids, exportGlobalGids, exportLocalGids;
320 int *importProcs, *importToPart, *exportProcs, *exportToPart;
337 TEUCHOS_ASSERT(rc == 0);
339 std::cout <<
"num parts requested = " << num_parts
340 <<
" changes= " << changes
341 <<
" num import = " << numImport
342 <<
" num export = " << numExport << std::endl;
367 Teuchos::Array< Teuchos::Array<CijkData> > part_list(num_parts);
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) {
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) {
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) {
381 int part = exportToPart[idx];
383 c.
i = i; c.
j =
j; c.
k = k; c.
gid = idx;
384 part_list[part].push_back(c);
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;
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;
413 std::cout <<
"max part dimension = " << dim_max
414 <<
" min part dimension = " << dim_min << std::endl;
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 << i <<
", " <<
j <<
", " << k <<
", "
431 << exportToPart[idx++] << std::endl;
439 Zoltan_LB_Free_Part(&importGlobalGids, &importLocalGids,
440 &importProcs, &importToPart);
441 Zoltan_LB_Free_Part(&exportGlobalGids, &exportLocalGids,
442 &exportProcs, &exportToPart);
448 catch (std::exception& e) {
449 std::cout << e.what() << std::endl;