106 MPI_Init(&argc,&
argv);
110 Teuchos::CommandLineProcessor
CLP;
112 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
114 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
116 CLP.setOption(
"order", &p,
"Polynomial order");
117 double drop = 1.0e-12;
118 CLP.setOption(
"drop", &drop,
"Drop tolerance");
119 std::string file =
"A.mm";
120 CLP.setOption(
"filename", &file,
"Matrix Market filename");
126 CLP.setOption(
"growth", &growth_type,
130 CLP.setOption(
"product_basis", &prod_basis_type,
133 "Product basis type");
135 CLP.setOption(
"ordering", &ordering_type,
138 "Product basis ordering");
140 CLP.setOption(
"alpha", &alpha,
"Jacobi alpha index");
142 CLP.setOption(
"beta", &beta,
"Jacobi beta index");
144 CLP.setOption(
"full",
"linear", &
full,
"Use full or linear expansion");
145 bool use_old =
false;
146 CLP.setOption(
"old",
"new", &use_old,
"Use old or new Cijk algorithm");
148 CLP.setOption(
"print",
"no-print", &print,
"Print Cijk to screen");
149 bool save_3tensor =
false;
150 CLP.setOption(
"save_3tensor",
"no-save_3tensor", &save_3tensor,
151 "Save full 3tensor to file");
152 std::string file_3tensor =
"Cijk.dat";
153 CLP.setOption(
"filename_3tensor", &file_3tensor,
154 "Filename to store full 3-tensor");
156 CLP.setOption(
"unique",
"no-unique", &unique,
157 "Only save the unique non-zeros");
159 CLP.setOption(
"rcm",
"no-rcm", &rcm,
"Reorder using RCM");
165 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
166 for (
int i=0; i<d; i++) {
169 p,
true, growth_type));
172 p,
true, growth_type));
183 p, 1.0,
true, growth_type));
186 p, alpha, beta,
true, growth_type));
188 Teuchos::RCP<const Stokhos::ProductBasis<int,double> > basis;
195 bases, drop, use_old));
196 else if (prod_basis_type ==
TENSOR) {
210 else if (prod_basis_type ==
TOTAL) {
224 else if (prod_basis_type ==
SMOLYAK) {
229 bases, index_set, drop));
233 bases, index_set, drop));
237 bases, index_set, drop));
242 Teuchos::RCP<Cijk_type> Cijk;
244 Cijk = basis->computeTripleProductTensor();
246 Cijk = basis->computeLinearTripleProductTensor();
248 std::cout <<
"basis size = " << basis->size()
249 <<
" num nonzero Cijk entries = " << Cijk->
num_entries()
259 Teuchos::RCP<Cijk_type> Cijk3 = Teuchos::rcp(
new Cijk_type);
261 Cijk_type::k_iterator k_begin = Cijk->k_begin();
262 Cijk_type::k_iterator k_end = Cijk->k_end();
263 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
265 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
266 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
267 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
269 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
270 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
271 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
273 double cijk = value(i_it);
274 if (i != 0 || (i == 0 &&
j == 0 && k == 0))
275 Cijk3->add_term(i,
j, k, cijk);
280 Cijk3->fillComplete();
282 Teuchos::RCP<Epetra_CrsGraph> graph =
287 Ifpack_RCMReordering ifpack_rcm;
288 ifpack_rcm.SetParameter(
"reorder: root node", basis->size()-1);
289 ifpack_rcm.Compute(mat);
291 Teuchos::RCP<Cijk_type> Cijk2 = Teuchos::rcp(
new Cijk_type);
292 Cijk_type::k_iterator k_begin = Cijk->k_begin();
293 Cijk_type::k_iterator k_end = Cijk->k_end();
294 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
295 int k = ifpack_rcm.Reorder(index(k_it));
296 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
297 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
298 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
299 int j = ifpack_rcm.Reorder(index(j_it));
300 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
301 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
302 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
303 int i = ifpack_rcm.Reorder(index(i_it));
304 double cijk = value(i_it);
305 Cijk2->add_term(i,
j, k, cijk);
309 Cijk2->fillComplete();
314 std::cout << *Cijk << std::endl;
322 std::ofstream cijk_file(file_3tensor.c_str());
323 cijk_file.precision(14);
324 cijk_file.setf(std::ios::scientific);
325 cijk_file <<
"i, j, k, cijk" << std::endl;
326 Cijk_type::k_iterator k_begin = Cijk->k_begin();
327 Cijk_type::k_iterator k_end = Cijk->k_end();
328 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
330 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
331 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
332 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
334 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
335 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
336 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
338 double cijk = value(i_it);
339 if (!unique || ( i >=
j &&
j >= k ))
340 cijk_file << i <<
", "
343 << cijk << std::endl;
350 Teuchos::TimeMonitor::summarize(std::cout);
353 catch (std::exception& e) {
354 std::cout << e.what() << std::endl;
Teuchos::RCP< Epetra_CrsGraph > sparse3Tensor2CrsGraph(const Stokhos::OrthogPolyBasis< ordinal_type, value_type > &basis, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, const Epetra_Comm &comm)
Build an Epetra_CrsGraph from a sparse 3 tensor.