84 MPI_Init(&argc,&
argv);
88 Teuchos::CommandLineProcessor
CLP;
90 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
92 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
94 CLP.setOption(
"order", &p,
"Polynomial order");
95 double drop = 1.0e-15;
96 CLP.setOption(
"drop", &drop,
"Drop tolerance");
97 std::string file_base =
"A";
98 CLP.setOption(
"base filename", &file_base,
"Base filename for matrix market files");
104 CLP.setOption(
"full",
"linear", &
full,
"Use full or linear expansion");
105 bool use_old =
false;
106 CLP.setOption(
"old",
"new", &use_old,
"Use old or new Cijk algorithm");
112 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
113 for (
int i=0; i<d; i++) {
122 Teuchos::RCP<const Stokhos::CompletePolynomialBasis<int,double> > basis =
128 Teuchos::RCP<Stokhos::Sparse3Tensor<int,double> > Cijk;
130 num_k = basis->size();
131 Cijk = basis->computeTripleProductTensor();
134 num_k = basis->dimension()+1;
135 Cijk = basis->computeLinearTripleProductTensor();
138 std::cout <<
"basis size = " << basis->size()
139 <<
" num nonzero Cijk entries = " << Cijk->num_entries()
149 int num_rows = basis->size();
158 for (Cijk_type::k_iterator k_it=Cijk->k_begin();
159 k_it!=Cijk->k_end(); ++k_it) {
162 for (Cijk_type::kj_iterator j_it = Cijk->j_begin(k_it);
163 j_it != Cijk->j_end(k_it); ++j_it) {
165 for (Cijk_type::kji_iterator i_it = Cijk->i_begin(j_it);
166 i_it != Cijk->i_end(j_it); ++i_it) {
174 std::stringstream ss;
175 ss << file_base <<
"_" << k <<
".mm";
176 std::string file = ss.str();
179 EpetraExt::RowMatrixToMatrixMarketFile(file.c_str(), mat);
183 catch (std::exception& e) {
184 std::cout << e.what() << std::endl;