43 #include "Teuchos_CommandLineProcessor.hpp" 46 #include "Epetra_MpiComm.h" 48 #include "Epetra_SerialComm.h" 51 #include "Ifpack_RCMReordering.h" 73 "hermite",
"legendre",
"clenshaw-curtis",
"gauss-patterson",
"rys",
"jacobi" };
87 "complete",
"tensor",
"total",
"smolyak" };
95 "total",
"lexicographic",
"morton-z" };
103 MPI_Init(&argc,&
argv);
107 Teuchos::CommandLineProcessor
CLP;
109 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
111 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
113 CLP.setOption(
"order", &p,
"Polynomial order");
114 double drop = 1.0e-12;
115 CLP.setOption(
"drop", &drop,
"Drop tolerance");
116 std::string file =
"A.mm";
117 CLP.setOption(
"filename", &file,
"Matrix Market filename");
123 CLP.setOption(
"growth", &growth_type,
127 CLP.setOption(
"product_basis", &prod_basis_type,
130 "Product basis type");
132 CLP.setOption(
"ordering", &ordering_type,
135 "Product basis ordering");
137 CLP.setOption(
"alpha", &alpha,
"Jacobi alpha index");
139 CLP.setOption(
"beta", &beta,
"Jacobi beta index");
141 CLP.setOption(
"full",
"linear", &
full,
"Use full or linear expansion");
142 bool use_old =
false;
143 CLP.setOption(
"old",
"new", &use_old,
"Use old or new Cijk algorithm");
145 CLP.setOption(
"print",
"no-print", &print,
"Print Cijk to screen");
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");
153 CLP.setOption(
"unique",
"no-unique", &unique,
154 "Only save the unique non-zeros");
156 CLP.setOption(
"rcm",
"no-rcm", &rcm,
"Reorder using RCM");
162 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
163 for (
int i=0; i<d; i++) {
166 p,
true, growth_type));
169 p,
true, growth_type));
180 p, 1.0,
true, growth_type));
183 p, alpha, beta,
true, growth_type));
185 Teuchos::RCP<const Stokhos::ProductBasis<int,double> > basis;
192 bases, drop, use_old));
193 else if (prod_basis_type ==
TENSOR) {
207 else if (prod_basis_type ==
TOTAL) {
221 else if (prod_basis_type ==
SMOLYAK) {
226 bases, index_set, drop));
230 bases, index_set, drop));
234 bases, index_set, drop));
239 Teuchos::RCP<Cijk_type> Cijk;
241 Cijk = basis->computeTripleProductTensor();
243 Cijk = basis->computeLinearTripleProductTensor();
245 std::cout <<
"basis size = " << basis->size()
246 <<
" num nonzero Cijk entries = " << Cijk->
num_entries()
250 Epetra_MpiComm comm(MPI_COMM_WORLD);
252 Epetra_SerialComm comm;
256 Teuchos::RCP<Cijk_type> Cijk3 = Teuchos::rcp(
new Cijk_type);
258 Cijk_type::k_iterator k_begin = Cijk->k_begin();
259 Cijk_type::k_iterator k_end = Cijk->k_end();
260 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
262 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
263 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
264 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
266 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
267 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
268 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
271 if (i != 0 || (i == 0 &&
j == 0 && k == 0))
272 Cijk3->add_term(i,
j, k,
cijk);
277 Cijk3->fillComplete();
279 Teuchos::RCP<Epetra_CrsGraph> graph =
281 Epetra_CrsMatrix mat(Copy, *graph);
284 Ifpack_RCMReordering ifpack_rcm;
285 ifpack_rcm.SetParameter(
"reorder: root node", basis->size()-1);
286 ifpack_rcm.Compute(mat);
288 Teuchos::RCP<Cijk_type> Cijk2 = Teuchos::rcp(
new Cijk_type);
289 Cijk_type::k_iterator k_begin = Cijk->k_begin();
290 Cijk_type::k_iterator k_end = Cijk->k_end();
291 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
292 int k = ifpack_rcm.Reorder(
index(k_it));
293 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
294 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
295 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
296 int j = ifpack_rcm.Reorder(
index(j_it));
297 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
298 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
299 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
300 int i = ifpack_rcm.Reorder(
index(i_it));
302 Cijk2->add_term(i,
j, k,
cijk);
306 Cijk2->fillComplete();
311 std::cout << *Cijk << std::endl;
319 std::ofstream cijk_file(file_3tensor.c_str());
320 cijk_file.precision(14);
321 cijk_file.setf(std::ios::scientific);
322 cijk_file <<
"i, j, k, cijk" << std::endl;
323 Cijk_type::k_iterator k_begin = Cijk->k_begin();
324 Cijk_type::k_iterator k_end = Cijk->k_end();
325 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
327 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
328 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
329 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
331 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
332 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
333 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
336 if (!unique || ( i >=
j &&
j >= k ))
337 cijk_file << i <<
", " 340 <<
cijk << std::endl;
347 Teuchos::TimeMonitor::summarize(std::cout);
350 catch (std::exception& e) {
351 std::cout << e.what() << std::endl;
Hermite polynomial basis.
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 ProductBasisType prod_basis_type_values[]
const char * basis_type_names[]
void sparse3Tensor2MatrixMarket(const Stokhos::OrthogPolyBasis< ordinal_type, value_type > &basis, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, const Epetra_Comm &comm, const std::string &file)
const char * ordering_type_names[]
const char * growth_type_names[]
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
const BasisType basis_type_values[]
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
const int num_growth_types
const int num_prod_basis_types
Legendre polynomial basis using Gauss-Patterson quadrature points.
const int num_ordering_types
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.
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Multivariate orthogonal polynomial basis generated from a tensor product of univariate polynomials...
Legendre polynomial basis.
Stokhos::Sparse3Tensor< int, double > Cijk_type
int main(int argc, char **argv)
An isotropic total order index set.
A comparison functor implementing a strict weak ordering based Morton Z-ordering. ...
Legendre polynomial basis using Clenshaw-Curtis quadrature points.
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.
SparseArrayIterator< index_iterator, value_iterator >::value_reference value(const SparseArrayIterator< index_iterator, value_iterator > &it)
const char * prod_basis_type_names[]
A comparison functor implementing a strict weak ordering based lexographic ordering.
const Stokhos::GrowthPolicy growth_type_values[]
const OrderingType ordering_type_values[]
const int num_basis_types