42#include "Teuchos_TimeMonitor.hpp"
43#include "Teuchos_TestForException.hpp"
45template <
typename ordinal_type,
typename value_type,
typename ordering_type>
46template <
typename index_set_type>
49 const Teuchos::Array< Teuchos::RCP<
const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases_,
50 const index_set_type& index_set,
51 const value_type& sparse_tol_,
52 const ordering_type& coeff_compare) :
57 sparse_tol(sparse_tol_),
59 basis_set(coeff_compare),
90 typedef std::map<multiindex_type,ordinal_type,index_compare> index_map_type;
91 ordinal_type dim = index_set.dimension();
92 alpha_set_type alpha_set(dim, 1);
93 typename alpha_set_type::iterator alpha_begin = alpha_set.begin();
94 typename alpha_set_type::iterator alpha_end = alpha_set.end();
95 typename index_set_type::iterator index_iterator = index_set.begin();
96 typename index_set_type::iterator index_end = index_set.end();
98 index_map_type index_map;
99 for (; index_iterator != index_end; ++index_iterator) {
100 for (
typename alpha_set_type::iterator alpha = alpha_begin;
101 alpha != alpha_end; ++alpha) {
102 bool valid_index =
true;
103 for (ordinal_type i=0; i<dim; ++i) {
104 diff[i] = (*index_iterator)[i] - (*alpha)[i];
111 ordinal_type alpha_order = alpha->order();
113 if (alpha_order % 2 == 0)
117 typename index_map_type::iterator index_map_iterator =
118 index_map.find(diff);
119 if (index_map_iterator == index_map.end())
120 index_map[diff] =
val;
122 index_map_iterator->second +=
val;
128 typename index_map_type::iterator index_map_iterator = index_map.begin();
129 typename index_map_type::iterator index_map_end = index_map.end();
130 for (; index_map_iterator != index_map_end; ++index_map_iterator) {
133 if (index_map_iterator->second == 0)
138 for (ordinal_type i=0; i<dim; ++i) {
139 coeff_growth_index[i] =
140 bases[i]->coefficientGrowth(index_map_iterator->first[i]);
144 Teuchos::RCP<tensor_product_basis_type> tp =
149 for (ordinal_type i=0; i<tp->size(); ++i)
150 basis_set[tp->term(i)] = ordinal_type(0);
154 TensorProductPredicate<ordinal_type>(coeff_growth_index) );
160 ordinal_type idx = 0;
162 for (
typename coeff_set_type::iterator i =
basis_set.begin();
171 for (ordinal_type i=0; i<
sz; ++i) {
172 for (ordinal_type
j=0;
j<dim; ++
j)
178 for (ordinal_type i=0; i<dim; i++)
184 for (ordinal_type i=0; i<
d; i++) {
192 for (ordinal_type k=0; k<
sz; k++) {
193 nrm = value_type(1.0);
194 for (ordinal_type i=0; i<
d; i++)
200 name =
"Smolyak basis (";
201 for (ordinal_type i=0; i<
d-1; i++)
207 for (ordinal_type
j=0;
j<
d;
j++)
211template <
typename ordinal_type,
typename value_type,
typename ordering_type>
217template <
typename ordinal_type,
typename value_type,
typename ordering_type>
225template <
typename ordinal_type,
typename value_type,
typename ordering_type>
233template <
typename ordinal_type,
typename value_type,
typename ordering_type>
241template <
typename ordinal_type,
typename value_type,
typename ordering_type>
242const Teuchos::Array<value_type>&
249template <
typename ordinal_type,
typename value_type,
typename ordering_type>
257template <
typename ordinal_type,
typename value_type,
typename ordering_type>
258Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
262#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
263 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Triple-Product Tensor Fill Time");
267 bases, basis_set, basis_map, sm_pred, sm_pred, sparse_tol);
270template <
typename ordinal_type,
typename value_type,
typename ordering_type>
271Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
275#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
276 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Triple-Product Tensor Fill Time");
279 SmolyakPredicate< TotalOrderPredicate<ordinal_type> > k_pred;
280 for (ordinal_type i=0; i<sm_pred.tp_preds.size(); ++i) {
281 k_pred.tp_preds.push_back(
282 TotalOrderPredicate<ordinal_type>(1, sm_pred.tp_preds[i].orders) );
286 bases, basis_set, basis_map, sm_pred, k_pred, sparse_tol);
289template <
typename ordinal_type,
typename value_type,
typename ordering_type>
296 value_type z = value_type(1.0);
297 for (ordinal_type
j=0;
j<d;
j++)
298 z = z * bases[
j]->evaluate(value_type(0.0), basis_map[i][
j]);
303template <
typename ordinal_type,
typename value_type,
typename ordering_type>
306evaluateBases(
const Teuchos::ArrayView<const value_type>& point,
307 Teuchos::Array<value_type>& basis_vals)
const
309 for (ordinal_type
j=0;
j<d;
j++)
310 bases[
j]->evaluateBases(point[
j], basis_eval_tmp[
j]);
313 for (ordinal_type i=0; i<sz; i++) {
314 value_type t = value_type(1.0);
315 for (ordinal_type
j=0;
j<d;
j++)
316 t *= basis_eval_tmp[
j][basis_map[i][
j]];
321template <
typename ordinal_type,
typename value_type,
typename ordering_type>
324print(std::ostream& os)
const
326 os <<
"Smolyak basis of order " << p <<
", dimension " << d
327 <<
", and size " << sz <<
". Component bases:\n";
328 for (ordinal_type i=0; i<d; i++)
330 os <<
"Basis vector norms (squared):\n\t";
331 for (ordinal_type i=0; i<static_cast<ordinal_type>(norms.size()); i++)
332 os << norms[i] <<
" ";
336template <
typename ordinal_type,
typename value_type,
typename ordering_type>
339term(ordinal_type i)
const
344template <
typename ordinal_type,
typename value_type,
typename ordering_type>
349 typename coeff_set_type::const_iterator it = basis_set.find(term);
350 TEUCHOS_TEST_FOR_EXCEPTION(it == basis_set.end(), std::logic_error,
351 "Invalid term " << term);
355template <
typename ordinal_type,
typename value_type,
typename ordering_type>
363template <
typename ordinal_type,
typename value_type,
typename ordering_type>
364Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<ordinal_type, value_type> > >
371template <
typename ordinal_type,
typename value_type,
typename ordering_type>
A comparison functor implementing a strict weak ordering based lexographic ordering.
void resize(ordinal_type d, ordinal_type v=ordinal_type(0))
Resize.
static Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor(const Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > &bases, const basis_set_type &basis_set, const basis_map_type &basis_map, const coeff_predicate_type &coeff_pred, const k_coeff_predicate_type &k_coeff_pred, const value_type sparse_tol=1.0e-12)
Teuchos::Array< Teuchos::RCP< tensor_product_basis_type > > tp_bases
Tensor product bases comprising Smolyak set.
virtual void evaluateBases(const Teuchos::ArrayView< const value_type > &point, Teuchos::Array< value_type > &basis_vals) const
Evaluate basis polynomials at given point point.
virtual const Teuchos::Array< value_type > & norm_squared() const
Return array storing norm-squared of each basis polynomial.
ordinal_type dimension() const
Return dimension of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeLinearTripleProductTensor() const
Compute linear triple product tensor where k = 0,1,..,d.
virtual ordinal_type index(const MultiIndex< ordinal_type > &term) const
Get index of the multivariate polynomial given orders of each coordinate.
Teuchos::Array< Teuchos::Array< value_type > > basis_eval_tmp
Temporary array used in basis evaluation.
Teuchos::Array< ordinal_type > smolyak_coeffs
Smolyak coefficients.
TensorProductBasis< ordinal_type, value_type, LexographicLess< coeff_type > > tensor_product_basis_type
std::string name
Name of basis.
virtual Teuchos::RCP< Stokhos::Sparse3Tensor< ordinal_type, value_type > > computeTripleProductTensor() const
Compute triple product tensor.
ordinal_type d
Total dimension of basis.
ordinal_type order() const
Return order of basis.
coeff_map_type basis_map
Basis map.
virtual ordinal_type size() const
Return total size of basis.
coeff_type max_orders
Maximum orders for each dimension.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > bases
Array of bases.
ordinal_type sz
Total size of basis.
value_type sparse_tol
Tolerance for computing sparse Cijk.
coeff_set_type basis_set
Basis set.
virtual const std::string & getName() const
Return string name of basis.
Teuchos::Array< value_type > norms
Norms.
virtual MultiIndex< ordinal_type > getMaxOrders() const
Return maximum order allowable for each coordinate basis.
SmolyakPredicate< TensorProductPredicate< ordinal_type > > sm_pred
Predicate for building sparse triple products.
Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const
Return coordinate bases.
SmolyakBasis(const Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > &bases, const index_set_type &index_set, const value_type &sparse_tol=1.0e-12, const coeff_compare_type &coeff_compare=coeff_compare_type())
Constructor.
virtual const MultiIndex< ordinal_type > & term(ordinal_type i) const
Get orders of each coordinate polynomial given an index i.
virtual void print(std::ostream &os) const
Print basis to stream os.
virtual ~SmolyakBasis()
Destructor.
ordinal_type p
Total order of basis.
virtual value_type evaluateZero(ordinal_type i) const
Evaluate basis polynomial i at zero.
A tensor product index set.