Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
48 #include "Stokhos_ProductBasis.hpp"
51 
52 #include "Kokkos_Core.hpp"
53 
54 //----------------------------------------------------------------------------
55 //----------------------------------------------------------------------------
56 
57 namespace Stokhos {
58 
67 template< typename ValueType, class ExecutionSpace, bool PackIndex >
69 
72 template< typename ValueType, class ExecutionSpace >
73 class CooProductTensor<ValueType,ExecutionSpace,true> {
74 public:
75 
76  typedef ExecutionSpace execution_space;
77  typedef typename execution_space::size_type size_type;
78  typedef ValueType value_type;
79 
80 private:
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 
114 public:
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
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 
281 template< typename ValueType, class ExecutionSpace>
282 class CooProductTensor<ValueType,ExecutionSpace,false> {
283 public:
284 
285  typedef ExecutionSpace execution_space;
286  typedef typename execution_space::size_type size_type;
287  typedef ValueType value_type;
288 
289 private:
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 
300 public:
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
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 
470 template< class Device, bool Pack, typename OrdinalType, typename ValueType >
471 CooProductTensor<ValueType, Device, Pack>
476 {
478  basis, Cijk, params );
479 }
480 
481 template <typename ValueType, typename Device, bool Pack>
482 std::ostream& operator << (
483  std::ostream& os, const CooProductTensor<ValueType,Device,Pack>& tensor)
484 {
485  tensor.print(os);
486  return os;
487 }
488 
489 template< typename ValueType, typename Device, bool Pack >
490 class BlockMultiply< CooProductTensor< ValueType, Device, Pack > >
491 {
492 public:
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 */
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::View< value_type[], execution_space > value_array_type
k_iterator k_begin() const
Iterator pointing to first k entry.
KOKKOS_INLINE_FUNCTION size_type dimension() const
Dimension of the tensor.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
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< size_type[][3], execution_space > coord_array_type
kj_iterator j_begin(const k_iterator &k) const
Iterator pointing to first j entry for given k.
KOKKOS_INLINE_FUNCTION size_type entry_count() const
Number of sparse entries.
static CooProductTensor create(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())
kj_iterator j_end(const k_iterator &k) const
Iterator pointing to last j entry for given k.
KOKKOS_INLINE_FUNCTION size_type num_non_zeros() const
Number of non-zero&#39;s.
KOKKOS_INLINE_FUNCTION size_type entry_count() const
Number of sparse entries.
std::ostream & operator<<(std::ostream &os, const ProductContainer< coeff_type > &vec)
Kokkos::View< size_type[], execution_space > coord_array_type
static size_type pack(const size_type i, const size_type j, const size_type k)
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
KOKKOS_INLINE_FUNCTION size_type num_non_zeros() const
Number of non-zero&#39;s.
k_iterator k_end() const
Iterator pointing to last k entry.
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)
Stokhos::Sparse3Tensor< int, double > Cijk_type
expr val()
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&#39;s per multiply-add.
KOKKOS_INLINE_FUNCTION size_type num_flops() const
Number flop&#39;s per multiply-add.
KOKKOS_INLINE_FUNCTION const value_type & value(const size_type entry) const
Value of an entry.
#define TEUCHOS_ASSERT(assertion_test)
static KOKKOS_INLINE_FUNCTION size_type vector_size(const tensor_type &tensor)
KOKKOS_INLINE_FUNCTION void unpack(size_type ijk, size_type &i, size_type &j, size_type &k) const
kji_iterator i_begin(const kj_iterator &j) const
Iterator pointing to first i entry for given j and k.
Sparse product tensor using &#39;COO&#39;-like storage format.
Kokkos::View< value_type[], execution_space > value_array_type
kji_iterator i_end(const kj_iterator &j) const
Iterator pointing to last i entry for given j and k.
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())
int n
virtual ordinal_type size() const =0
Return total size of basis.
KOKKOS_INLINE_FUNCTION const value_type & value(const size_type entry) const
Value of an entry.
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
static KOKKOS_INLINE_FUNCTION void apply(const tensor_type &tensor, const MatrixValue *const a, const VectorValue *const x, VectorValue *const y)