80 value_type> > >&
bases,
82 const coeff_compare_type& coeff_compare = coeff_compare_type());
91 ordinal_type
order()
const;
97 virtual ordinal_type
size()
const;
104 virtual const Teuchos::Array<value_type>&
norm_squared()
const;
107 virtual const value_type&
norm_squared(ordinal_type i)
const;
117 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
122 Teuchos::RCP< Stokhos::Sparse3Tensor<ordinal_type, value_type> >
134 const Teuchos::ArrayView<const value_type>& point,
135 Teuchos::Array<value_type>& basis_vals)
const;
138 virtual void print(std::ostream& os)
const;
141 virtual const std::string&
getName()
const;
154 virtual const MultiIndex<ordinal_type>&
term(ordinal_type i)
const;
161 virtual ordinal_type
index(
const MultiIndex<ordinal_type>&
term)
const;
203 Teuchos::Array< Teuchos::RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >
bases;
239 bool symmetric =
false) {
240#ifdef STOKHOS_TEUCHOS_TIME_MONITOR
241 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Triple-Product Tensor Time");
245 using Teuchos::Array;
247 typedef MultiIndex<ordinal_type> coeff_type;
248 const Array< RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases = product_basis.getCoordinateBases();
249 ordinal_type d = bases.size();
251 Array<ordinal_type> basis_orders(d);
252 for (
int i=0; i<d; ++i)
253 basis_orders[i] = bases[i]->order();
256 Array< RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
257 for (ordinal_type i=0; i<d; i++) {
259 bases[i]->computeSparseTripleProductTensor(bases[i]->order()+1);
262 RCP< Sparse3Tensor<ordinal_type, value_type> > Cijk =
263 rcp(
new Sparse3Tensor<ordinal_type, value_type>);
267 Array<Cijk_Iterator> Cijk_1d_iterators(d);
268 coeff_type terms_i(d,0), terms_j(d,0), terms_k(d,0);
269 Array<ordinal_type> sum_i(d,0), sum_j(d,0), sum_k(d,0);
270 for (ordinal_type dim=0; dim<d; dim++) {
271 Cijk_1d_iterators[dim] = Cijk_Iterator(bases[dim]->order(), symmetric);
277 ordinal_type cnt = 0;
282 for (ordinal_type dim=0; dim<d; ++dim) {
283 terms_i[dim] = Cijk_1d_iterators[dim].i;
284 terms_j[dim] = Cijk_1d_iterators[dim].j;
285 terms_k[dim] = Cijk_1d_iterators[dim].k;
302 value_type c = value_type(1.0);
303 for (ordinal_type dim=0; dim<d; dim++) {
304 c *= Cijk_1d[dim]->getValue(Cijk_1d_iterators[dim].i,
305 Cijk_1d_iterators[dim].
j,
306 Cijk_1d_iterators[dim].k);
309 TEUCHOS_TEST_FOR_EXCEPTION(
310 std::abs(c) <= 1.0e-12,
312 "Got 0 triple product value " << c
313 <<
", I = " << I <<
" = " << terms_i
314 <<
", J = " << J <<
" = " << terms_j
315 <<
", K = " << K <<
" = " << terms_k
319 Cijk->add_term(I,J,K,c);
327 ordinal_type cdim = d-1;
329 while (inc && cdim >= 0) {
330 ordinal_type delta_i, delta_j, delta_k;
332 Cijk_1d_iterators[cdim].increment(delta_i, delta_j, delta_k);
342 for (ordinal_type ii=0; ii<delta_i; ++ii)
345 basis_orders[cdim+1]-sum_i[cdim] -
346 (Cijk_1d_iterators[cdim].i-ii)+d-cdim,
350 for (ordinal_type ii=0; ii<-delta_i; ++ii)
353 basis_orders[cdim+1]-sum_i[cdim] -
354 (Cijk_1d_iterators[cdim].i+ii)+d-cdim-1,
359 for (ordinal_type jj=0; jj<delta_j; ++jj)
362 basis_orders[cdim+1]-sum_j[cdim] -
363 (Cijk_1d_iterators[cdim].
j-jj)+d-cdim,
367 for (ordinal_type jj=0; jj<-delta_j; ++jj)
370 basis_orders[cdim+1]-sum_j[cdim] -
371 (Cijk_1d_iterators[cdim].
j+jj)+d-cdim-1,
376 for (ordinal_type kk=0; kk<delta_k; ++kk)
379 basis_orders[cdim+1]-sum_k[cdim] -
380 (Cijk_1d_iterators[cdim].k-kk)+d-cdim,
384 for (ordinal_type kk=0; kk<-delta_k; ++kk)
387 basis_orders[cdim+1]-sum_k[cdim] -
388 (Cijk_1d_iterators[cdim].k+kk)+d-cdim-1,
402 for (ordinal_type dim=cdim+1; dim<d; ++dim) {
405 sum_i[dim] = sum_i[dim-1] + Cijk_1d_iterators[dim-1].i;
406 sum_j[dim] = sum_j[dim-1] + Cijk_1d_iterators[dim-1].j;
407 sum_k[dim] = sum_k[dim-1] + Cijk_1d_iterators[dim-1].k;
410 Cijk_1d_iterators[dim] =
411 Cijk_Iterator(basis_orders[dim]-sum_i[dim],
412 basis_orders[dim]-sum_j[dim],
413 basis_orders[dim]-sum_k[dim],
425 Cijk->fillComplete();