42 #ifndef STOKHOS_LTB_SPARSE_3_TENSOR_HPP 43 #define STOKHOS_LTB_SPARSE_3_TENSOR_HPP 56 template <
typename ordinal_type,
typename value_type>
62 Teuchos::Array< Teuchos::RCP<CijkNode> >
children;
85 void print(std::ostream& os)
const {}
93 if (
head != Teuchos::null)
94 return head->total_num_entries;
100 if (
head != Teuchos::null)
101 return head->total_num_leafs;
126 template <
typename ordinal_type,
typename value_type>
134 template <
typename ordinal_type>
136 Teuchos::Array< Teuchos::RCP<LexicographicTreeBasisNode> >
children;
137 Teuchos::Array< Stokhos::MultiIndex<ordinal_type> >
terms;
170 template <
typename ordinal_type>
171 Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> >
173 const Teuchos::ArrayView<const ordinal_type>& basis_orders,
187 Teuchos::RCP<node_type> node = Teuchos::rcp(
new node_type);
188 node->index_begin = index_begin;
189 node->terms.resize(my_p+1);
191 node->terms[i].resize(prev_d+1);
193 node->terms[i][
j] = term_prefix[
j];
194 node->terms[i][prev_d] = i;
199 node->block_size = my_p+1;
204 Teuchos::ArrayView<const ordinal_type> bo =
205 Teuchos::arrayView(&basis_orders[1],my_d-1);
206 node->children.resize(my_p+1);
207 node->index_begin = index_begin;
208 node->block_size = 0;
211 bo, total_order, index_begin+node->block_size, order_sum+i,
213 node->children[i] = child;
214 node->block_size += child->block_size;
225 Teuchos::RCP< LTBSparse3Tensor<ordinal_type, value_type> >
228 bool symmetric =
false) {
229 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 230 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Triple-Product Tensor Time");
234 using Teuchos::Array;
235 using Teuchos::ArrayView;
237 const Array< RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases = product_basis.getCoordinateBases();
240 Array<ordinal_type> basis_orders(d);
241 for (
int i=0; i<d; ++i)
242 basis_orders[i] = bases[i]->order();
243 ArrayView<const ordinal_type> basis_orders_view = basis_orders();
246 Array< RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
249 bases[i]->computeSparseTripleProductTensor(bases[i]->order()+1);
251 ArrayView<const RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d_view
255 Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> > ltb =
260 RCP<Cijk_type> Cijk = rcp(
new Cijk_type(symmetric));
261 RCP<typename Cijk_type::CijkNode> head =
263 basis_orders_view, Cijk_1d_view, ltb, ltb, ltb, p, symmetric);
264 Cijk->setHeadNode(head);
271 Teuchos::RCP<typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode>
273 const Teuchos::ArrayView<const ordinal_type>& basis_orders,
279 const bool symmetric,
287 using Teuchos::ArrayView;
288 using Teuchos::arrayView;
297 Cijk_Iterator cijk_iterator(p_i, p_j, p_k, symmetric);
299 RCP<node_type> node = rcp(
new node_type);
303 node->i_begin = i_ltb->index_begin;
304 node->j_begin = j_ltb->index_begin;
305 node->k_begin = k_ltb->index_begin;
306 node->i_size = i_ltb->block_size;
307 node->j_size = j_ltb->block_size;
308 node->k_size = k_ltb->block_size;
312 node->is_leaf =
true;
316 Cijk_1d[0]->getValue(cijk_iterator.i,
319 node->values.push_back(
cijk*cijk_base);
320 more = cijk_iterator.increment();
322 node->my_num_entries = node->values.size();
323 node->total_num_entries = node->values.size();
324 node->total_num_leafs = 1;
329 node->is_leaf =
false;
330 ArrayView<const ordinal_type> bo = arrayView(&basis_orders[1], my_d-1);
331 ArrayView<const RCP<Sparse3Tensor<ordinal_type,value_type> > > c1d =
332 arrayView(&Cijk_1d[1], my_d-1);
333 node->total_num_entries = 0;
334 node->total_num_leafs = 0;
338 Cijk_1d[0]->getValue(cijk_iterator.i,
341 RCP<node_type> child =
343 i_ltb->children[cijk_iterator.i],
344 j_ltb->children[cijk_iterator.j],
345 k_ltb->children[cijk_iterator.k],
346 total_order, symmetric,
347 sum_i + cijk_iterator.i,
348 sum_j + cijk_iterator.j,
349 sum_k + cijk_iterator.k,
351 node->children.push_back(child);
352 node->total_num_entries += child->total_num_entries;
353 node->total_num_leafs += child->total_num_leafs;
354 more = cijk_iterator.increment();
356 node->my_num_entries = node->children.size();
367 Teuchos::RCP< LTBSparse3Tensor<ordinal_type, value_type> >
370 bool symmetric =
false) {
371 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 372 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Total Triple-Product Tensor Time");
376 using Teuchos::Array;
377 using Teuchos::ArrayView;
379 const Array< RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases = product_basis.getCoordinateBases();
382 Array<ordinal_type> basis_orders(d);
383 for (
int i=0; i<d; ++i)
384 basis_orders[i] = bases[i]->order();
385 ArrayView<const ordinal_type> basis_orders_view = basis_orders();
388 Array< RCP<Dense3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
390 Cijk_1d[i] = bases[i]->computeTripleProductTensor();
392 ArrayView<const RCP<Dense3Tensor<ordinal_type,value_type> > > Cijk_1d_view
396 Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> > ltb =
401 RCP<Cijk_type> Cijk = rcp(
new Cijk_type(symmetric));
402 RCP<typename Cijk_type::CijkNode> head =
404 basis_orders_view, Cijk_1d_view, ltb, ltb, ltb, p, symmetric);
405 Cijk->setHeadNode(head);
412 Teuchos::RCP<typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode>
414 const Teuchos::ArrayView<const ordinal_type>& basis_orders,
420 const bool symmetric,
425 const bool parent_j_equals_k =
true) {
429 using Teuchos::ArrayView;
430 using Teuchos::arrayView;
439 RCP<node_type> node = rcp(
new node_type);
443 node->i_begin = i_ltb->index_begin;
444 node->j_begin = j_ltb->index_begin;
445 node->k_begin = k_ltb->index_begin;
446 node->i_size = i_ltb->block_size;
447 node->j_size = j_ltb->block_size;
448 node->k_size = k_ltb->block_size;
449 node->parent_j_equals_k = parent_j_equals_k;
456 node->is_leaf =
true;
457 node->values.reserve((p_i+1)*(p_j+1)*(p_k+1));
466 if (symmetric && (k0%2 != (i+
j)%2)) ++k0;
471 if (
j+node->j_begin == k+node->k_begin)
cijk *= 0.5;
472 node->values.push_back(
cijk*cijk_base);
476 node->my_num_entries = node->values.size();
477 node->total_num_entries = node->values.size();
478 node->total_num_leafs = 1;
483 node->is_leaf =
false;
484 ArrayView<const ordinal_type> bo = arrayView(&basis_orders[1], my_d-1);
485 ArrayView<const RCP<Dense3Tensor<ordinal_type,value_type> > > c1d =
486 arrayView(&Cijk_1d[1], my_d-1);
487 node->total_num_entries = 0;
488 node->total_num_leafs = 0;
493 if (symmetric && (k0%2 != (i+
j)%2)) ++k0;
498 RCP<node_type> child =
503 total_order, symmetric,
508 parent_j_equals_k &&
j == k);
509 node->children.push_back(child);
510 node->total_num_entries += child->total_num_entries;
511 node->total_num_leafs += child->total_num_leafs;
515 node->my_num_entries = node->children.size();
521 template <
typename ordinal_type>
528 template <
typename ordinal_type,
typename value_type>
543 template <
typename ordinal_type,
typename value_type>
544 Teuchos::RCP< FlatLTBSparse3Tensor<ordinal_type,value_type> >
547 bool symmetric =
false) {
548 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR 549 TEUCHOS_FUNC_TIME_MONITOR(
"Stokhos: Flat Triple-Product Tensor Time");
556 RCP<Cijk_LTB_type> ltb_tensor =
560 RCP< FlatLTBSparse3Tensor<ordinal_type,value_type> > flat_tensor =
562 flat_tensor->node.reserve(ltb_tensor->num_leafs());
563 flat_tensor->cijk.reserve(ltb_tensor->num_entries());
564 flat_tensor->symmetric = symmetric;
567 typedef typename Cijk_LTB_type::CijkNode
node_type;
568 Teuchos::Array< Teuchos::RCP<const node_type> > node_stack;
569 Teuchos::Array< ordinal_type > index_stack;
570 node_stack.push_back(ltb_tensor->getHeadNode());
571 index_stack.push_back(0);
572 Teuchos::RCP<const node_type> node;
574 while (node_stack.size() > 0) {
575 node = node_stack.back();
576 child_index = index_stack.back();
584 leaf.
p_i = node->p_i;
585 leaf.
p_j = node->p_j;
586 leaf.
p_k = node->p_k;
588 flat_tensor->node.push_back(leaf);
589 flat_tensor->cijk.insert(flat_tensor->cijk.end(),
590 node->values.begin(),
592 node_stack.pop_back();
593 index_stack.pop_back();
597 else if (child_index < node->children.size()) {
598 ++index_stack.back();
599 node = node->children[child_index];
600 node_stack.push_back(node);
601 index_stack.push_back(0);
606 node_stack.pop_back();
607 index_stack.pop_back();
615 template <
int max_size,
typename ordinal_type,
typename value_type>
633 typedef typename cijk_type::node_const_iterator node_iterator;
634 typedef typename cijk_type::cijk_const_iterator cijk_iterator;
635 node_iterator ni =
cijk.node.begin();
636 node_iterator ni_end =
cijk.node.end();
637 cijk_iterator ci =
cijk.cijk.begin();
638 for (; ni != ni_end; ++ni) {
640 const value_type *a_j_block = ca + ni->j_begin;
641 const value_type *b_k_block = cb + ni->k_begin;
642 const value_type *a_k_block = ca + ni->k_begin;
643 const value_type *b_j_block = cb + ni->j_begin;
649 ab[
j][k] = a_j_block[
j]*b_k_block[k] + a_k_block[k]*b_j_block[
j];
661 if (
cijk.symmetric && (k0%2 != (i+
j)%2)) ++k0;
665 tmp += (*ci)*ab[
j][k];
679 #endif // STOKHOS_LTB_SPARSE_3_TENSOR_HPP
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
ordinal_type total_num_leafs
void init(const value_type &v)
Initialize coefficients to value.
LTBSparse3Tensor(const bool symm)
Constructor.
cijk_array_type::iterator cijk_iterator
Teuchos::Array< FlatLTBSparse3TensorNode< ordinal_type > > node_array_type
LexicographicTreeBasisNode(const LexicographicTreeBasisNode &node)
LexicographicTreeBasisNode()
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
pointer coeff()
Return coefficient array.
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTBBlockLeaf(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
~LTBSparse3Tensor()
Destructor.
ordinal_type my_num_entries
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
LexicographicTreeBasisNode & operator=(const LexicographicTreeBasisNode &node)
node_array_type::const_iterator node_const_iterator
cijk_array_type::const_iterator cijk_const_iterator
Data structure storing a sparse 3-tensor C(i,j,k) in a a tree-based format for lexicographically orde...
node_array_type::iterator node_iterator
Teuchos::RCP< LexicographicTreeBasisNode< ordinal_type > > build_lexicographic_basis_tree(const Teuchos::ArrayView< const ordinal_type > &basis_orders, const ordinal_type total_order, const ordinal_type index_begin=ordinal_type(0), const ordinal_type order_sum=ordinal_type(0), const Stokhos::MultiIndex< ordinal_type > &term_prefix=Stokhos::MultiIndex< ordinal_type >())
Teuchos::RCP< FlatLTBSparse3Tensor< ordinal_type, value_type > > computeFlatTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
std::ostream & operator<<(std::ostream &os, const ProductContainer< coeff_type > &vec)
Teuchos::Array< Teuchos::RCP< LexicographicTreeBasisNode > > children
void flatLTB3TensorMultiply(OrthogPolyApprox< ordinal_type, value_type > &c, const OrthogPolyApprox< ordinal_type, value_type > &a, const OrthogPolyApprox< ordinal_type, value_type > &b, const FlatLTBSparse3Tensor< ordinal_type, value_type > &cijk)
bool symmetric() const
Return if symmetric.
KOKKOS_INLINE_FUNCTION PCE< Storage > max(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
Top-level namespace for Stokhos classes and functions.
Teuchos::RCP< typename LTBSparse3Tensor< ordinal_type, value_type >::CijkNode > computeCijkLTBNode(const Teuchos::ArrayView< const ordinal_type > &basis_orders, const Teuchos::ArrayView< const Teuchos::RCP< Sparse3Tensor< ordinal_type, value_type > > > &Cijk_1d, const Teuchos::RCP< LexicographicTreeBasisNode< ordinal_type > > &i_ltb, const Teuchos::RCP< LexicographicTreeBasisNode< ordinal_type > > &j_ltb, const Teuchos::RCP< LexicographicTreeBasisNode< ordinal_type > > &k_ltb, const ordinal_type total_order, const bool symmetric, const ordinal_type sum_i=ordinal_type(0), const ordinal_type sum_j=ordinal_type(0), const ordinal_type sum_k=ordinal_type(0), const value_type cijk_base=value_type(1))
Teuchos::RCP< typename LTBSparse3Tensor< ordinal_type, value_type >::CijkNode > computeCijkLTBNodeBlockLeaf(const Teuchos::ArrayView< const ordinal_type > &basis_orders, const Teuchos::ArrayView< const Teuchos::RCP< Dense3Tensor< ordinal_type, value_type > > > &Cijk_1d, const Teuchos::RCP< LexicographicTreeBasisNode< ordinal_type > > &i_ltb, const Teuchos::RCP< LexicographicTreeBasisNode< ordinal_type > > &j_ltb, const Teuchos::RCP< LexicographicTreeBasisNode< ordinal_type > > &k_ltb, const ordinal_type total_order, const bool symmetric, const ordinal_type sum_i=ordinal_type(0), const ordinal_type sum_j=ordinal_type(0), const ordinal_type sum_k=ordinal_type(0), const value_type cijk_base=value_type(1), const bool parent_j_equals_k=true)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
value_type getValue(ordinal_type i, ordinal_type j, ordinal_type k) const
Get Cijk value for a given i, j, k indices.
Teuchos::RCP< CijkNode > head
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)
Stokhos::Sparse3Tensor< int, double > Cijk_type
Node type used in constructing the tree.
ordinal_type num_entries() const
Return number of non-zero entries.
Teuchos::RCP< const CijkNode > getHeadNode() const
Get the head node.
void print(std::ostream &os) const
Print tensor.
A comparison functor implementing a strict weak ordering based lexographic ordering.
Data structure storing a dense 3-tensor C(i,j,k).
ordinal_type total_num_entries
Teuchos::Array< value_type > values
void setHeadNode(const Teuchos::RCP< CijkNode > &h)
Set the head node.
Teuchos::Array< value_type > cijk_array_type
Teuchos::Array< Teuchos::RCP< CijkNode > > children
LTBSparse3Tensor & operator=(const LTBSparse3Tensor &b)
ordinal_type num_leafs() const
Return number of nodes.