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 //
4 // Stokhos Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Eric T. Phipps (etphipp@sandia.gov).
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #ifndef STOKHOS_LTB_SPARSE_3_TENSOR_HPP
43 #define STOKHOS_LTB_SPARSE_3_TENSOR_HPP
44 
45 #include <ostream>
46 
47 #include "Stokhos_TotalOrderBasis.hpp"
48 #include "Stokhos_OrthogPolyApprox.hpp"
49 
50 namespace Stokhos {
51 
56  template <typename ordinal_type, typename value_type>
57  class LTBSparse3Tensor {
58  public:
59 
61  struct CijkNode {
62  Teuchos::Array< Teuchos::RCP<CijkNode> > children;
63  Teuchos::Array<value_type> values;
64  ordinal_type p_i, p_j, p_k;
65  ordinal_type i_begin, j_begin, k_begin;
66  ordinal_type i_size, j_size, k_size;
67  ordinal_type my_num_entries, total_num_entries;
68  ordinal_type total_num_leafs;
69  bool is_leaf, parent_j_equals_k;
70  };
71 
73  LTBSparse3Tensor(const bool symm) : is_symmetric(symm) {}
74 
76  ~LTBSparse3Tensor() {}
77 
79  void setHeadNode(const Teuchos::RCP<CijkNode>& h) { head = h; }
80 
82  Teuchos::RCP<const CijkNode> getHeadNode() const { return head; }
83 
85  void print(std::ostream& os) const {}
86 
88  value_type getValue(ordinal_type i, ordinal_type j, ordinal_type k) const
89  {}
90 
92  ordinal_type num_entries() const {
93  if (head != Teuchos::null)
94  return head->total_num_entries;
95  return 0;
96  }
97 
99  ordinal_type num_leafs() const {
100  if (head != Teuchos::null)
101  return head->total_num_leafs;
102  return 0;
103  }
104 
106  bool symmetric() const { return is_symmetric; }
107 
108  private:
109 
110  // Prohibit copying
111  LTBSparse3Tensor(const LTBSparse3Tensor&);
112 
113  // Prohibit Assignment
114  LTBSparse3Tensor& operator=(const LTBSparse3Tensor& b);
115 
116  protected:
117 
118  Teuchos::RCP<CijkNode> head;
119  bool is_symmetric;
120 
121  }; // class LTBSparse3Tensor
122 
126  template <typename ordinal_type, typename value_type>
127  std::ostream&
128  operator << (std::ostream& os,
129  const LTBSparse3Tensor<ordinal_type, value_type>& Cijk) {
130  Cijk.print(os);
131  return os;
132  }
133 
134  template <typename ordinal_type>
135  struct LexicographicTreeBasisNode {
136  Teuchos::Array< Teuchos::RCP<LexicographicTreeBasisNode> > children;
137  Teuchos::Array< Stokhos::MultiIndex<ordinal_type> > terms;
138  ordinal_type index_begin;
139  ordinal_type block_size;
140 
141  // Default constructor
142  LexicographicTreeBasisNode() :
143  children(), terms(), index_begin(0), block_size(0) {}
144 
145  // Copy constructor
146  LexicographicTreeBasisNode(const LexicographicTreeBasisNode& node) :
147  children(node.children.size()), terms(node.terms),
148  index_begin(node.index_begin), block_size(node.block_size) {
149  for (ordinal_type i=0; i<children.size(); ++i)
150  children[i] =
151  Teuchos::rcp(new LexicographicTreeBasisNode(*(node->children[i])));
152  }
153 
154  // Assignment operator
155  LexicographicTreeBasisNode&
156  operator=(const LexicographicTreeBasisNode& node) {
157  if (this != &node) {
158  children.resize(node.children.size());
159  for (ordinal_type i=0; i<children.size(); ++i)
160  children[i] =
161  Teuchos::rcp(new LexicographicTreeBasisNode(*(node->children[i])));
162  terms = node.terms;
163  index_begin = node.index_begin;
164  block_size = node.block_size;
165  }
166  return *this;
167  }
168  };
169 
170  template <typename ordinal_type>
171  Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> >
172  build_lexicographic_basis_tree(
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),
177  const Stokhos::MultiIndex<ordinal_type>& term_prefix =
178  Stokhos::MultiIndex<ordinal_type>()) {
179 
180  typedef LexicographicTreeBasisNode<ordinal_type> node_type;
181 
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);
186 
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;
195  }
196 
197  // Base case for dimension = 1
198  if (my_d == 1) {
199  node->block_size = my_p+1;
200  }
201 
202  // General case -- call recursively stripping off first dimension
203  else {
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) {
210  Teuchos::RCP<node_type> child = build_lexicographic_basis_tree(
211  bo, total_order, index_begin+node->block_size, order_sum+i,
212  node->terms[i]);
213  node->children[i] = child;
214  node->block_size += child->block_size;
215  }
216  }
217 
218  return node;
219  }
220 
221  // An approach for building a sparse 3-tensor only for lexicographically
222  // ordered total order basis using a tree-based format
223  template <typename ordinal_type,
224  typename value_type>
225  Teuchos::RCP< LTBSparse3Tensor<ordinal_type, value_type> >
226  computeTripleProductTensorLTB(
227  const TotalOrderBasis<ordinal_type, value_type,LexographicLess<MultiIndex<ordinal_type> > >& product_basis,
228  bool symmetric = false) {
229 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
230  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Time");
231 #endif
232  using Teuchos::RCP;
233  using Teuchos::rcp;
234  using Teuchos::Array;
235  using Teuchos::ArrayView;
236 
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();
244 
245  // Create 1-D triple products
246  Array< RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d(d);
247  for (ordinal_type i=0; i<d; i++) {
248  Cijk_1d[i] =
249  bases[i]->computeSparseTripleProductTensor(bases[i]->order()+1);
250  }
251  ArrayView<const RCP<Sparse3Tensor<ordinal_type,value_type> > > Cijk_1d_view
252  = Cijk_1d();
253 
254  // Create lexicographic tree basis
255  Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> > ltb =
256  build_lexicographic_basis_tree(basis_orders_view, p);
257 
258  // Current implementation is recursive on the dimension d
259  typedef LTBSparse3Tensor<ordinal_type, value_type> Cijk_type;
260  RCP<Cijk_type> Cijk = rcp(new Cijk_type(symmetric));
261  RCP<typename Cijk_type::CijkNode> head =
262  computeCijkLTBNode(
263  basis_orders_view, Cijk_1d_view, ltb, ltb, ltb, p, symmetric);
264  Cijk->setHeadNode(head);
265 
266  return Cijk;
267  }
268 
269  template <typename ordinal_type,
270  typename value_type>
271  Teuchos::RCP<typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode>
272  computeCijkLTBNode(
273  const Teuchos::ArrayView<const ordinal_type>& basis_orders,
274  const Teuchos::ArrayView<const Teuchos::RCP<Sparse3Tensor<ordinal_type,value_type> > >& Cijk_1d,
275  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& i_ltb,
276  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& j_ltb,
277  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& k_ltb,
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)) {
284 
285  using Teuchos::RCP;
286  using Teuchos::rcp;
287  using Teuchos::ArrayView;
288  using Teuchos::arrayView;
289  typedef typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode node_type;
290  typedef ProductBasisUtils::Cijk_1D_Iterator<ordinal_type> Cijk_Iterator;
291 
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);
298 
299  RCP<node_type> node = rcp(new node_type);
300  node->p_i = p_i;
301  node->p_j = p_j;
302  node->p_k = p_k;
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;
309 
310  // Base case -- compute the actual cijk values
311  if (my_d == 1) {
312  node->is_leaf = true;
313  bool more = true;
314  while (more) {
315  value_type cijk =
316  Cijk_1d[0]->getValue(cijk_iterator.i,
317  cijk_iterator.j,
318  cijk_iterator.k);
319  node->values.push_back(cijk*cijk_base);
320  more = cijk_iterator.increment();
321  }
322  node->my_num_entries = node->values.size();
323  node->total_num_entries = node->values.size();
324  node->total_num_leafs = 1;
325  }
326 
327  // General case -- call recursively stripping off first dimension
328  else {
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;
335  bool more = true;
336  while (more) {
337  value_type cijk =
338  Cijk_1d[0]->getValue(cijk_iterator.i,
339  cijk_iterator.j,
340  cijk_iterator.k);
341  RCP<node_type> child =
342  computeCijkLTBNode(bo, c1d,
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,
350  cijk_base*cijk);
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();
355  }
356  node->my_num_entries = node->children.size();
357  }
358 
359  return node;
360  }
361 
362  // An approach for building a sparse 3-tensor only for lexicographically
363  // ordered total order basis using a tree-based format -- the leaf nodes
364  // are replaced by a dense p_i x p_j x p_k block
365  template <typename ordinal_type,
366  typename value_type>
367  Teuchos::RCP< LTBSparse3Tensor<ordinal_type, value_type> >
368  computeTripleProductTensorLTBBlockLeaf(
369  const TotalOrderBasis<ordinal_type, value_type,LexographicLess<MultiIndex<ordinal_type> > >& product_basis,
370  bool symmetric = false) {
371 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
372  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Total Triple-Product Tensor Time");
373 #endif
374  using Teuchos::RCP;
375  using Teuchos::rcp;
376  using Teuchos::Array;
377  using Teuchos::ArrayView;
378 
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();
386 
387  // Create 1-D triple products
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();
391  }
392  ArrayView<const RCP<Dense3Tensor<ordinal_type,value_type> > > Cijk_1d_view
393  = Cijk_1d();
394 
395  // Create lexicographic tree basis
396  Teuchos::RCP< LexicographicTreeBasisNode<ordinal_type> > ltb =
397  build_lexicographic_basis_tree(basis_orders_view, p);
398 
399  // Current implementation is recursive on the dimension d
400  typedef LTBSparse3Tensor<ordinal_type, value_type> Cijk_type;
401  RCP<Cijk_type> Cijk = rcp(new Cijk_type(symmetric));
402  RCP<typename Cijk_type::CijkNode> head =
403  computeCijkLTBNodeBlockLeaf(
404  basis_orders_view, Cijk_1d_view, ltb, ltb, ltb, p, symmetric);
405  Cijk->setHeadNode(head);
406 
407  return Cijk;
408  }
409 
410  template <typename ordinal_type,
411  typename value_type>
412  Teuchos::RCP<typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode>
413  computeCijkLTBNodeBlockLeaf(
414  const Teuchos::ArrayView<const ordinal_type>& basis_orders,
415  const Teuchos::ArrayView<const Teuchos::RCP<Dense3Tensor<ordinal_type,value_type> > >& Cijk_1d,
416  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& i_ltb,
417  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& j_ltb,
418  const Teuchos::RCP<LexicographicTreeBasisNode<ordinal_type> >& k_ltb,
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) {
426 
427  using Teuchos::RCP;
428  using Teuchos::rcp;
429  using Teuchos::ArrayView;
430  using Teuchos::arrayView;
431  typedef typename LTBSparse3Tensor<ordinal_type, value_type>::CijkNode node_type;
432 
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);
438 
439  RCP<node_type> node = rcp(new node_type);
440  node->p_i = p_i;
441  node->p_j = p_j;
442  node->p_k = p_k;
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;
450 
451  // Base case -- compute the actual cijk values as a "brick"
452  // Could store values as a "pyramid" using commented out code below
453  // However that seems to result in very bad performance, e.g., a lot
454  // of register spill in multiply code based on this tensor
455  if (my_d == 1) {
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) {
460  // ordinal_type k0 = parent_j_equals_k ? std::max(j,std::abs(i-j)) :
461  // std::abs(i-j);
462  // if (symmetric && (k0%2 != (i+j)%2)) ++k0;
463  // const ordinal_type k_end = std::min(p_k,i+j);
464  // const ordinal_type k_inc = symmetric ? 2 : 1;
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);
473  }
474  }
475  }
476  node->my_num_entries = node->values.size();
477  node->total_num_entries = node->values.size();
478  node->total_num_leafs = 1;
479  }
480 
481  // General case -- call recursively stripping off first dimension
482  else {
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)) :
492  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 =
499  computeCijkLTBNodeBlockLeaf(bo, c1d,
500  i_ltb->children[i],
501  j_ltb->children[j],
502  k_ltb->children[k],
503  total_order, symmetric,
504  sum_i + i,
505  sum_j + j,
506  sum_k + k,
507  cijk_base*cijk,
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;
512  }
513  }
514  }
515  node->my_num_entries = node->children.size();
516  }
517 
518  return node;
519  }
520 
521  template <typename ordinal_type>
522  struct FlatLTBSparse3TensorNode {
523  ordinal_type i_begin, j_begin, k_begin;
524  ordinal_type p_i, p_j, p_k;
525  bool parent_j_equals_k;
526  };
527 
528  template <typename ordinal_type, typename value_type>
529  struct FlatLTBSparse3Tensor {
530  typedef Teuchos::Array< FlatLTBSparse3TensorNode<ordinal_type> > node_array_type;
531  typedef Teuchos::Array< value_type > cijk_array_type;
532 
533  typedef typename node_array_type::iterator node_iterator;
534  typedef typename node_array_type::const_iterator node_const_iterator;
535  typedef typename cijk_array_type::iterator cijk_iterator;
536  typedef typename cijk_array_type::const_iterator cijk_const_iterator;
537 
538  node_array_type node;
539  cijk_array_type cijk;
540  bool symmetric;
541  };
542 
543  template <typename ordinal_type, typename value_type>
544  Teuchos::RCP< FlatLTBSparse3Tensor<ordinal_type,value_type> >
545  computeFlatTripleProductTensorLTB(
546  const TotalOrderBasis<ordinal_type, value_type,LexographicLess<MultiIndex<ordinal_type> > >& product_basis,
547  bool symmetric = false) {
548 #ifdef STOKHOS_TEUCHOS_TIME_MONITOR
549  TEUCHOS_FUNC_TIME_MONITOR("Stokhos: Flat Triple-Product Tensor Time");
550 #endif
551  using Teuchos::RCP;
552  using Teuchos::rcp;
553 
554  // Build LTB 3 tensor
555  typedef LTBSparse3Tensor<ordinal_type, value_type> Cijk_LTB_type;
556  RCP<Cijk_LTB_type> ltb_tensor =
557  computeTripleProductTensorLTBBlockLeaf(product_basis, symmetric);
558 
559  // Create flat LTB 3 tensor
560  RCP< FlatLTBSparse3Tensor<ordinal_type,value_type> > flat_tensor =
561  rcp(new FlatLTBSparse3Tensor<ordinal_type,value_type>);
562  flat_tensor->node.reserve(ltb_tensor->num_leafs());
563  flat_tensor->cijk.reserve(ltb_tensor->num_entries());
564  flat_tensor->symmetric = symmetric;
565 
566  // Fill flat 3 tensor
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();
577 
578  // Leaf
579  if (node->is_leaf) {
580  FlatLTBSparse3TensorNode<ordinal_type> leaf;
581  leaf.i_begin = node->i_begin;
582  leaf.j_begin = node->j_begin;
583  leaf.k_begin = node->k_begin;
584  leaf.p_i = node->p_i;
585  leaf.p_j = node->p_j;
586  leaf.p_k = node->p_k;
587  leaf.parent_j_equals_k = node->parent_j_equals_k;
588  flat_tensor->node.push_back(leaf);
589  flat_tensor->cijk.insert(flat_tensor->cijk.end(),
590  node->values.begin(),
591  node->values.end());
592  node_stack.pop_back();
593  index_stack.pop_back();
594  }
595 
596  // More children to process -- process them first
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);
602  }
603 
604  // No more children
605  else {
606  node_stack.pop_back();
607  index_stack.pop_back();
608  }
609 
610  }
611 
612  return flat_tensor;
613  }
614 
615  template <int max_size, typename ordinal_type, typename value_type>
616  void
617  flatLTB3TensorMultiply(
618  OrthogPolyApprox<ordinal_type,value_type>& c,
619  const OrthogPolyApprox<ordinal_type,value_type>& a,
620  const OrthogPolyApprox<ordinal_type,value_type>& b,
621  const FlatLTBSparse3Tensor<ordinal_type,value_type>& cijk) {
622  value_type ab[max_size][max_size];
623 
624  // Set coefficients to 0
625  c.init(value_type(0));
626 
627  // Get pointers to coefficients
628  const value_type *ca = a.coeff();
629  const value_type *cb = b.coeff();
630  value_type *cc = c.coeff();
631 
632  typedef FlatLTBSparse3Tensor<ordinal_type,value_type> cijk_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) {
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) {
653  // This is for pyramid instead of brick
654  // ordinal_type k0 = ni->parent_j_equals_k ? std::max(j,std::abs(i-j)) :
655  // std::abs(i-j);
656  // if (cijk.symmetric && (k0%2 != (i+j)%2)) ++k0;
657  // const ordinal_type k_end = std::min(p_k,i+j);
658  // const ordinal_type k_inc = cijk.symmetric ? 2 : 1;
659 
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];
666  ++ci;
667  }
668  }
669  c_block[i] += tmp;
670  }
671  }
672  }
673 
674 } // namespace Stokhos
675 
676 // Include template definitions
677 //#include "Stokhos_LTBSparse3TensorImp.hpp"
678 
679 #endif // STOKHOS_LTB_SPARSE_3_TENSOR_HPP
Stokhos::LTBSparse3Tensor::getHeadNode
Teuchos::RCP< const CijkNode > getHeadNode() const
Get the head node.
Definition: Stokhos_LTBSparse3Tensor.hpp:82
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:226
Stokhos::TotalOrderBasis
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Definition: Stokhos_TotalOrderBasis.hpp:68
Stokhos::Sparse3Tensor
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Definition: Stokhos_Sparse3Tensor.hpp:56
Stokhos::ProductBasisUtils::Cijk_1D_Iterator
Definition: Stokhos_ProductBasisUtils.hpp:1293
Stokhos::OrthogPolyApprox::init
void init(const value_type &v)
Initialize coefficients to value.
Definition: Stokhos_OrthogPolyApproxImp.hpp:106
Teuchos::ArrayView::size
size_type size() const
Stokhos::LTBSparse3Tensor::LTBSparse3Tensor
LTBSparse3Tensor(const bool symm)
Constructor.
Definition: Stokhos_LTBSparse3Tensor.hpp:73
Stokhos::FlatLTBSparse3Tensor::cijk_iterator
cijk_array_type::iterator cijk_iterator
Definition: Stokhos_LTBSparse3Tensor.hpp:535
Stokhos::FlatLTBSparse3Tensor::node_array_type
Teuchos::Array< FlatLTBSparse3TensorNode< ordinal_type > > node_array_type
Definition: Stokhos_LTBSparse3Tensor.hpp:530
Stokhos::LexicographicTreeBasisNode::LexicographicTreeBasisNode
LexicographicTreeBasisNode(const LexicographicTreeBasisNode &node)
Definition: Stokhos_LTBSparse3Tensor.hpp:146
Stokhos::LexicographicTreeBasisNode::terms
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
Definition: Stokhos_LTBSparse3Tensor.hpp:137
Stokhos::OrthogPolyApprox::coeff
pointer coeff()
Return coefficient array.
Definition: Stokhos_OrthogPolyApproxImp.hpp:168
Stokhos::LTBSparse3Tensor::num_entries
ordinal_type num_entries() const
Return number of non-zero entries.
Definition: Stokhos_LTBSparse3Tensor.hpp:92
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:368
Stokhos::LTBSparse3Tensor::num_leafs
ordinal_type num_leafs() const
Return number of nodes.
Definition: Stokhos_LTBSparse3Tensor.hpp:99
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:1247
node_type
Kokkos::Serial node_type
Definition: Stokhos_SacadoTraitsUnitTest.cpp:59
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Stokhos::LexicographicTreeBasisNode::operator=
LexicographicTreeBasisNode & operator=(const LexicographicTreeBasisNode &node)
Definition: Stokhos_LTBSparse3Tensor.hpp:156
Stokhos::FlatLTBSparse3Tensor::node_const_iterator
node_array_type::const_iterator node_const_iterator
Definition: Stokhos_LTBSparse3Tensor.hpp:534
Stokhos::FlatLTBSparse3Tensor::cijk_const_iterator
cijk_array_type::const_iterator cijk_const_iterator
Definition: Stokhos_LTBSparse3Tensor.hpp:536
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:57
Stokhos::FlatLTBSparse3Tensor::node_iterator
node_array_type::iterator node_iterator
Definition: Stokhos_LTBSparse3Tensor.hpp:533
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:172
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:545
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:184
Stokhos::LexicographicTreeBasisNode::children
Teuchos::Array< Teuchos::RCP< LexicographicTreeBasisNode > > children
Definition: Stokhos_LTBSparse3Tensor.hpp:136
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:617
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:1211
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:272
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:413
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1181
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:217
Sparse3TensorUnitTest::Cijk_type
Stokhos::Sparse3Tensor< int, double > Cijk_type
Definition: Stokhos_Sparse3TensorUnitTest.cpp:52
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:106
Stokhos::LTBSparse3Tensor::CijkNode
Node type used in constructing the tree.
Definition: Stokhos_LTBSparse3Tensor.hpp:61
Teuchos::Array::size
size_type size() const
Stokhos::LTBSparse3Tensor::print
void print(std::ostream &os) const
Print tensor.
Definition: Stokhos_LTBSparse3Tensor.hpp:85
Stokhos::LexographicLess
A comparison functor implementing a strict weak ordering based lexographic ordering.
Definition: Stokhos_ProductBasisUtils.hpp:799
Stokhos::Dense3Tensor
Data structure storing a dense 3-tensor C(i,j,k).
Definition: Stokhos_Dense3Tensor.hpp:57
Stokhos::LTBSparse3Tensor::setHeadNode
void setHeadNode(const Teuchos::RCP< CijkNode > &h)
Set the head node.
Definition: Stokhos_LTBSparse3Tensor.hpp:79
Stokhos::FlatLTBSparse3Tensor::cijk_array_type
Teuchos::Array< value_type > cijk_array_type
Definition: Stokhos_LTBSparse3Tensor.hpp:531
Stokhos::LTBSparse3Tensor::CijkNode::children
Teuchos::Array< Teuchos::RCP< CijkNode > > children
Definition: Stokhos_LTBSparse3Tensor.hpp:62
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:88
Stokhos::LTBSparse3Tensor::operator=
LTBSparse3Tensor & operator=(const LTBSparse3Tensor &b)
Generated by   doxygen 1.8.5