Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_LTBSparse3Tensor.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef STOKHOS_LTB_SPARSE_3_TENSOR_HPP
11 #define STOKHOS_LTB_SPARSE_3_TENSOR_HPP
12 
13 #include <ostream>
14 
15 #include "Stokhos_TotalOrderBasis.hpp"
16 #include "Stokhos_OrthogPolyApprox.hpp"
17 
18 namespace Stokhos {
19 
24  template <typename ordinal_type, typename value_type>
25  class LTBSparse3Tensor {
26  public:
27 
29  struct CijkNode {
30  Teuchos::Array< Teuchos::RCP<CijkNode> > children;
31  Teuchos::Array<value_type> values;
32  ordinal_type p_i, p_j, p_k;
33  ordinal_type i_begin, j_begin, k_begin;
34  ordinal_type i_size, j_size, k_size;
35  ordinal_type my_num_entries, total_num_entries;
36  ordinal_type total_num_leafs;
37  bool is_leaf, parent_j_equals_k;
38  };
39 
41  LTBSparse3Tensor(const bool symm) : is_symmetric(symm) {}
42 
44  ~LTBSparse3Tensor() {}
45 
47  void setHeadNode(const Teuchos::RCP<CijkNode>& h) { head = h; }
48 
50  Teuchos::RCP<const CijkNode> getHeadNode() const { return head; }
51 
53  void print(std::ostream& os) const {}
54 
56  value_type getValue(ordinal_type i, ordinal_type j, ordinal_type k) const
57  {}
58 
60  ordinal_type num_entries() const {
61  if (head != Teuchos::null)
62  return head->total_num_entries;
63  return 0;
64  }
65 
67  ordinal_type num_leafs() const {
68  if (head != Teuchos::null)
69  return head->total_num_leafs;
70  return 0;
71  }
72 
74  bool symmetric() const { return is_symmetric; }
75 
76  private:
77 
78  // Prohibit copying
79  LTBSparse3Tensor(const LTBSparse3Tensor&);
80 
81  // Prohibit Assignment
82  LTBSparse3Tensor& operator=(const LTBSparse3Tensor& b);
83 
84  protected:
85 
86  Teuchos::RCP<CijkNode> head;
87  bool is_symmetric;
88 
89  }; // class LTBSparse3Tensor
90 
94  template <typename ordinal_type, typename value_type>
95  std::ostream&
96  operator << (std::ostream& os,
97  const LTBSparse3Tensor<ordinal_type, value_type>& Cijk) {
98  Cijk.print(os);
99  return os;
100  }
101 
102  template <typename ordinal_type>
103  struct LexicographicTreeBasisNode {
104  Teuchos::Array< Teuchos::RCP<LexicographicTreeBasisNode> > children;
105  Teuchos::Array< Stokhos::MultiIndex<ordinal_type> > terms;
106  ordinal_type index_begin;
107  ordinal_type block_size;
108 
109  // Default constructor
110  LexicographicTreeBasisNode() :
111  children(), terms(), index_begin(0), block_size(0) {}
112 
113  // Copy constructor
114  LexicographicTreeBasisNode(const LexicographicTreeBasisNode& node) :
115  children(node.children.size()), terms(node.terms),
116  index_begin(node.index_begin), block_size(node.block_size) {
117  for (ordinal_type i=0; i<children.size(); ++i)
118  children[i] =
119  Teuchos::rcp(new LexicographicTreeBasisNode(*(node->children[i])));
120  }
121 
122  // Assignment operator
123  LexicographicTreeBasisNode&
124  operator=(const LexicographicTreeBasisNode& node) {
125  if (this != &node) {
126  children.resize(node.children.size());
127  for (ordinal_type i=0; i<children.size(); ++i)
128  children[i] =
129  Teuchos::rcp(new LexicographicTreeBasisNode(*(node->children[i])));
130  terms = node.terms;
131  index_begin = node.index_begin;
132  block_size = node.block_size;
133  }
134  return *this;
135  }
136  };
137 
138  template <typename ordinal_type>
139  Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> >
140  build_lexicographic_basis_tree(
141  const Teuchos::ArrayView<const ordinal_type>& basis_orders,
142  const ordinal_type total_order,
143  const ordinal_type index_begin = ordinal_type(0),
144  const ordinal_type order_sum = ordinal_type(0),
145  const Stokhos::MultiIndex<ordinal_type>& term_prefix =
146  Stokhos::MultiIndex<ordinal_type>()) {
147 
148  typedef LexicographicTreeBasisNode<ordinal_type> node_type;
149 
150  ordinal_type my_d = basis_orders.size();
151  ordinal_type prev_d = term_prefix.dimension();
152  ordinal_type p = basis_orders[0];
153  ordinal_type my_p = std::min(total_order-order_sum, p);
154 
155  Teuchos::RCP<node_type> node = Teuchos::rcp(new node_type);
156  node->index_begin = index_begin;
157  node->terms.resize(my_p+1);
158  for (ordinal_type i=0; i<=my_p; ++i) {
159  node->terms[i].resize(prev_d+1);
160  for (ordinal_type j=0; j<prev_d; ++j)
161  node->terms[i][j] = term_prefix[j];
162  node->terms[i][prev_d] = i;
163  }
164 
165  // Base case for dimension = 1
166  if (my_d == 1) {
167  node->block_size = my_p+1;
168  }
169 
170  // General case -- call recursively stripping off first dimension
171  else {
172  Teuchos::ArrayView<const ordinal_type> bo =
173  Teuchos::arrayView(&basis_orders[1],my_d-1);
174  node->children.resize(my_p+1);
175  node->index_begin = index_begin;
176  node->block_size = 0;
177  for (ordinal_type i=0; i<=my_p; ++i) {
178  Teuchos::RCP<node_type> child = build_lexicographic_basis_tree(
179  bo, total_order, index_begin+node->block_size, order_sum+i,
180  node->terms[i]);
181  node->children[i] = child;
182  node->block_size += child->block_size;
183  }
184  }
185 
186  return node;
187  }
188 
189  // An approach for building a sparse 3-tensor only for lexicographically
190  // ordered total order basis using a tree-based format
191  template <typename ordinal_type,
192  typename value_type>
193  Teuchos::RCP< LTBSparse3Tensor<ordinal_type, value_type> >
194  computeTripleProductTensorLTB(
195  const TotalOrderBasis<ordinal_type, value_type,LexographicLess<MultiIndex<ordinal_type> > >& product_basis,
196  bool symmetric = false) {
197 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
198  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Time");
199 #endif
200  using Teuchos::RCP;
201  using Teuchos::rcp;
202  using Teuchos::Array;
203  using Teuchos::ArrayView;
204 
205  const Array< RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases = product_basis.getCoordinateBases();
206  ordinal_type d = bases.size();
207  ordinal_type p = product_basis.order();
208  Array<ordinal_type> basis_orders(d);
209  for (int i=0; i<d; ++i)
210  basis_orders[i] = bases[i]->order();
211  ArrayView<const ordinal_type> basis_orders_view = basis_orders();
212 
213  // Create 1-D triple products
214  Array< RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
215  for (ordinal_type i=0; i<d; i++) {
216  Cijk_1d[i] =
217  bases[i]->computeSparseTripleProductTensor(bases[i]->order()+1);
218  }
219  ArrayView<const RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d_view
220  = Cijk_1d();
221 
222  // Create lexicographic tree basis
223  Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> > ltb =
224  build_lexicographic_basis_tree(basis_orders_view, p);
225 
226  // Current implementation is recursive on the dimension d
227  typedef LTBSparse3Tensor<ordinal_type, value_type> Cijk_type;
228  RCP<Cijk_type> Cijk = rcp(new Cijk_type(symmetric));
229  RCP<typename Cijk_type::CijkNode> head =
230  computeCijkLTBNode(
231  basis_orders_view, Cijk_1d_view, ltb, ltb, ltb, p, symmetric);
232  Cijk->setHeadNode(head);
233 
234  return Cijk;
235  }
236 
237  template <typename ordinal_type,
238  typename value_type>
239  Teuchos::RCP<typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode>
240  computeCijkLTBNode(
241  const Teuchos::ArrayView<const ordinal_type>& basis_orders,
242  const Teuchos::ArrayView<const Teuchos::RCP<Sparse3Tensor<ordinal_type,value_type> > >& Cijk_1d,
243  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& i_ltb,
244  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& j_ltb,
245  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& k_ltb,
246  const ordinal_type total_order,
247  const bool symmetric,
248  const ordinal_type sum_i = ordinal_type(0),
249  const ordinal_type sum_j = ordinal_type(0),
250  const ordinal_type sum_k = ordinal_type(0),
251  const value_type cijk_base = value_type(1)) {
252 
253  using Teuchos::RCP;
254  using Teuchos::rcp;
255  using Teuchos::ArrayView;
256  using Teuchos::arrayView;
257  typedef typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode node_type;
258  typedef ProductBasisUtils::Cijk_1D_Iterator<ordinal_type> Cijk_Iterator;
259 
260  ordinal_type my_d = basis_orders.size();
261  ordinal_type p = basis_orders[0];
262  ordinal_type p_i = std::min(total_order-sum_i, p);
263  ordinal_type p_j = std::min(total_order-sum_j, p);
264  ordinal_type p_k = std::min(total_order-sum_k, p);
265  Cijk_Iterator cijk_iterator(p_i, p_j, p_k, symmetric);
266 
267  RCP<node_type> node = rcp(new node_type);
268  node->p_i = p_i;
269  node->p_j = p_j;
270  node->p_k = p_k;
271  node->i_begin = i_ltb->index_begin;
272  node->j_begin = j_ltb->index_begin;
273  node->k_begin = k_ltb->index_begin;
274  node->i_size = i_ltb->block_size;
275  node->j_size = j_ltb->block_size;
276  node->k_size = k_ltb->block_size;
277 
278  // Base case -- compute the actual cijk values
279  if (my_d == 1) {
280  node->is_leaf = true;
281  bool more = true;
282  while (more) {
283  value_type cijk =
284  Cijk_1d[0]->getValue(cijk_iterator.i,
285  cijk_iterator.j,
286  cijk_iterator.k);
287  node->values.push_back(cijk*cijk_base);
288  more = cijk_iterator.increment();
289  }
290  node->my_num_entries = node->values.size();
291  node->total_num_entries = node->values.size();
292  node->total_num_leafs = 1;
293  }
294 
295  // General case -- call recursively stripping off first dimension
296  else {
297  node->is_leaf = false;
298  ArrayView<const ordinal_type> bo = arrayView(&basis_orders[1], my_d-1);
299  ArrayView<const RCP<Sparse3Tensor<ordinal_type,value_type> > > c1d =
300  arrayView(&Cijk_1d[1], my_d-1);
301  node->total_num_entries = 0;
302  node->total_num_leafs = 0;
303  bool more = true;
304  while (more) {
305  value_type cijk =
306  Cijk_1d[0]->getValue(cijk_iterator.i,
307  cijk_iterator.j,
308  cijk_iterator.k);
309  RCP<node_type> child =
310  computeCijkLTBNode(bo, c1d,
311  i_ltb->children[cijk_iterator.i],
312  j_ltb->children[cijk_iterator.j],
313  k_ltb->children[cijk_iterator.k],
314  total_order, symmetric,
315  sum_i + cijk_iterator.i,
316  sum_j + cijk_iterator.j,
317  sum_k + cijk_iterator.k,
318  cijk_base*cijk);
319  node->children.push_back(child);
320  node->total_num_entries += child->total_num_entries;
321  node->total_num_leafs += child->total_num_leafs;
322  more = cijk_iterator.increment();
323  }
324  node->my_num_entries = node->children.size();
325  }
326 
327  return node;
328  }
329 
330  // An approach for building a sparse 3-tensor only for lexicographically
331  // ordered total order basis using a tree-based format -- the leaf nodes
332  // are replaced by a dense p_i x p_j x p_k block
333  template <typename ordinal_type,
334  typename value_type>
335  Teuchos::RCP< LTBSparse3Tensor<ordinal_type, value_type> >
336  computeTripleProductTensorLTBBlockLeaf(
337  const TotalOrderBasis<ordinal_type, value_type,LexographicLess<MultiIndex<ordinal_type> > >& product_basis,
338  bool symmetric = false) {
339 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
340  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Time");
341 #endif
342  using Teuchos::RCP;
343  using Teuchos::rcp;
344  using Teuchos::Array;
345  using Teuchos::ArrayView;
346 
347  const Array< RCP<const OneDOrthogPolyBasis<ordinal_type, value_type> > >& bases = product_basis.getCoordinateBases();
348  ordinal_type d = bases.size();
349  ordinal_type p = product_basis.order();
350  Array<ordinal_type> basis_orders(d);
351  for (int i=0; i<d; ++i)
352  basis_orders[i] = bases[i]->order();
353  ArrayView<const ordinal_type> basis_orders_view = basis_orders();
354 
355  // Create 1-D triple products
356  Array< RCP<Dense3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
357  for (ordinal_type i=0; i<d; i++) {
358  Cijk_1d[i] = bases[i]->computeTripleProductTensor();
359  }
360  ArrayView<const RCP<Dense3Tensor<ordinal_type,value_type> > > Cijk_1d_view
361  = Cijk_1d();
362 
363  // Create lexicographic tree basis
364  Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> > ltb =
365  build_lexicographic_basis_tree(basis_orders_view, p);
366 
367  // Current implementation is recursive on the dimension d
368  typedef LTBSparse3Tensor<ordinal_type, value_type> Cijk_type;
369  RCP<Cijk_type> Cijk = rcp(new Cijk_type(symmetric));
370  RCP<typename Cijk_type::CijkNode> head =
371  computeCijkLTBNodeBlockLeaf(
372  basis_orders_view, Cijk_1d_view, ltb, ltb, ltb, p, symmetric);
373  Cijk->setHeadNode(head);
374 
375  return Cijk;
376  }
377 
378  template <typename ordinal_type,
379  typename value_type>
380  Teuchos::RCP<typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode>
381  computeCijkLTBNodeBlockLeaf(
382  const Teuchos::ArrayView<const ordinal_type>& basis_orders,
383  const Teuchos::ArrayView<const Teuchos::RCP<Dense3Tensor<ordinal_type,value_type> > >& Cijk_1d,
384  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& i_ltb,
385  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& j_ltb,
386  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& k_ltb,
387  const ordinal_type total_order,
388  const bool symmetric,
389  const ordinal_type sum_i = ordinal_type(0),
390  const ordinal_type sum_j = ordinal_type(0),
391  const ordinal_type sum_k = ordinal_type(0),
392  const value_type cijk_base = value_type(1),
393  const bool parent_j_equals_k = true) {
394 
395  using Teuchos::RCP;
396  using Teuchos::rcp;
397  using Teuchos::ArrayView;
398  using Teuchos::arrayView;
399  typedef typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode node_type;
400 
401  ordinal_type my_d = basis_orders.size();
402  ordinal_type p = basis_orders[0];
403  ordinal_type p_i = std::min(total_order-sum_i, p);
404  ordinal_type p_j = std::min(total_order-sum_j, p);
405  ordinal_type p_k = std::min(total_order-sum_k, p);
406 
407  RCP<node_type> node = rcp(new node_type);
408  node->p_i = p_i;
409  node->p_j = p_j;
410  node->p_k = p_k;
411  node->i_begin = i_ltb->index_begin;
412  node->j_begin = j_ltb->index_begin;
413  node->k_begin = k_ltb->index_begin;
414  node->i_size = i_ltb->block_size;
415  node->j_size = j_ltb->block_size;
416  node->k_size = k_ltb->block_size;
417  node->parent_j_equals_k = parent_j_equals_k;
418 
419  // Base case -- compute the actual cijk values as a "brick"
420  // Could store values as a "pyramid" using commented out code below
421  // However that seems to result in very bad performance, e.g., a lot
422  // of register spill in multiply code based on this tensor
423  if (my_d == 1) {
424  node->is_leaf = true;
425  node->values.reserve((p_i+1)*(p_j+1)*(p_k+1));
426  for (ordinal_type i=0; i<=p_i; ++i) {
427  for (ordinal_type j=0; j<=p_j; ++j) {
428  // ordinal_type k0 = parent_j_equals_k ? std::max(j,std::abs(i-j)) :
429  // std::abs(i-j);
430  // if (symmetric && (k0%2 != (i+j)%2)) ++k0;
431  // const ordinal_type k_end = std::min(p_k,i+j);
432  // const ordinal_type k_inc = symmetric ? 2 : 1;
433  ordinal_type k0 = parent_j_equals_k ? j : 0;
434  if (symmetric && (k0%2 != (i+j)%2)) ++k0;
435  const ordinal_type k_end = p_k;
436  const ordinal_type k_inc = symmetric ? 2 : 1;
437  for (ordinal_type k=k0; k<=k_end; k+=k_inc) {
438  value_type cijk = (*Cijk_1d[0])(i, j, k);
439  if (j+node->j_begin == k+node->k_begin) cijk *= 0.5;
440  node->values.push_back(cijk*cijk_base);
441  }
442  }
443  }
444  node->my_num_entries = node->values.size();
445  node->total_num_entries = node->values.size();
446  node->total_num_leafs = 1;
447  }
448 
449  // General case -- call recursively stripping off first dimension
450  else {
451  node->is_leaf = false;
452  ArrayView<const ordinal_type> bo = arrayView(&basis_orders[1], my_d-1);
453  ArrayView<const RCP<Dense3Tensor<ordinal_type,value_type> > > c1d =
454  arrayView(&Cijk_1d[1], my_d-1);
455  node->total_num_entries = 0;
456  node->total_num_leafs = 0;
457  for (ordinal_type i=0; i<=p_i; ++i) {
458  for (ordinal_type j=0; j<=p_j; ++j) {
459  ordinal_type k0 = parent_j_equals_k ? std::max(j,std::abs(i-j)) :
460  std::abs(i-j);
461  if (symmetric && (k0%2 != (i+j)%2)) ++k0;
462  const ordinal_type k_end = std::min(p_k,i+j);
463  const ordinal_type k_inc = symmetric ? 2 : 1;
464  for (ordinal_type k=k0; k<=k_end; k+=k_inc) {
465  value_type cijk = (*Cijk_1d[0])(i, j, k);
466  RCP<node_type> child =
467  computeCijkLTBNodeBlockLeaf(bo, c1d,
468  i_ltb->children[i],
469  j_ltb->children[j],
470  k_ltb->children[k],
471  total_order, symmetric,
472  sum_i + i,
473  sum_j + j,
474  sum_k + k,
475  cijk_base*cijk,
476  parent_j_equals_k && j == k);
477  node->children.push_back(child);
478  node->total_num_entries += child->total_num_entries;
479  node->total_num_leafs += child->total_num_leafs;
480  }
481  }
482  }
483  node->my_num_entries = node->children.size();
484  }
485 
486  return node;
487  }
488 
489  template <typename ordinal_type>
490  struct FlatLTBSparse3TensorNode {
491  ordinal_type i_begin, j_begin, k_begin;
492  ordinal_type p_i, p_j, p_k;
493  bool parent_j_equals_k;
494  };
495 
496  template <typename ordinal_type, typename value_type>
497  struct FlatLTBSparse3Tensor {
498  typedef Teuchos::Array< FlatLTBSparse3TensorNode<ordinal_type> > node_array_type;
499  typedef Teuchos::Array< value_type > cijk_array_type;
500 
501  typedef typename node_array_type::iterator node_iterator;
502  typedef typename node_array_type::const_iterator node_const_iterator;
503  typedef typename cijk_array_type::iterator cijk_iterator;
504  typedef typename cijk_array_type::const_iterator cijk_const_iterator;
505 
506  node_array_type node;
507  cijk_array_type cijk;
508  bool symmetric;
509  };
510 
511  template <typename ordinal_type, typename value_type>
512  Teuchos::RCP< FlatLTBSparse3Tensor<ordinal_type,value_type> >
513  computeFlatTripleProductTensorLTB(
514  const TotalOrderBasis<ordinal_type, value_type,LexographicLess<MultiIndex<ordinal_type> > >& product_basis,
515  bool symmetric = false) {
516 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
517  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Flat Triple-Product Tensor Time");
518 #endif
519  using Teuchos::RCP;
520  using Teuchos::rcp;
521 
522  // Build LTB 3 tensor
523  typedef LTBSparse3Tensor<ordinal_type, value_type> Cijk_LTB_type;
524  RCP<Cijk_LTB_type> ltb_tensor =
525  computeTripleProductTensorLTBBlockLeaf(product_basis, symmetric);
526 
527  // Create flat LTB 3 tensor
528  RCP< FlatLTBSparse3Tensor<ordinal_type,value_type> > flat_tensor =
529  rcp(new FlatLTBSparse3Tensor<ordinal_type,value_type>);
530  flat_tensor->node.reserve(ltb_tensor->num_leafs());
531  flat_tensor->cijk.reserve(ltb_tensor->num_entries());
532  flat_tensor->symmetric = symmetric;
533 
534  // Fill flat 3 tensor
535  typedef typename Cijk_LTB_type::CijkNode node_type;
536  Teuchos::Array< Teuchos::RCP<const node_type> > node_stack;
537  Teuchos::Array< ordinal_type > index_stack;
538  node_stack.push_back(ltb_tensor->getHeadNode());
539  index_stack.push_back(0);
540  Teuchos::RCP<const node_type> node;
541  ordinal_type child_index;
542  while (node_stack.size() > 0) {
543  node = node_stack.back();
544  child_index = index_stack.back();
545 
546  // Leaf
547  if (node->is_leaf) {
548  FlatLTBSparse3TensorNode<ordinal_type> leaf;
549  leaf.i_begin = node->i_begin;
550  leaf.j_begin = node->j_begin;
551  leaf.k_begin = node->k_begin;
552  leaf.p_i = node->p_i;
553  leaf.p_j = node->p_j;
554  leaf.p_k = node->p_k;
555  leaf.parent_j_equals_k = node->parent_j_equals_k;
556  flat_tensor->node.push_back(leaf);
557  flat_tensor->cijk.insert(flat_tensor->cijk.end(),
558  node->values.begin(),
559  node->values.end());
560  node_stack.pop_back();
561  index_stack.pop_back();
562  }
563 
564  // More children to process -- process them first
565  else if (child_index < node->children.size()) {
566  ++index_stack.back();
567  node = node->children[child_index];
568  node_stack.push_back(node);
569  index_stack.push_back(0);
570  }
571 
572  // No more children
573  else {
574  node_stack.pop_back();
575  index_stack.pop_back();
576  }
577 
578  }
579 
580  return flat_tensor;
581  }
582 
583  template <int max_size, typename ordinal_type, typename value_type>
584  void
585  flatLTB3TensorMultiply(
586  OrthogPolyApprox<ordinal_type,value_type>& c,
587  const OrthogPolyApprox<ordinal_type,value_type>& a,
588  const OrthogPolyApprox<ordinal_type,value_type>& b,
589  const FlatLTBSparse3Tensor<ordinal_type,value_type>& cijk) {
590  value_type ab[max_size][max_size];
591 
592  // Set coefficients to 0
593  c.init(value_type(0));
594 
595  // Get pointers to coefficients
596  const value_type *ca = a.coeff();
597  const value_type *cb = b.coeff();
598  value_type *cc = c.coeff();
599 
600  typedef FlatLTBSparse3Tensor<ordinal_type,value_type> cijk_type;
601  typedef typename cijk_type::node_const_iterator node_iterator;
602  typedef typename cijk_type::cijk_const_iterator cijk_iterator;
603  node_iterator ni = cijk.node.begin();
604  node_iterator ni_end = cijk.node.end();
605  cijk_iterator ci = cijk.cijk.begin();
606  for (; ni != ni_end; ++ni) {
607  value_type *c_block = cc + ni->i_begin;
608  const value_type *a_j_block = ca + ni->j_begin;
609  const value_type *b_k_block = cb + ni->k_begin;
610  const value_type *a_k_block = ca + ni->k_begin;
611  const value_type *b_j_block = cb + ni->j_begin;
612  const ordinal_type p_i = ni->p_i;
613  const ordinal_type p_j = ni->p_j;
614  const ordinal_type p_k = ni->p_k;
615  for (ordinal_type j=0; j<=p_j; ++j)
616  for (ordinal_type k=0; k<=p_k; ++k)
617  ab[j][k] = a_j_block[j]*b_k_block[k] + a_k_block[k]*b_j_block[j];
618  for (ordinal_type i=0; i<=p_i; ++i) {
619  value_type tmp = value_type(0);
620  for (ordinal_type j=0; j<=p_j; ++j) {
621  // This is for pyramid instead of brick
622  // ordinal_type k0 = ni->parent_j_equals_k ? std::max(j,std::abs(i-j)) :
623  // std::abs(i-j);
624  // if (cijk.symmetric && (k0%2 != (i+j)%2)) ++k0;
625  // const ordinal_type k_end = std::min(p_k,i+j);
626  // const ordinal_type k_inc = cijk.symmetric ? 2 : 1;
627 
628  ordinal_type k0 = ni->parent_j_equals_k ? j : 0;
629  if (cijk.symmetric && (k0%2 != (i+j)%2)) ++k0;
630  const ordinal_type k_end = p_k;
631  const ordinal_type k_inc = cijk.symmetric ? 2 : 1;
632  for (ordinal_type k=k0; k<=k_end; k+=k_inc) {
633  tmp += (*ci)*ab[j][k];
634  ++ci;
635  }
636  }
637  c_block[i] += tmp;
638  }
639  }
640  }
641 
642 } // namespace Stokhos
643 
644 // Include template definitions
645 //#include "Stokhos_LTBSparse3TensorImp.hpp"
646 
647 #endif // STOKHOS_LTB_SPARSE_3_TENSOR_HPP
Stokhos::LTBSparse3Tensor::getHeadNode
Teuchos::RCP< const CijkNode > getHeadNode() const
Get the head node.
Definition: Stokhos_LTBSparse3Tensor.hpp:50
TEUCHOS_FUNC_TIME_MONITOR
#define TEUCHOS_FUNC_TIME_MONITOR(FUNCNAME)
Stokhos::computeTripleProductTensorLTB
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Definition: Stokhos_LTBSparse3Tensor.hpp:194
Stokhos::TotalOrderBasis
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Definition: Stokhos_TotalOrderBasis.hpp:36
Stokhos::Sparse3Tensor
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Definition: Stokhos_Sparse3Tensor.hpp:24
Stokhos::ProductBasisUtils::Cijk_1D_Iterator
Definition: Stokhos_ProductBasisUtils.hpp:1261
Stokhos::OrthogPolyApprox::init
void init(const value_type &v)
Initialize coefficients to value.
Definition: Stokhos_OrthogPolyApproxImp.hpp:72
Teuchos::ArrayView::size
size_type size() const
Stokhos::LTBSparse3Tensor::LTBSparse3Tensor
LTBSparse3Tensor(const bool symm)
Constructor.
Definition: Stokhos_LTBSparse3Tensor.hpp:41
Stokhos::FlatLTBSparse3Tensor::cijk_iterator
cijk_array_type::iterator cijk_iterator
Definition: Stokhos_LTBSparse3Tensor.hpp:503
Stokhos::FlatLTBSparse3Tensor::node_array_type
Teuchos::Array< FlatLTBSparse3TensorNode< ordinal_type > > node_array_type
Definition: Stokhos_LTBSparse3Tensor.hpp:498
Stokhos::LexicographicTreeBasisNode::LexicographicTreeBasisNode
LexicographicTreeBasisNode(const LexicographicTreeBasisNode &node)
Definition: Stokhos_LTBSparse3Tensor.hpp:114
Stokhos::LexicographicTreeBasisNode::terms
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
Definition: Stokhos_LTBSparse3Tensor.hpp:105
Stokhos::OrthogPolyApprox::coeff
pointer coeff()
Return coefficient array.
Definition: Stokhos_OrthogPolyApproxImp.hpp:134
Stokhos::LTBSparse3Tensor::num_entries
ordinal_type num_entries() const
Return number of non-zero entries.
Definition: Stokhos_LTBSparse3Tensor.hpp:60
Stokhos::computeTripleProductTensorLTBBlockLeaf
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTBBlockLeaf(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Definition: Stokhos_LTBSparse3Tensor.hpp:336
Stokhos::LTBSparse3Tensor::num_leafs
ordinal_type num_leafs() const
Return number of nodes.
Definition: Stokhos_LTBSparse3Tensor.hpp:67
Sacado::UQ::min
KOKKOS_INLINE_FUNCTION PCE< Storage > min(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:1217
node_type
Kokkos::Serial node_type
Definition: Stokhos_SacadoTraitsUnitTest.cpp:29
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:495
Stokhos::LexicographicTreeBasisNode::operator=
LexicographicTreeBasisNode & operator=(const LexicographicTreeBasisNode &node)
Definition: Stokhos_LTBSparse3Tensor.hpp:124
Stokhos::FlatLTBSparse3Tensor::node_const_iterator
node_array_type::const_iterator node_const_iterator
Definition: Stokhos_LTBSparse3Tensor.hpp:502
Stokhos::FlatLTBSparse3Tensor::cijk_const_iterator
cijk_array_type::const_iterator cijk_const_iterator
Definition: Stokhos_LTBSparse3Tensor.hpp:504
Stokhos::LTBSparse3Tensor
Data structure storing a sparse 3-tensor C(i,j,k) in a a tree-based format for lexicographically orde...
Definition: Stokhos_LTBSparse3Tensor.hpp:25
Stokhos::FlatLTBSparse3Tensor::node_iterator
node_array_type::iterator node_iterator
Definition: Stokhos_LTBSparse3Tensor.hpp:501
Stokhos::build_lexicographic_basis_tree
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 >())
Definition: Stokhos_LTBSparse3Tensor.hpp:140
Stokhos::computeFlatTripleProductTensorLTB
Teuchos::RCP< FlatLTBSparse3Tensor< ordinal_type, value_type > > computeFlatTripleProductTensorLTB(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Definition: Stokhos_LTBSparse3Tensor.hpp:513
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Stokhos::operator<<
std::ostream & operator<<(std::ostream &os, const ProductContainer< coeff_type > &vec)
Definition: Stokhos_ProductContainer.hpp:152
Stokhos::LexicographicTreeBasisNode::children
Teuchos::Array< Teuchos::RCP< LexicographicTreeBasisNode > > children
Definition: Stokhos_LTBSparse3Tensor.hpp:104
Stokhos::flatLTB3TensorMultiply
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)
Definition: Stokhos_LTBSparse3Tensor.hpp:585
Sacado::UQ::max
KOKKOS_INLINE_FUNCTION PCE< Storage > max(const typename PCE< Storage >::value_type &a, const PCE< Storage > &b)
Definition: Sacado_UQ_PCE_Imp.hpp:1180
Stokhos_OrthogPolyApprox.hpp
Stokhos::computeCijkLTBNode
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))
Definition: Stokhos_LTBSparse3Tensor.hpp:240
Stokhos::computeCijkLTBNodeBlockLeaf
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)
Definition: Stokhos_LTBSparse3Tensor.hpp:381
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1150
Teuchos::Array::end
iterator end()
Teuchos::Array< FlatLTBSparse3TensorNode< ordinal_type > >::const_iterator
std::vector< FlatLTBSparse3TensorNode< ordinal_type > >::const_iterator const_iterator
Kokkos::cijk
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)
Definition: KokkosExp_View_UQ_PCE_Contiguous.hpp:88
Sparse3TensorUnitTest::Cijk_type
Stokhos::Sparse3Tensor< int, double > Cijk_type
Definition: Stokhos_Sparse3TensorUnitTest.cpp:20
Teuchos::Array::back
reference back()
Teuchos::Array::push_back
void push_back(const value_type &x)
Stokhos::LTBSparse3Tensor::symmetric
bool symmetric() const
Return if symmetric.
Definition: Stokhos_LTBSparse3Tensor.hpp:74
Stokhos::LTBSparse3Tensor::CijkNode
Node type used in constructing the tree.
Definition: Stokhos_LTBSparse3Tensor.hpp:29
Teuchos::Array::size
size_type size() const
Stokhos::LTBSparse3Tensor::print
void print(std::ostream &os) const
Print tensor.
Definition: Stokhos_LTBSparse3Tensor.hpp:53
Stokhos::LexographicLess
A comparison functor implementing a strict weak ordering based lexographic ordering.
Definition: Stokhos_ProductBasisUtils.hpp:767
Stokhos::Dense3Tensor
Data structure storing a dense 3-tensor C(i,j,k).
Definition: Stokhos_Dense3Tensor.hpp:23
Stokhos::LTBSparse3Tensor::setHeadNode
void setHeadNode(const Teuchos::RCP< CijkNode > &h)
Set the head node.
Definition: Stokhos_LTBSparse3Tensor.hpp:47
Stokhos::FlatLTBSparse3Tensor::cijk_array_type
Teuchos::Array< value_type > cijk_array_type
Definition: Stokhos_LTBSparse3Tensor.hpp:499
Stokhos::LTBSparse3Tensor::CijkNode::children
Teuchos::Array< Teuchos::RCP< CijkNode > > children
Definition: Stokhos_LTBSparse3Tensor.hpp:30
Teuchos::Array< FlatLTBSparse3TensorNode< ordinal_type > >::iterator
std::vector< FlatLTBSparse3TensorNode< ordinal_type > >::iterator iterator
Teuchos::Array::begin
iterator begin()
Stokhos_TotalOrderBasis.hpp
Stokhos::LTBSparse3Tensor::getValue
value_type getValue(ordinal_type i, ordinal_type j, ordinal_type k) const
Get Cijk value for a given i, j, k indices.
Definition: Stokhos_LTBSparse3Tensor.hpp:56
Stokhos::LTBSparse3Tensor::operator=
LTBSparse3Tensor & operator=(const LTBSparse3Tensor &b)
Generated on Fri Nov 22 2024 09:18:13 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5