45 #include "Teuchos_CommandLineProcessor.hpp" 59 "complete",
"tensor",
"total",
"smolyak" };
67 "total",
"lexicographic" };
77 "none",
"2-way",
"6-way" };
82 using Teuchos::ArrayView;
83 using Teuchos::ArrayRCP;
90 Teuchos::CommandLineProcessor
CLP;
92 "This example generates the sparsity pattern for the block stochastic Galerkin matrix.\n");
94 CLP.setOption(
"dimension", &d,
"Stochastic dimension");
96 CLP.setOption(
"order", &p,
"Polynomial order");
97 double drop = 1.0e-12;
98 CLP.setOption(
"drop", &drop,
"Drop tolerance");
99 bool symmetric =
true;
100 CLP.setOption(
"symmetric",
"asymmetric", &symmetric,
"Use basis polynomials with symmetric PDF");
102 CLP.setOption(
"growth", &growth_type,
106 CLP.setOption(
"product_basis", &prod_basis_type,
109 "Product basis type");
111 CLP.setOption(
"ordering", &ordering_type,
114 "Product basis ordering");
116 CLP.setOption(
"symmetry", &symmetry_type,
119 "Cijk symmetry type");
121 CLP.setOption(
"full",
"linear", &
full,
"Use full or linear expansion");
123 CLP.setOption(
"tile_size", &tile_size,
"Tile size");
124 bool save_3tensor =
false;
125 CLP.setOption(
"save_3tensor",
"no-save_3tensor", &save_3tensor,
126 "Save full 3tensor to file");
127 std::string file_3tensor =
"Cijk.dat";
128 CLP.setOption(
"filename_3tensor", &file_3tensor,
129 "Filename to store full 3-tensor");
135 Array< RCP<const Stokhos::OneDOrthogPolyBasis<int,double> > > bases(d);
136 const double alpha = 1.0;
137 const double beta = symmetric ? 1.0 : 2.0 ;
138 for (
int i=0; i<d; i++) {
140 p, alpha, beta,
true, growth_type));
142 RCP<const Stokhos::ProductBasis<int,double> > basis;
149 else if (prod_basis_type ==
TENSOR) {
159 else if (prod_basis_type ==
TOTAL) {
169 else if (prod_basis_type ==
SMOLYAK) {
174 bases, index_set, drop));
178 bases, index_set, drop));
185 Cijk = basis->computeTripleProductTensor();
187 Cijk = basis->computeLinearTripleProductTensor();
189 int basis_size = basis->size();
190 std::cout <<
"basis size = " << basis_size
191 <<
" num nonzero Cijk entries = " << Cijk->
num_entries()
195 std::ofstream cijk_file;
197 cijk_file.open(file_3tensor.c_str());
198 cijk_file.precision(14);
199 cijk_file.setf(std::ios::scientific);
200 cijk_file <<
"i, j, k, part" << std::endl;
204 Teuchos::ArrayRCP< Stokhos::CijkData<int,double> > coordinate_list =
209 rcb_type rcb(tile_size, 10000, coordinate_list());
210 int num_parts = rcb.get_num_parts();
211 std::cout <<
"num parts = " << num_parts << std::endl;
214 RCP< Array< RCP<rcb_type::Box> > > parts = rcb.get_parts();
215 for (
int i=0; i<num_parts; ++i) {
216 RCP<rcb_type::Box> box = (*parts)[i];
217 std::cout <<
"part " << i <<
" bounding box =" 218 <<
" [ " << box->delta_x <<
", " << box->delta_y <<
", " 219 << box->delta_z <<
" ]" <<
" nnz = " 220 << box->coords.size() << std::endl;
224 ArrayRCP<int> part_ids = rcb.get_part_IDs();
226 Cijk_type::k_iterator k_begin = Cijk->k_begin();
227 Cijk_type::k_iterator k_end = Cijk->k_end();
228 for (Cijk_type::k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
230 Cijk_type::kj_iterator j_begin = Cijk->j_begin(k_it);
231 Cijk_type::kj_iterator j_end = Cijk->j_end(k_it);
232 for (Cijk_type::kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
234 Cijk_type::kji_iterator i_begin = Cijk->i_begin(j_it);
235 Cijk_type::kji_iterator i_end = Cijk->i_end(j_it);
236 for (Cijk_type::kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
241 cijk_file << i <<
", " <<
j <<
", " << k <<
", " 242 << part_ids[idx++] << std::endl;
253 catch (std::exception& e) {
254 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...
const char * ordering_type_names[]
const char * symmetry_type_names[]
GrowthPolicy
Enumerated type for determining Smolyak growth policies.
A comparison functor implementing a strict weak ordering based total-order ordering, recursive on the dimension.
const OrderingType ordering_type_values[]
const ProductBasisType prod_basis_type_values[]
const Stokhos::CijkSymmetryType symmetry_type_values[]
const int num_symmetry_types
const Stokhos::GrowthPolicy growth_type_values[]
const int num_ordering_types
ordinal_type num_entries() const
Return number of non-zero entries.
const int num_growth_types
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.
const char * prod_basis_type_names[]
int main(int argc, char **argv)
A comparison functor implementing a strict weak ordering based lexographic ordering.
const char * growth_type_names[]
Teuchos::ArrayRCP< CijkData< ordinal_type, scalar_type > > build_cijk_coordinate_list(const Sparse3Tensor< ordinal_type, scalar_type > &Cijk, CijkSymmetryType symmetry_type)
const int num_prod_basis_types