43 #include "Teuchos_CommandLineProcessor.hpp" 44 #include "Teuchos_ParameterList.hpp" 45 #include "Teuchos_toString.hpp" 63 "complete",
"tensor",
"total",
"smolyak" };
71 "total",
"lexicographic" };
75 using Teuchos::ParameterList;
77 using Teuchos::toString;
81 RCP<const Stokhos::ProductBasis<int,double> >
basis;
82 RCP<const Stokhos::Sparse3Tensor<int,double> >
Cijk;
98 return td->
Cijk->num_entries();
103 ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID,
104 int wgt_dim,
float *obj_wgts,
int *ierr) {
108 int n = td->
Cijk->num_entries();
109 for (
int i=0; i<n; ++i) {
127 ZOLTAN_ID_PTR globalID, ZOLTAN_ID_PTR localID,
128 int num_dim,
double *geom_vec,
int *ierr)
149 geom_vec[3*idx ] = i;
150 geom_vec[3*idx+1] =
j;
151 geom_vec[3*idx+2] = k;
167 int rc = Zoltan_Initialize(argc,
argv,&version);
168 TEUCHOS_ASSERT(rc == 0);
171 Teuchos::CommandLineProcessor
CLP;
173 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
175 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
177 CLP.setOption(
"order", &p,
"Polynomial order");
178 double drop = 1.0e-12;
179 CLP.setOption(
"drop", &drop,
"Drop tolerance");
180 bool symmetric =
true;
181 CLP.setOption(
"symmetric",
"asymmetric", &symmetric,
"Use basis polynomials with symmetric PDF");
183 CLP.setOption(
"growth", &growth_type,
187 CLP.setOption(
"product_basis", &prod_basis_type,
190 "Product basis type");
192 CLP.setOption(
"ordering", &ordering_type,
195 "Product basis ordering");
196 double imbalance_tol = 1.2;
197 CLP.setOption(
"imbalance", &imbalance_tol,
"Imbalance tolerance");
199 CLP.setOption(
"full",
"linear", &
full,
"Use full or linear expansion");
201 CLP.setOption(
"tile_size", &tile_size,
"Tile size");
202 bool save_3tensor =
false;
203 CLP.setOption(
"save_3tensor",
"no-save_3tensor", &save_3tensor,
204 "Save full 3tensor to file");
205 std::string file_3tensor =
"Cijk.dat";
206 CLP.setOption(
"filename_3tensor", &file_3tensor,
207 "Filename to store full 3-tensor");
213 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
214 const double alpha = 1.0;
215 const double beta = symmetric ? 1.0 : 2.0 ;
216 for (
int i=0; i<d; i++) {
218 p, alpha, beta,
true, growth_type));
220 RCP<const Stokhos::ProductBasis<int,double> > basis;
227 else if (prod_basis_type ==
TENSOR) {
237 else if (prod_basis_type ==
TOTAL) {
247 else if (prod_basis_type ==
SMOLYAK) {
252 bases, index_set, drop));
256 bases, index_set, drop));
263 Cijk = basis->computeTripleProductTensor();
265 Cijk = basis->computeLinearTripleProductTensor();
267 int basis_size = basis->size();
268 std::cout <<
"basis size = " << basis_size
269 <<
" num nonzero Cijk entries = " << Cijk->
num_entries()
273 std::ofstream cijk_file;
275 cijk_file.open(file_3tensor.c_str());
276 cijk_file.precision(14);
277 cijk_file.setf(std::ios::scientific);
278 cijk_file <<
"i, j, k, part" << std::endl;
282 Zoltan_Struct *zz = Zoltan_Create(MPI_COMM_WORLD);
285 Zoltan_Set_Param(zz,
"DEBUG_LEVEL",
"2");
288 Zoltan_Set_Param(zz,
"LB_METHOD",
"RCB");
290 Zoltan_Set_Param(zz,
"NUM_GID_ENTRIES",
"1");
291 Zoltan_Set_Param(zz,
"NUM_LID_ENTRIES",
"1");
293 Zoltan_Set_Param(zz,
"RETURN_LISTS",
"PARTS");
294 Zoltan_Set_Param(zz,
"OBJ_WEIGHT_DIM",
"0");
295 int num_parts = basis_size / tile_size;
296 if (num_parts * tile_size < basis_size) ++num_parts;
297 Zoltan_Set_Param(zz,
"NUM_GLOBAL_PARTS", toString(num_parts).c_str());
298 Zoltan_Set_Param(zz,
"NUM_LOCAL_PARTS", toString(num_parts).c_str());
299 Zoltan_Set_Param(zz,
"IMBALANCE_TOL", toString(imbalance_tol).c_str());
301 Zoltan_Set_Param(zz,
"KEEP_CUTS",
"1");
302 Zoltan_Set_Param(zz,
"RCB_RECTILINEAR_BLOCKS",
"1");
303 Zoltan_Set_Param(zz,
"RCB_RECOMPUTE_BOX",
"1");
304 Zoltan_Set_Param(zz,
"RCB_MAX_ASPECT_RATIO",
"3");
315 int changes, numGidEntries, numLidEntries, numImport, numExport;
316 ZOLTAN_ID_PTR importGlobalGids, importLocalGids, exportGlobalGids, exportLocalGids;
317 int *importProcs, *importToPart, *exportProcs, *exportToPart;
334 TEUCHOS_ASSERT(rc == 0);
336 std::cout <<
"num parts requested = " << num_parts
337 <<
" changes= " << changes
338 <<
" num import = " << numImport
339 <<
" num export = " << numExport << std::endl;
364 Teuchos::Array< Teuchos::Array<CijkData> > part_list(num_parts);
366 Cijk_type::k_iterator k_begin = Cijk->k_begin();
367 Cijk_type::k_iterator k_end = Cijk->k_end();
368 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
370 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
371 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
372 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
374 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
375 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
376 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
378 int part = exportToPart[idx];
380 c.
i = i; c.
j =
j; c.
k = k; c.
gid = idx;
381 part_list[part].push_back(c);
388 int dim_max = -1, dim_min = basis_size;
389 for (
int part=0; part<num_parts; ++part) {
390 int imin = basis_size, jmin = basis_size, kmin = basis_size;
391 int imax = -1, jmax = -1, kmax = -1;
392 for (
int idx=0; idx<part_list[part].size(); ++idx) {
393 if (part_list[part][idx].i < imin) imin = part_list[part][idx].i;
394 if (part_list[part][idx].
j < jmin) jmin = part_list[part][idx].j;
395 if (part_list[part][idx].k < kmin) kmin = part_list[part][idx].k;
396 if (part_list[part][idx].i > imax) imax = part_list[part][idx].i;
397 if (part_list[part][idx].
j > jmax) jmax = part_list[part][idx].j;
398 if (part_list[part][idx].k > kmax) kmax = part_list[part][idx].k;
400 int delta_i = imax - imin;
401 int delta_j = jmax - jmin;
402 int delta_k = kmax - kmin;
403 if (delta_i < dim_min) dim_min = delta_i;
404 if (delta_j < dim_min) dim_min = delta_j;
405 if (delta_k < dim_min) dim_min = delta_k;
406 if (delta_i > dim_max) dim_max = delta_i;
407 if (delta_j > dim_max) dim_max = delta_j;
408 if (delta_k > dim_max) dim_max = delta_k;
410 std::cout <<
"max part dimension = " << dim_max
411 <<
" min part dimension = " << dim_min << std::endl;
415 Cijk_type::k_iterator k_begin = Cijk->k_begin();
416 Cijk_type::k_iterator k_end = Cijk->k_end();
417 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
419 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
420 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
421 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
423 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
424 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
425 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
427 cijk_file << i <<
", " <<
j <<
", " << k <<
", " 428 << exportToPart[idx++] << std::endl;
436 Zoltan_LB_Free_Part(&importGlobalGids, &importLocalGids,
437 &importProcs, &importToPart);
438 Zoltan_LB_Free_Part(&exportGlobalGids, &exportLocalGids,
439 &exportProcs, &exportToPart);
445 catch (std::exception& e) {
446 std::cout << e.what() << std::endl;
const Stokhos::GrowthPolicy growth_type_values[]
SparseArrayIterator< index_iterator, value_iterator >::value_type index(const SparseArrayIterator< index_iterator, value_iterator > &it)
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
const char * prod_basis_type_names[]
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)
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 main(int argc, char **argv)
int get_number_of_objects(void *data, int *ierr)
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
const OrderingType ordering_type_values[]
Bi-directional iterator for traversing a sparse array.
RCP< const Stokhos::ProductBasis< int, double > > basis
ordinal_type num_entries() const
Return number of non-zero entries.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Multivariate orthogonal polynomial basis generated from a Smolyak sparse grid.
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials...
Stokhos::Sparse3Tensor< int, double > Cijk_type
An isotropic total order index set.
const int num_ordering_types
const int num_growth_types
const int num_prod_basis_types
A comparison functor implementing a strict weak ordering based lexographic ordering.
const char * growth_type_names[]
int get_num_geometry(void *data, int *ierr)
Stokhos::Sparse3Tensor< int, double > Cijk_type
RCP< const Stokhos::Sparse3Tensor< int, double > > Cijk
const char * ordering_type_names[]
const ProductBasisType prod_basis_type_values[]