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 {}
88 value_type
getValue(ordinal_type i, ordinal_type
j, ordinal_type k)
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;
149 for (ordinal_type i=0; i<
children.size(); ++i)
159 for (ordinal_type i=0; i<
children.size(); ++i)
170 template <
typename ordinal_type>
171 Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> >
173 const Teuchos::ArrayView<const ordinal_type>& basis_orders,
174 const ordinal_type total_order,
175 const ordinal_type index_begin = ordinal_type(0),
176 const ordinal_type order_sum = ordinal_type(0),
182 ordinal_type my_d = basis_orders.size();
183 ordinal_type prev_d = term_prefix.dimension();
184 ordinal_type p = basis_orders[0];
185 ordinal_type my_p = std::min(total_order-order_sum, p);
187 Teuchos::RCP<node_type> node = Teuchos::rcp(
new node_type);
188 node->index_begin = index_begin;
189 node->terms.resize(my_p+1);
190 for (ordinal_type i=0; i<=my_p; ++i) {
191 node->terms[i].resize(prev_d+1);
192 for (ordinal_type
j=0;
j<prev_d; ++
j)
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;
209 for (ordinal_type i=0; i<=my_p; ++i) {
211 bo, total_order, index_begin+node->block_size, order_sum+i,
213 node->children[i] = child;
214 node->block_size += child->block_size;
223 template <
typename ordinal_type,
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();
238 ordinal_type d = bases.size();
239 ordinal_type p = product_basis.order();
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);
247 for (ordinal_type i=0; i<d; i++) {
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);
269 template <
typename ordinal_type,
271 Teuchos::RCP<typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode>
273 const Teuchos::ArrayView<const ordinal_type>& basis_orders,
278 const ordinal_type total_order,
279 const bool symmetric,
280 const ordinal_type sum_i = ordinal_type(0),
281 const ordinal_type sum_j = ordinal_type(0),
282 const ordinal_type sum_k = ordinal_type(0),
283 const value_type cijk_base = value_type(1)) {
287 using Teuchos::ArrayView;
288 using Teuchos::arrayView;
292 ordinal_type my_d = basis_orders.size();
293 ordinal_type p = basis_orders[0];
294 ordinal_type p_i = std::min(total_order-sum_i, p);
295 ordinal_type p_j = std::min(total_order-sum_j, p);
296 ordinal_type p_k = std::min(total_order-sum_k, p);
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();
365 template <
typename ordinal_type,
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();
380 ordinal_type d = bases.size();
381 ordinal_type p = product_basis.order();
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);
389 for (ordinal_type i=0; i<d; i++) {
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);
410 template <
typename ordinal_type,
412 Teuchos::RCP<typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode>
414 const Teuchos::ArrayView<const ordinal_type>& basis_orders,
419 const ordinal_type total_order,
420 const bool symmetric,
421 const ordinal_type sum_i = ordinal_type(0),
422 const ordinal_type sum_j = ordinal_type(0),
423 const ordinal_type sum_k = ordinal_type(0),
424 const value_type cijk_base = value_type(1),
425 const bool parent_j_equals_k =
true) {
429 using Teuchos::ArrayView;
430 using Teuchos::arrayView;
433 ordinal_type my_d = basis_orders.size();
434 ordinal_type p = basis_orders[0];
435 ordinal_type p_i = std::min(total_order-sum_i, p);
436 ordinal_type p_j = std::min(total_order-sum_j, p);
437 ordinal_type p_k = std::min(total_order-sum_k, p);
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));
458 for (ordinal_type i=0; i<=p_i; ++i) {
459 for (ordinal_type
j=0;
j<=p_j; ++
j) {
465 ordinal_type k0 = parent_j_equals_k ?
j : 0;
466 if (symmetric && (k0%2 != (i+
j)%2)) ++k0;
467 const ordinal_type k_end = p_k;
468 const ordinal_type k_inc = symmetric ? 2 : 1;
469 for (ordinal_type k=k0; k<=k_end; k+=k_inc) {
470 value_type cijk = (*Cijk_1d[0])(i,
j, k);
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;
489 for (ordinal_type i=0; i<=p_i; ++i) {
490 for (ordinal_type
j=0;
j<=p_j; ++
j) {
491 ordinal_type k0 = parent_j_equals_k ? std::max(
j,std::abs(i-
j)) :
493 if (symmetric && (k0%2 != (i+
j)%2)) ++k0;
494 const ordinal_type k_end = std::min(p_k,i+
j);
495 const ordinal_type k_inc = symmetric ? 2 : 1;
496 for (ordinal_type k=k0; k<=k_end; k+=k_inc) {
497 value_type cijk = (*Cijk_1d[0])(i,
j, k);
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;
573 ordinal_type child_index;
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>
622 value_type ab[max_size][max_size];
625 c.
init(value_type(0));
628 const value_type *ca = a.
coeff();
629 const value_type *cb = b.
coeff();
630 value_type *cc = c.
coeff();
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) {
639 value_type *c_block = cc + ni->i_begin;
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;
644 const ordinal_type p_i = ni->p_i;
645 const ordinal_type p_j = ni->p_j;
646 const ordinal_type p_k = ni->p_k;
647 for (ordinal_type
j=0;
j<=p_j; ++
j)
648 for (ordinal_type k=0; k<=p_k; ++k)
649 ab[
j][k] = a_j_block[
j]*b_k_block[k] + a_k_block[k]*b_j_block[
j];
650 for (ordinal_type i=0; i<=p_i; ++i) {
651 value_type tmp = value_type(0);
652 for (ordinal_type
j=0;
j<=p_j; ++
j) {
660 ordinal_type k0 = ni->parent_j_equals_k ?
j : 0;
661 if (cijk.symmetric && (k0%2 != (i+
j)%2)) ++k0;
662 const ordinal_type k_end = p_k;
663 const ordinal_type k_inc = cijk.symmetric ? 2 : 1;
664 for (ordinal_type k=k0; k<=k_end; k+=k_inc) {
665 tmp += (*ci)*ab[
j][k];
Data structure storing a dense 3-tensor C(i,j,k).
Data structure storing a sparse 3-tensor C(i,j,k) in a a tree-based format for lexicographically orde...
LTBSparse3Tensor(const LTBSparse3Tensor &)
~LTBSparse3Tensor()
Destructor.
ordinal_type num_leafs() const
Return number of nodes.
bool symmetric() const
Return if symmetric.
Teuchos::RCP< CijkNode > head
value_type getValue(ordinal_type i, ordinal_type j, ordinal_type k) const
Get Cijk value for a given i, j, k indices.
void setHeadNode(const Teuchos::RCP< CijkNode > &h)
Set the head node.
void print(std::ostream &os) const
Print tensor.
LTBSparse3Tensor & operator=(const LTBSparse3Tensor &b)
LTBSparse3Tensor(const bool symm)
Constructor.
Teuchos::RCP< const CijkNode > getHeadNode() const
Get the head node.
ordinal_type num_entries() const
Return number of non-zero entries.
A comparison functor implementing a strict weak ordering based lexographic ordering.
A multidimensional index.
Class to store coefficients of a projection onto an orthogonal polynomial basis.
void init(const value_type &v)
Initialize coefficients to value.
pointer coeff()
Return coefficient array.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Top-level namespace for Stokhos classes and functions.
Teuchos::RCP< FlatLTBSparse3Tensor< ordinal_type, value_type > > computeFlatTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTBBlockLeaf(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
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)
std::ostream & operator<<(std::ostream &os, const ProductContainer< coeff_type > &vec)
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
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< 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)
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::Array< FlatLTBSparse3TensorNode< ordinal_type > > node_array_type
cijk_array_type::iterator cijk_iterator
cijk_array_type::const_iterator cijk_const_iterator
node_array_type::const_iterator node_const_iterator
Teuchos::Array< value_type > cijk_array_type
node_array_type::iterator node_iterator
Node type used in constructing the tree.
Teuchos::Array< Teuchos::RCP< CijkNode > > children
ordinal_type total_num_entries
ordinal_type total_num_leafs
ordinal_type my_num_entries
Teuchos::Array< value_type > values
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
LexicographicTreeBasisNode()
LexicographicTreeBasisNode & operator=(const LexicographicTreeBasisNode &node)
Teuchos::Array< Teuchos::RCP< LexicographicTreeBasisNode > > children
LexicographicTreeBasisNode(const LexicographicTreeBasisNode &node)