43#include "Teuchos_CommandLineProcessor.hpp"
44#include "Teuchos_ParameterList.hpp"
45#include "Isorropia_Epetra.hpp"
46#include "Isorropia_EpetraPartitioner.hpp"
74 "hermite",
"legendre",
"clenshaw-curtis",
"gauss-patterson",
"rys",
"jacobi" };
88 "complete",
"tensor",
"total",
"smolyak" };
96 "total",
"lexicographic" };
104 "rcb",
"hg_matrix" };
108using Teuchos::ParameterList;
117 MPI_Init(&argc,&
argv);
121 Teuchos::CommandLineProcessor
CLP;
123 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
125 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
127 CLP.setOption(
"order", &p,
"Polynomial order");
128 double drop = 1.0e-12;
129 CLP.setOption(
"drop", &drop,
"Drop tolerance");
130 std::string file =
"A.mm";
131 CLP.setOption(
"filename", &file,
"Matrix Market filename");
137 CLP.setOption(
"growth", &growth_type,
141 CLP.setOption(
"product_basis", &prod_basis_type,
144 "Product basis type");
146 CLP.setOption(
"ordering", &ordering_type,
149 "Product basis ordering");
151 CLP.setOption(
"partitioning", &partitioning_type,
154 "Partitioning Method");
155 double imbalance_tol = 1.0;
156 CLP.setOption(
"imbalance", &imbalance_tol,
"Imbalance tolerance");
158 CLP.setOption(
"alpha", &alpha,
"Jacobi alpha index");
160 CLP.setOption(
"beta", &beta,
"Jacobi beta index");
162 CLP.setOption(
"full",
"linear", &
full,
"Use full or linear expansion");
164 CLP.setOption(
"tile_size", &tile_size,
"Tile size");
165 bool save_3tensor =
false;
166 CLP.setOption(
"save_3tensor",
"no-save_3tensor", &save_3tensor,
167 "Save full 3tensor to file");
168 std::string file_3tensor =
"Cijk.dat";
169 CLP.setOption(
"filename_3tensor", &file_3tensor,
170 "Filename to store full 3-tensor");
176 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
177 for (
int i=0; i<d; i++) {
180 p,
true, growth_type));
183 p,
true, growth_type));
194 p, 1.0,
true, growth_type));
197 p, alpha, beta,
true, growth_type));
199 RCP<const Stokhos::ProductBasis<int,double> > basis;
206 else if (prod_basis_type ==
TENSOR) {
217 else if (prod_basis_type ==
TOTAL) {
227 else if (prod_basis_type ==
SMOLYAK) {
232 bases, index_set, drop));
236 bases, index_set, drop));
243 Cijk = basis->computeTripleProductTensor();
245 Cijk = basis->computeLinearTripleProductTensor();
247 int basis_size = basis->size();
248 std::cout <<
"basis size = " << basis_size
249 <<
" num nonzero Cijk entries = " << Cijk->
num_entries()
259 std::ofstream cijk_file;
261 cijk_file.open(file_3tensor.c_str());
262 cijk_file.precision(14);
263 cijk_file.setf(std::ios::scientific);
264 cijk_file <<
"i, j, k, part" << std::endl;
267 Teuchos::Array<int> parts;
268 if (partitioning_type ==
RCB) {
270 int num_cijk_entries = Cijk->num_entries();
274 Cijk_type::k_iterator k_begin = Cijk->k_begin();
275 Cijk_type::k_iterator k_end = Cijk->k_end();
276 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
278 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
279 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
280 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
282 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
283 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
284 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
286 ijk_triples[0][idx] = i;
287 ijk_triples[1][idx] =
j;
288 ijk_triples[2][idx] = k;
295 ParameterList params;
296 params.set(
"partitioning method",
"rcb");
297 int num_parts = num_cijk_entries / tile_size;
298 if (num_cijk_entries % tile_size > 0)
300 std::stringstream ss;
302 params.set<std::string>(
"num parts", ss.str());
303 RCP<const Epetra_MultiVector> ijk_triples_rcp =
304 rcp(&ijk_triples,
false);
305 Isorropia::Epetra::Partitioner partitioner(ijk_triples_rcp, params);
306 partitioner.compute();
308 std::cout <<
"num parts requested = " << num_parts
309 <<
" num parts computed = " << partitioner.numProperties() << std::endl;
311 parts.resize(num_cijk_entries);
313 partitioner.extractPartsCopy(num_cijk_entries, sz, &parts[0]);
314 TEUCHOS_ASSERT(sz == num_cijk_entries);
318 for (
int i=0; i<num_cijk_entries; ++i) {
319 cijk_file << ijk_triples[0][i] <<
", "
320 << ijk_triples[1][i] <<
", "
321 << ijk_triples[2][i] <<
", "
322 << parts[i] << std::endl;
327 else if (partitioning_type ==
HG_MATRIX) {
329 RCP<const EpetraExt::MultiComm> multiComm =
337 ParameterList params;
338 params.set(
"partitioning method",
"hypergraph");
340 int num_parts = comm.
NumProc();
341 std::stringstream ss;
343 params.set<std::string>(
"num parts", ss.str());
344 std::stringstream ss2;
345 ss2 << imbalance_tol;
346 params.set<std::string>(
"imbalance tol", ss2.str());
367 Teuchos::RCP<const Epetra_CrsGraph> rebalanced_graph =
368 Teuchos::rcp(Isorropia::Epetra::createBalancedCopy(
371 std::cout << permuted_map << std::endl;
377 EpetraExt::RowMatrixToMatrixMarketFile(file.c_str(), mat);
381 Cijk_type::k_iterator k_begin = Cijk->k_begin();
382 Cijk_type::k_iterator k_end = Cijk->k_end();
383 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
385 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
386 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
387 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
389 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
390 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
391 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it){
393 cijk_file << i <<
", " <<
j <<
", " << k <<
", " << parts[i]
408 catch (std::exception& e) {
409 std::cout << e.what() << std::endl;
const PartitioningType partitioning_type_values[]
int main(int argc, char **argv)
const int num_ordering_types
const BasisType basis_type_values[]
const char * ordering_type_names[]
const OrderingType ordering_type_values[]
const int num_basis_types
const int num_partitioning_types
const int num_growth_types
const char * basis_type_names[]
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[]
int FillComplete(bool OptimizeDataStorage=true)
int PutScalar(double ScalarConstant)
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Teuchos::RCP< const Epetra_CrsGraph > getStochasticGraph() const
Get stochastic graph.
Legendre polynomial basis using Gauss-Patterson quadrature points.
Hermite polynomial basis.
Legendre 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,...
Teuchos::RCP< const EpetraExt::MultiComm > buildMultiComm(const Epetra_Comm &globalComm, int num_global_stochastic_blocks, int num_spatial_procs=-1)
GrowthPolicy
Enumerated type for determining Smolyak growth policies.