43 #include "Teuchos_CommandLineProcessor.hpp" 44 #include "Teuchos_ParameterList.hpp" 76 "complete",
"tensor",
"total",
"smolyak" };
84 "total",
"lexicographic" };
88 using Teuchos::ParameterList;
96 template <
typename coord_t>
112 MPI_Init(&argc,&
argv);
116 Teuchos::CommandLineProcessor
CLP;
118 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
120 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
122 CLP.setOption(
"order", &p,
"Polynomial order");
123 double drop = 1.0e-12;
124 CLP.setOption(
"drop", &drop,
"Drop tolerance");
125 std::string file =
"A.mm";
126 CLP.setOption(
"filename", &file,
"Matrix Market filename");
127 bool symmetric =
true;
128 CLP.setOption(
"symmetric",
"asymmetric", &symmetric,
129 "Use basis polynomials with symmetric PDF");
131 CLP.setOption(
"growth", &growth_type,
135 CLP.setOption(
"product_basis", &prod_basis_type,
138 "Product basis type");
140 CLP.setOption(
"ordering", &ordering_type,
143 "Product basis ordering");
144 int i_tile_size = 128;
145 CLP.setOption(
"tile_size", &i_tile_size,
"Tile size");
146 bool save_3tensor =
false;
147 CLP.setOption(
"save_3tensor",
"no-save_3tensor", &save_3tensor,
148 "Save full 3tensor to file");
149 std::string file_3tensor =
"Cijk.dat";
150 CLP.setOption(
"filename_3tensor", &file_3tensor,
151 "Filename to store full 3-tensor");
157 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
158 const double alpha = 1.0;
159 const double beta = symmetric ? 1.0 : 2.0 ;
160 for (
int i=0; i<d; i++) {
162 p, alpha, beta,
true, growth_type));
164 RCP<const Stokhos::ProductBasis<int,double> > basis;
171 else if (prod_basis_type ==
TENSOR) {
182 else if (prod_basis_type ==
TOTAL) {
192 else if (prod_basis_type ==
SMOLYAK) {
197 bases, index_set, drop));
201 bases, index_set, drop));
206 RCP<Cijk_type> Cijk = basis->computeTripleProductTensor();
208 int basis_size = basis->size();
209 std::cout <<
"basis size = " << basis_size
210 <<
" num nonzero Cijk entries = " << Cijk->
num_entries()
214 RCP<Cijk_type> Cijk_sym = rcp(
new Cijk_type);
215 Cijk_type::i_iterator i_begin = Cijk->i_begin();
216 Cijk_type::i_iterator i_end = Cijk->i_end();
217 for (Cijk_type::i_iterator i_it=i_begin; i_it!=i_end; ++i_it) {
219 Cijk_type::ik_iterator k_begin = Cijk->k_begin(i_it);
220 Cijk_type::ik_iterator k_end = Cijk->k_end(i_it);
221 for (Cijk_type::ik_iterator k_it = k_begin; k_it != k_end; ++k_it) {
223 Cijk_type::ikj_iterator j_begin = Cijk->j_begin(k_it);
224 Cijk_type::ikj_iterator j_end = Cijk->j_end(k_it);
225 for (Cijk_type::ikj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
228 double c =
value(j_it);
229 Cijk_sym->add_term(i,
j, k, c);
234 Cijk_sym->fillComplete();
237 int j_tile_size = i_tile_size / 2;
238 int num_i_parts = (basis_size + i_tile_size-1) / i_tile_size;
239 int its = basis_size / num_i_parts;
240 Array<ITile> i_tiles(num_i_parts);
241 for (
int i=0; i<num_i_parts; ++i) {
242 i_tiles[i].lower = i*its;
243 i_tiles[i].upper = (i+1)*its;
244 i_tiles[i].parts.resize(1);
245 i_tiles[i].parts[0].lower = basis_size;
246 i_tiles[i].parts[0].upper = 0;
250 for (Cijk_type::i_iterator i_it=Cijk_sym->i_begin();
251 i_it!=Cijk_sym->i_end(); ++i_it) {
256 while (idx < num_i_parts && i >= i_tiles[idx].lower) ++idx;
258 TEUCHOS_ASSERT(idx >= 0 && idx < num_i_parts);
260 Cijk_type::ik_iterator k_begin = Cijk_sym->k_begin(i_it);
261 Cijk_type::ik_iterator k_end = Cijk_sym->k_end(i_it);
262 for (Cijk_type::ik_iterator k_it = k_begin; k_it != k_end; ++k_it) {
265 if (
j < i_tiles[idx].parts[0].lower)
266 i_tiles[idx].parts[0].lower =
j;
267 if (
j > i_tiles[idx].parts[0].upper)
268 i_tiles[idx].parts[0].upper =
j;
271 for (
int idx=0; idx<num_i_parts; ++idx) {
272 int lower = i_tiles[idx].parts[0].lower;
273 int upper = i_tiles[idx].parts[0].upper;
274 int range = upper - lower + 1;
275 int num_j_parts = (range + j_tile_size-1) / j_tile_size;
276 int jts = range / num_j_parts;
277 Array<JTile> j_tiles(num_j_parts);
278 for (
int j=0;
j<num_j_parts; ++
j) {
279 j_tiles[
j].lower = lower +
j*jts;
280 j_tiles[
j].upper = lower + (
j+1)*jts;
281 j_tiles[
j].parts.resize(1);
282 j_tiles[
j].parts[0].lower = basis_size;
283 j_tiles[
j].parts[0].upper = 0;
285 i_tiles[idx].parts.swap(j_tiles);
289 for (Cijk_type::i_iterator i_it=Cijk_sym->i_begin();
290 i_it!=Cijk_sym->i_end(); ++i_it) {
295 while (idx < num_i_parts && i >= i_tiles[idx].lower) ++idx;
297 TEUCHOS_ASSERT(idx >= 0 && idx < num_i_parts);
299 Cijk_type::ik_iterator k_begin = Cijk_sym->k_begin(i_it);
300 Cijk_type::ik_iterator k_end = Cijk_sym->k_end(i_it);
301 for (Cijk_type::ik_iterator k_it = k_begin; k_it != k_end; ++k_it) {
305 int num_j_parts = i_tiles[idx].parts.size();
307 while (jdx < num_j_parts && j >= i_tiles[idx].parts[jdx].lower) ++jdx;
309 TEUCHOS_ASSERT(jdx >= 0 && jdx < num_j_parts);
311 Cijk_type::ikj_iterator j_begin = Cijk_sym->j_begin(k_it);
312 Cijk_type::ikj_iterator j_end = Cijk_sym->j_end(k_it);
313 for (Cijk_type::ikj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
318 coord.
i = i; coord.
j =
j; coord.
k = k;
319 i_tiles[idx].parts[jdx].parts[0].parts.push_back(coord);
320 if (k < i_tiles[idx].parts[jdx].parts[0].lower)
321 i_tiles[idx].parts[jdx].parts[0].lower = k;
322 if (k > i_tiles[idx].parts[jdx].parts[0].upper)
323 i_tiles[idx].parts[jdx].parts[0].upper = k;
332 for (
int idx=0; idx<num_i_parts; ++idx) {
333 int num_j_parts = i_tiles[idx].parts.size();
334 for (
int jdx=0; jdx<num_j_parts; ++jdx) {
335 int lower = i_tiles[idx].parts[jdx].parts[0].lower;
336 int upper = i_tiles[idx].parts[jdx].parts[0].upper;
337 int range = upper - lower + 1;
338 int num_k_parts = (range + j_tile_size-1) / j_tile_size;
339 int kts = range / num_k_parts;
340 Array<KTile> k_tiles(num_k_parts);
341 for (
int k=0; k<num_k_parts; ++k) {
342 k_tiles[k].lower = lower + k*kts;
343 k_tiles[k].upper = lower + (k+1)*kts;
345 int num_k = i_tiles[idx].parts[jdx].parts[0].parts.size();
346 for (
int l=0; l<num_k; ++l) {
347 int i = i_tiles[idx].parts[jdx].parts[0].parts[l].i;
348 int j = i_tiles[idx].parts[jdx].parts[0].parts[l].j;
349 int k = i_tiles[idx].parts[jdx].parts[0].parts[l].k;
353 while (kdx < num_k_parts && k >= k_tiles[kdx].lower) ++kdx;
355 TEUCHOS_ASSERT(kdx >= 0 && kdx < num_k_parts);
358 coord.
i = i; coord.
j =
j; coord.
k = k;
359 k_tiles[kdx].parts.push_back(coord);
361 if (
j != k) ++num_coord;
365 Array<KTile> k_tiles2;
366 for (
int k=0; k<num_k_parts; ++k) {
367 if (k_tiles[k].parts.size() > 0)
368 k_tiles2.push_back(k_tiles[k]);
370 num_parts += k_tiles2.size();
371 i_tiles[idx].parts[jdx].parts.swap(k_tiles2);
374 TEUCHOS_ASSERT(num_coord == Cijk->num_entries());
376 std::cout <<
"num parts requested = " << num_parts << std::endl;
379 Teuchos::Array<int> part_IDs(num_parts);
380 for (
int i=0; i<num_parts; ++i)
382 std::random_shuffle(part_IDs.begin(), part_IDs.end());
386 for (
int idx=0; idx<num_i_parts; ++idx) {
387 int num_j_parts = i_tiles[idx].parts.size();
388 for (
int jdx=0; jdx<num_j_parts; ++jdx) {
389 int num_k_parts = i_tiles[idx].parts[jdx].parts.size();
390 for (
int kdx=0; kdx<num_k_parts; ++kdx) {
391 int num_k = i_tiles[idx].parts[jdx].parts[kdx].parts.size();
392 for (
int l=0; l<num_k; ++l) {
393 i_tiles[idx].parts[jdx].parts[kdx].parts[l].gid = part_IDs[pp];
401 for (
int idx=0; idx<num_i_parts; ++idx) {
402 int num_j_parts = i_tiles[idx].parts.size();
403 for (
int jdx=0; jdx<num_j_parts; ++jdx) {
404 int num_k_parts = i_tiles[idx].parts[jdx].parts.size();
405 for (
int kdx=0; kdx<num_k_parts; ++kdx) {
406 std::cout << part++ <<
" : [" 407 << i_tiles[idx].lower <<
"," 408 << i_tiles[idx].upper <<
") x [" 409 << i_tiles[idx].parts[jdx].lower <<
"," 410 << i_tiles[idx].parts[jdx].upper <<
") x [" 411 << i_tiles[idx].parts[jdx].parts[kdx].lower <<
"," 412 << i_tiles[idx].parts[jdx].parts[kdx].upper <<
") : " 413 << i_tiles[idx].parts[jdx].parts[kdx].parts.size()
420 std::ofstream cijk_file;
422 cijk_file.open(file_3tensor.c_str());
423 cijk_file.precision(14);
424 cijk_file.setf(std::ios::scientific);
425 cijk_file <<
"i, j, k, part" << std::endl;
426 for (
int idx=0; idx<num_i_parts; ++idx) {
427 int num_j_parts = i_tiles[idx].parts.size();
428 for (
int jdx=0; jdx<num_j_parts; ++jdx) {
429 int num_k_parts = i_tiles[idx].parts[jdx].parts.size();
430 for (
int kdx=0; kdx<num_k_parts; ++kdx) {
431 int num_k = i_tiles[idx].parts[jdx].parts[kdx].parts.size();
432 for (
int l=0; l<num_k; ++l) {
433 Coord c = i_tiles[idx].parts[jdx].parts[kdx].parts[l];
434 cijk_file << c.
i <<
", " << c.
j <<
", " << c.
k <<
", " << c.
gid 437 cijk_file << c.
i <<
", " << c.
k <<
", " << c.
j <<
", " << c.
gid 448 catch (std::exception& e) {
449 std::cout << e.what() << std::endl;
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...
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
const ProductBasisType prod_basis_type_values[]
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
int main(int argc, char **argv)
const char * ordering_type_names[]
const int num_growth_types
const char * growth_type_names[]
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.
SparseArrayIterator< index_iterator, value_iterator >::value_reference value(const SparseArrayIterator< index_iterator, value_iterator > &it)
const Stokhos::GrowthPolicy growth_type_values[]
A comparison functor implementing a strict weak ordering based lexographic ordering.
const int num_prod_basis_types
const char * prod_basis_type_names[]
const int num_ordering_types
const OrderingType ordering_type_values[]