Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
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
49
50namespace Stokhos {
51
56 template <typename ordinal_type, typename value_type>
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;
68 ordinal_type total_num_leafs;
70 };
71
73 LTBSparse3Tensor(const bool symm) : is_symmetric(symm) {}
74
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
112
113 // Prohibit Assignment
115
116 protected:
117
118 Teuchos::RCP<CijkNode> head;
120
121 }; // class LTBSparse3Tensor
122
126 template <typename ordinal_type, typename value_type>
127 std::ostream&
128 operator << (std::ostream& os,
130 Cijk.print(os);
131 return os;
132 }
133
134 template <typename ordinal_type>
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
143 children(), terms(), index_begin(0), block_size(0) {}
144
145 // Copy constructor
147 children(node.children.size()), terms(node.terms),
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
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;
164 block_size = node.block_size;
165 }
166 return *this;
167 }
168 };
169
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),
177 const Stokhos::MultiIndex<ordinal_type>& term_prefix =
179
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> >
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
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);
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>
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;
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> >
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
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);
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>
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;
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 =
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>
523 ordinal_type i_begin, j_begin, k_begin;
524 ordinal_type p_i, p_j, p_k;
526 };
527
528 template <typename ordinal_type, typename value_type>
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
541 };
542
543 template <typename ordinal_type, typename value_type>
544 Teuchos::RCP< FlatLTBSparse3Tensor<ordinal_type,value_type> >
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 =
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) {
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
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
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
Kokkos::Serial node_type
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 &)
ordinal_type num_leafs() const
Return number of nodes.
bool symmetric() const
Return if symmetric.
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
Teuchos::Array< Stokhos::MultiIndex< ordinal_type > > terms
LexicographicTreeBasisNode & operator=(const LexicographicTreeBasisNode &node)
Teuchos::Array< Teuchos::RCP< LexicographicTreeBasisNode > > children
LexicographicTreeBasisNode(const LexicographicTreeBasisNode &node)