Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_CooProductTensor.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_COO_PRODUCT_TENSOR_HPP
43#define STOKHOS_COO_PRODUCT_TENSOR_HPP
44
45#include "Kokkos_Core.hpp"
46
47#include "Stokhos_Multiply.hpp"
50#include "Teuchos_ParameterList.hpp"
51
52#include "Kokkos_Core.hpp"
53
54//----------------------------------------------------------------------------
55//----------------------------------------------------------------------------
56
57namespace Stokhos {
58
67template< typename ValueType, class ExecutionSpace, bool PackIndex >
69
72template< typename ValueType, class ExecutionSpace >
73class CooProductTensor<ValueType,ExecutionSpace,true> {
74public:
75
76 typedef ExecutionSpace execution_space;
77 typedef typename execution_space::size_type size_type;
78 typedef ValueType value_type;
79
80private:
81
82 // Number of bits available for each index
83 static const size_type bits = (sizeof(size_type)*8) / 3;
84
85 // Mask for packing index
86 static const size_type mask = (1 << bits)-1;
87
88 typedef Kokkos::View< value_type[], execution_space > vec_type;
89 typedef Kokkos::View< size_type[], execution_space > coord_array_type;
90 typedef Kokkos::View< value_type[], execution_space > value_array_type;
91
96
97 // Pack (i,j,k) into a single integer
98 static size_type
99 pack( const size_type i, const size_type j, const size_type k ) {
100 const size_type ii = i & mask;
101 const size_type jj = j & mask;
102 const size_type kk = k & mask;
103 size_type ijk = ii | (jj << bits) | (kk << 2*bits);
104 return ijk;
105 }
106
107 KOKKOS_INLINE_FUNCTION
108 void unpack( size_type ijk, size_type& i, size_type& j, size_type& k ) const {
109 i = ijk & mask; ijk >>= bits;
110 j = ijk & mask;
111 k = ijk >> bits;
112 }
113
114public:
115
117 static const size_type max_index = 1 << bits;
118
119 inline
121
122 inline
124 m_coord(),
125 m_value(),
126 m_dim(0),
127 m_flops(0) {}
128
129 inline
131 m_coord( rhs.m_coord ),
132 m_value( rhs.m_value ),
133 m_dim( rhs.m_dim ),
134 m_flops( rhs.m_flops ) {}
135
136 inline
137 CooProductTensor & operator = ( const CooProductTensor & rhs )
138 {
139 m_coord = rhs.m_coord;
140 m_value = rhs.m_value;
141 m_dim = rhs.m_dim;
142 m_flops = rhs.m_flops;
143 return *this;
144 }
145
147 KOKKOS_INLINE_FUNCTION
148 size_type dimension() const { return m_dim; }
149
151 KOKKOS_INLINE_FUNCTION
152 size_type entry_count() const { return m_coord.extent(0); }
153
155 KOKKOS_INLINE_FUNCTION
156 void coord( const size_type entry,
157 size_type& i, size_type& j, size_type& k ) const {
158 unpack(m_coord(entry), i, j, k);
159 }
160
162 KOKKOS_INLINE_FUNCTION
163 const value_type & value( const size_type entry ) const
164 { return m_value(entry); }
165
167 KOKKOS_INLINE_FUNCTION
168 size_type num_non_zeros() const { return m_coord.extent(0); }
169
171 KOKKOS_INLINE_FUNCTION
172 size_type num_flops() const { return m_flops; }
173
174 template <typename OrdinalType>
175 static CooProductTensor
178 const Teuchos::ParameterList& params = Teuchos::ParameterList())
179 {
181 typedef typename Cijk_type::i_iterator i_iterator;
182 typedef typename Cijk_type::ik_iterator k_iterator;
183 typedef typename Cijk_type::ikj_iterator j_iterator;
184
185 // Compute entry count
186 size_type entry_count = 0;
187 for (i_iterator i_it = Cijk.i_begin(); i_it!=Cijk.i_end(); ++i_it) {
188 for (k_iterator k_it = Cijk.k_begin(i_it); k_it != Cijk.k_end(i_it);
189 ++k_it) {
190 OrdinalType k = index(k_it);
191 for (j_iterator j_it = Cijk.j_begin(k_it); j_it != Cijk.j_end(k_it);
192 ++j_it) {
193 OrdinalType j = index(j_it);
194 if (j >= k) {
195 ++entry_count;
196 }
197 }
198 }
199 }
200
201 // Align entry_count
202#if defined( KOKKOS_ENABLE_CUDA )
203 enum { Align = std::is_same<ExecutionSpace,Kokkos::Cuda>::value ? 32 : 1 };
204#else
205 enum { Align = 1 };
206#endif
207
208 entry_count = (entry_count+Align-1) & ~(Align-1);
209 TEUCHOS_ASSERT(entry_count % Align == 0);
210
211 // Allocate tensor data
212 CooProductTensor tensor;
213 tensor.m_coord = coord_array_type( "tensor_coord", entry_count );
214 tensor.m_value = value_array_type( "tensor_value", entry_count );
215 tensor.m_dim = basis.size();
216 tensor.m_flops = 5*entry_count + tensor.m_dim;
217
218 // Create mirror, is a view if is host memory
219 typename coord_array_type::HostMirror
220 host_coord = Kokkos::create_mirror_view( tensor.m_coord );
221 typename value_array_type::HostMirror
222 host_value = Kokkos::create_mirror_view( tensor.m_value );
223
224 size_type n = 0;
225 OrdinalType i=0, j=0, k=0;
226 for (i_iterator i_it = Cijk.i_begin(); i_it!=Cijk.i_end(); ++i_it) {
227 i = index(i_it);
228 for (k_iterator k_it = Cijk.k_begin(i_it); k_it != Cijk.k_end(i_it);
229 ++k_it) {
230 k = index(k_it);
231 for (j_iterator j_it = Cijk.j_begin(k_it); j_it != Cijk.j_end(k_it);
232 ++j_it) {
233 j = index(j_it);
234 ValueType c = Stokhos::value(j_it);
235 if (j >= k) {
236 host_value(n) = (j != k) ? c : 0.5*c;
237 host_coord(n) = pack(i,j,k);
238 ++n;
239 }
240 }
241 }
242 }
243 for (; n < entry_count; ++n) {
244 host_value(n) = 0.0;
245 host_coord(n) = pack(i,j,k);
246 }
247
248 // Copy data to device if necessary
249 Kokkos::deep_copy( tensor.m_coord, host_coord );
250 Kokkos::deep_copy( tensor.m_value, host_value );
251
252 return tensor;
253 }
254
255 void print(std::ostream& os) const {
256 size_type num_entry = entry_count();
257 typename coord_array_type::HostMirror
258 host_coord = Kokkos::create_mirror_view( m_coord );
259 typename value_array_type::HostMirror
260 host_value = Kokkos::create_mirror_view( m_value );
261 Kokkos::deep_copy( host_coord, m_coord );
262 Kokkos::deep_copy( host_value, m_value );
263
264 os << "CooProductTensor: dim = "
265 << dimension() << ", entry_count = "
266 << num_entry << std::endl
267 << "Entries: i j k : cijk" << std::endl;
268 for (size_type l=0; l<num_entry; ++l) {
269 size_type i,j,k;
270 unpack(host_coord(l), i, j, k);
271 ValueType cijk = host_value(l);
272 os << "\t " << l << " : " << i << " " << j << " " << k << " = " << cijk
273 << std::endl;
274 if (l > 0 && l % 32 == 0) os << std::endl;
275 }
276 }
277};
278
281template< typename ValueType, class ExecutionSpace>
282class CooProductTensor<ValueType,ExecutionSpace,false> {
283public:
284
285 typedef ExecutionSpace execution_space;
286 typedef typename execution_space::size_type size_type;
287 typedef ValueType value_type;
288
289private:
290
291 typedef Kokkos::View< value_type[], execution_space > vec_type;
292 typedef Kokkos::View< size_type[][3], execution_space > coord_array_type;
293 typedef Kokkos::View< value_type[], execution_space > value_array_type;
294
299
300public:
301
302 inline
304
305 inline
307 m_coord(),
308 m_value(),
309 m_dim(0),
310 m_flops(0) {}
311
312 inline
314 m_coord( rhs.m_coord ),
315 m_value( rhs.m_value ),
316 m_dim( rhs.m_dim ),
317 m_flops( rhs.m_flops ) {}
318
319 inline
320 CooProductTensor & operator = ( const CooProductTensor & rhs )
321 {
322 m_coord = rhs.m_coord;
323 m_value = rhs.m_value;
324 m_dim = rhs.m_dim;
325 m_flops = rhs.m_flops;
326 return *this;
327 }
328
330 KOKKOS_INLINE_FUNCTION
331 size_type dimension() const { return m_dim; }
332
334 KOKKOS_INLINE_FUNCTION
335 size_type entry_count() const { return m_coord.extent(0); }
336
338 KOKKOS_INLINE_FUNCTION
339 void coord( const size_type entry,
340 size_type& i, size_type& j, size_type& k ) const {
341 i = m_coord(entry,0);
342 j = m_coord(entry,1);
343 k = m_coord(entry,2);
344 }
345
347 KOKKOS_INLINE_FUNCTION
348 const value_type & value( const size_type entry ) const
349 { return m_value( entry ); }
350
352 KOKKOS_INLINE_FUNCTION
353 size_type num_non_zeros() const { return m_coord.extent(0); }
354
356 KOKKOS_INLINE_FUNCTION
357 size_type num_flops() const { return m_flops; }
358
359 template <typename OrdinalType>
360 static CooProductTensor
363 const Teuchos::ParameterList& params = Teuchos::ParameterList())
364 {
366 typedef typename Cijk_type::i_iterator i_iterator;
367 typedef typename Cijk_type::ik_iterator k_iterator;
368 typedef typename Cijk_type::ikj_iterator j_iterator;
369
370 // Compute entry count
371 size_type entry_count = 0;
372 for (i_iterator i_it = Cijk.i_begin(); i_it!=Cijk.i_end(); ++i_it) {
373 for (k_iterator k_it = Cijk.k_begin(i_it); k_it != Cijk.k_end(i_it);
374 ++k_it) {
375 OrdinalType k = index(k_it);
376 for (j_iterator j_it = Cijk.j_begin(k_it); j_it != Cijk.j_end(k_it);
377 ++j_it) {
378 OrdinalType j = index(j_it);
379 if (j >= k) {
380 ++entry_count;
381 }
382 }
383 }
384 }
385
386 // Align entry_count
387#if defined( KOKKOS_ENABLE_CUDA )
388 enum { Align = std::is_same<ExecutionSpace,Kokkos::Cuda>::value ? 32 : 1 };
389#else
390 enum { Align = 1 };
391#endif
392
393 entry_count = (entry_count+Align-1) & ~(Align-1);
394 TEUCHOS_ASSERT(entry_count % Align == 0);
395
396 // Allocate tensor data
397 CooProductTensor tensor;
398 tensor.m_coord = coord_array_type( "tensor_coord", entry_count );
399 tensor.m_value = value_array_type( "tensor_value", entry_count );
400 tensor.m_dim = basis.size();
401 tensor.m_flops = 5*entry_count + tensor.m_dim;
402
403 // Create mirror, is a view if is host memory
404 typename coord_array_type::HostMirror
405 host_coord = Kokkos::create_mirror_view( tensor.m_coord );
406 typename value_array_type::HostMirror
407 host_value = Kokkos::create_mirror_view( tensor.m_value );
408
409 // Set entries
410 size_type n = 0;
411 OrdinalType i=0, j=0, k=0;
412 for (i_iterator i_it = Cijk.i_begin(); i_it!=Cijk.i_end(); ++i_it) {
413 i = index(i_it);
414 for (k_iterator k_it = Cijk.k_begin(i_it); k_it != Cijk.k_end(i_it);
415 ++k_it) {
416 k = index(k_it);
417 for (j_iterator j_it = Cijk.j_begin(k_it); j_it != Cijk.j_end(k_it);
418 ++j_it) {
419 j = index(j_it);
420 ValueType c = Stokhos::value(j_it);
421 if (j >= k) {
422 host_value(n) = (j != k) ? c : 0.5*c;
423 host_coord(n,0) = i;
424 host_coord(n,1) = j;
425 host_coord(n,2) = k;
426 ++n;
427 }
428 }
429 }
430 }
431 for (; n < entry_count; ++n) {
432 host_value(n) = 0.0;
433 host_coord(n,0) = i;
434 host_coord(n,1) = j;
435 host_coord(n,2) = k;
436 }
437
438 // Copy data to device if necessary
439 Kokkos::deep_copy( tensor.m_coord, host_coord );
440 Kokkos::deep_copy( tensor.m_value, host_value );
441
442 return tensor;
443 }
444
445 void print(std::ostream& os) const {
446 size_type num_entry = entry_count();
447 typename coord_array_type::HostMirror
448 host_coord = Kokkos::create_mirror_view( m_coord );
449 typename value_array_type::HostMirror
450 host_value = Kokkos::create_mirror_view( m_value );
451 Kokkos::deep_copy( host_coord, m_coord );
452 Kokkos::deep_copy( host_value, m_value );
453
454 os << "CooProductTensor: dim = "
455 << dimension() << ", entry_count = "
456 << num_entry << std::endl
457 << "Entries: i j k : cijk" << std::endl;
458 for (size_type l=0; l<num_entry; ++l) {
459 size_type i = host_coord(l,0);
460 size_type j = host_coord(l,1);
461 size_type k = host_coord(l,2);
462 ValueType cijk = host_value(l);
463 os << "\t " << l << " : " << i << " " << j << " " << k << " = " << cijk
464 << std::endl;
465 if (l > 0 && l % 32 == 0) os << std::endl;
466 }
467 }
468};
469
470template< class Device, bool Pack, typename OrdinalType, typename ValueType >
471CooProductTensor<ValueType, Device, Pack>
475 const Teuchos::ParameterList& params = Teuchos::ParameterList())
476{
478 basis, Cijk, params );
479}
480
481template <typename ValueType, typename Device, bool Pack>
482std::ostream& operator << (
483 std::ostream& os, const CooProductTensor<ValueType,Device,Pack>& tensor)
484{
485 tensor.print(os);
486 return os;
487}
488
489template< typename ValueType, typename Device, bool Pack >
490class BlockMultiply< CooProductTensor< ValueType, Device, Pack > >
491{
492public:
493
494 typedef typename Device::size_type size_type;
496
497 template< typename MatrixValue , typename VectorValue >
498 KOKKOS_INLINE_FUNCTION
499 static void apply( const tensor_type & tensor,
500 const MatrixValue * const a,
501 const VectorValue * const x,
502 VectorValue * const y )
503 {
504 const size_type nEntry = tensor.entry_count();
505 size_type i = 0, j = 0, k = 0, i_prev = -1;
506 VectorValue val = 0.0, carry_val = 0.0;
507 for ( size_type entry = 0 ; entry < nEntry ; ++entry ) {
508 tensor.coord(entry, i, j, k);
509 val = tensor.value(entry) * ( a[j] * x[k] + a[k] * x[j] );
510 if (i == i_prev)
511 carry_val += val;
512 else {
513 y[i_prev] += carry_val;
514 carry_val = val;
515 }
516 i_prev = i;
517 }
518 y[i] += carry_val;
519 }
520
521 KOKKOS_INLINE_FUNCTION
522 static size_type matrix_size( const tensor_type & tensor )
523 { return tensor.dimension(); }
524
525 KOKKOS_INLINE_FUNCTION
526 static size_type vector_size( const tensor_type & tensor )
527 { return tensor.dimension(); }
528};
529
530} /* namespace Stokhos */
531
532//----------------------------------------------------------------------------
533//----------------------------------------------------------------------------
534
535#endif /* #ifndef STOKHOS_COO_PRODUCT_TENSOR_HPP */
expr val()
static KOKKOS_INLINE_FUNCTION void apply(const tensor_type &tensor, const MatrixValue *const a, const VectorValue *const x, VectorValue *const y)
static KOKKOS_INLINE_FUNCTION size_type vector_size(const tensor_type &tensor)
static KOKKOS_INLINE_FUNCTION size_type matrix_size(const tensor_type &tensor)
KOKKOS_INLINE_FUNCTION size_type dimension() const
Dimension of the tensor.
Kokkos::View< value_type[], execution_space > value_array_type
KOKKOS_INLINE_FUNCTION size_type entry_count() const
Number of sparse entries.
KOKKOS_INLINE_FUNCTION size_type num_flops() const
Number flop's per multiply-add.
KOKKOS_INLINE_FUNCTION const value_type & value(const size_type entry) const
Value of an entry.
static CooProductTensor create(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())
KOKKOS_INLINE_FUNCTION size_type num_non_zeros() const
Number of non-zero's.
Kokkos::View< size_type[][3], execution_space > coord_array_type
KOKKOS_INLINE_FUNCTION void coord(const size_type entry, size_type &i, size_type &j, size_type &k) const
Get (i,j,k) coordinates of an entry.
static CooProductTensor create(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())
KOKKOS_INLINE_FUNCTION void coord(const size_type entry, size_type &i, size_type &j, size_type &k) const
Get (i,j,k) coordinates of an entry.
KOKKOS_INLINE_FUNCTION size_type num_flops() const
Number flop's per multiply-add.
KOKKOS_INLINE_FUNCTION size_type entry_count() const
Number of sparse entries.
KOKKOS_INLINE_FUNCTION size_type num_non_zeros() const
Number of non-zero's.
Kokkos::View< size_type[], execution_space > coord_array_type
Kokkos::View< value_type[], execution_space > value_array_type
KOKKOS_INLINE_FUNCTION void unpack(size_type ijk, size_type &i, size_type &j, size_type &k) const
static size_type pack(const size_type i, const size_type j, const size_type k)
KOKKOS_INLINE_FUNCTION size_type dimension() const
Dimension of the tensor.
KOKKOS_INLINE_FUNCTION const value_type & value(const size_type entry) const
Value of an entry.
Sparse product tensor using 'COO'-like storage format.
virtual ordinal_type size() const =0
Return total size of basis.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
k_iterator k_begin() const
Iterator pointing to first k entry.
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
k_iterator k_end() const
Iterator pointing to last k entry.
Top-level namespace for Stokhos classes and functions.
CooProductTensor< ValueType, Device, Pack > create_coo_product_tensor(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())
std::ostream & operator<<(std::ostream &os, const ProductContainer< coeff_type > &vec)