Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_FlatSparse3Tensor.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_FLAT_SPARSE_3_TENSOR_HPP
43#define STOKHOS_FLAT_SPARSE_3_TENSOR_HPP
44
45#include "Kokkos_Core.hpp"
46
47#include "Stokhos_Multiply.hpp"
50#include "Teuchos_ParameterList.hpp"
51
52
53//----------------------------------------------------------------------------
54//----------------------------------------------------------------------------
55
56namespace Stokhos {
57
61template< typename ValueType , class ExecutionSpace >
63public:
64
65 typedef ExecutionSpace execution_space ;
66 typedef typename execution_space::size_type size_type ;
67 typedef ValueType value_type ;
68
69private:
70
71 typedef Kokkos::View< size_type[] , execution_space > coord_array_type ;
72 typedef Kokkos::View< value_type[], execution_space > value_array_type ;
73 typedef Kokkos::View< size_type[], execution_space > entry_array_type ;
74 typedef Kokkos::View< size_type[], execution_space > row_map_array_type ;
75
85
86public:
87
88 inline
90
91 inline
93 m_k_coord() ,
94 m_j_coord() ,
95 m_value() ,
96 m_num_k() ,
97 m_num_j() ,
98 m_k_row_map() ,
99 m_j_row_map() ,
100 m_nnz(0) ,
101 m_flops(0) {}
102
103 inline
105 m_k_coord( rhs.m_k_coord ) ,
106 m_j_coord( rhs.m_j_coord ) ,
107 m_value( rhs.m_value ) ,
108 m_num_k( rhs.m_num_k ) ,
109 m_num_j( rhs.m_num_j ) ,
110 m_k_row_map( rhs.m_k_row_map ) ,
111 m_j_row_map( rhs.m_j_row_map ) ,
112 m_nnz( rhs.m_nnz ) ,
113 m_flops( rhs.m_flops ) {}
114
115 inline
117 {
118 m_k_coord = rhs.m_k_coord ;
119 m_j_coord = rhs.m_j_coord ;
120 m_value = rhs.m_value ;
121 m_num_k = rhs.m_num_k ;
122 m_num_j = rhs.m_num_j ;
123 m_k_row_map = rhs.m_k_row_map ;
124 m_j_row_map = rhs.m_j_row_map ;
125 m_nnz = rhs.m_nnz;
126 m_flops = rhs.m_flops;
127 return *this ;
128 }
129
131 KOKKOS_INLINE_FUNCTION
132 size_type dimension() const { return m_k_row_map.extent(0) - 1 ; }
133
135 KOKKOS_INLINE_FUNCTION
137 { return m_j_coord.extent(0); }
138
140 KOKKOS_INLINE_FUNCTION
142 { return m_k_row_map[i]; }
143
145 KOKKOS_INLINE_FUNCTION
147 { return m_k_row_map[i] + m_num_k(i); }
148
150 KOKKOS_INLINE_FUNCTION
152 { return m_num_k(i); }
153
155 KOKKOS_INLINE_FUNCTION
156 const size_type& k_coord( const size_type kEntry ) const
157 { return m_k_coord( kEntry ); }
158
160 KOKKOS_INLINE_FUNCTION
162 { return m_j_row_map[kEntry]; }
163
165 KOKKOS_INLINE_FUNCTION
166 size_type j_end( size_type kEntry ) const
167 { return m_j_row_map[kEntry] + m_num_j(kEntry); }
168
170 KOKKOS_INLINE_FUNCTION
171 size_type num_j( size_type kEntry ) const
172 { return m_num_j(kEntry); }
173
175 KOKKOS_INLINE_FUNCTION
176 const size_type& j_coord( const size_type jEntry ) const
177 { return m_j_coord( jEntry ); }
178
180 KOKKOS_INLINE_FUNCTION
181 const value_type & value( const size_type jEntry ) const
182 { return m_value( jEntry ); }
183
185 KOKKOS_INLINE_FUNCTION
187 { return m_nnz; }
188
190 KOKKOS_INLINE_FUNCTION
192 { return m_flops; }
193
194 template <typename OrdinalType>
195 static FlatSparse3Tensor
198 const Teuchos::ParameterList& params = Teuchos::ParameterList())
199 {
201
202 // Compute number of k's for each i
203 const size_type dimension = basis.size();
204 std::vector< size_t > k_coord_work( dimension , (size_t) 0 );
205 size_type k_entry_count = 0 ;
206 for (typename Cijk_type::i_iterator i_it=Cijk.i_begin();
207 i_it!=Cijk.i_end(); ++i_it) {
208 OrdinalType i = index(i_it);
209 k_coord_work[i] = Cijk.num_k(i_it);
210 k_entry_count += Cijk.num_k(i_it);
211 }
212
213 // Compute number of j's for each i and k
214 std::vector< size_t > j_coord_work( k_entry_count , (size_t) 0 );
215 size_type j_entry_count = 0 ;
216 size_type k_entry = 0 ;
217 for (typename Cijk_type::i_iterator i_it=Cijk.i_begin();
218 i_it!=Cijk.i_end(); ++i_it) {
219 for (typename Cijk_type::ik_iterator k_it = Cijk.k_begin(i_it);
220 k_it != Cijk.k_end(i_it); ++k_it, ++k_entry) {
221 OrdinalType k = index(k_it);
222 for (typename Cijk_type::ikj_iterator j_it = Cijk.j_begin(k_it);
223 j_it != Cijk.j_end(k_it); ++j_it) {
224 OrdinalType j = index(j_it);
225 if (j >= k) {
226 ++j_coord_work[k_entry];
227 ++j_entry_count;
228 }
229 }
230 }
231 }
232
233 /*
234 // Pad each row to have size divisible by alignment size
235 enum { Align = std::is_same<ExecutionSpace,Kokkos::Cuda>::value ? 32 : 2 };
236 for ( size_type i = 0 ; i < dimension ; ++i ) {
237 const size_t rem = coord_work[i] % Align;
238 if (rem > 0) {
239 const size_t pad = Align - rem;
240 coord_work[i] += pad;
241 entry_count += pad;
242 }
243 }
244 */
245
246 // Allocate tensor data
247 FlatSparse3Tensor tensor ;
248 tensor.m_k_coord = coord_array_type( "k_coord" , k_entry_count );
249 tensor.m_j_coord = coord_array_type( "j_coord" , j_entry_count );
250 tensor.m_value = value_array_type( "value" , j_entry_count );
251 tensor.m_num_k = entry_array_type( "num_k" , dimension );
252 tensor.m_num_j = entry_array_type( "num_j" , k_entry_count );
253 tensor.m_k_row_map = row_map_array_type( "k_row_map" , dimension+1 );
254 tensor.m_j_row_map = row_map_array_type( "j_row_map" , k_entry_count+1 );
255
256 // Create mirror, is a view if is host memory
257 typename coord_array_type::HostMirror
258 host_k_coord = Kokkos::create_mirror_view( tensor.m_k_coord );
259 typename coord_array_type::HostMirror
260 host_j_coord = Kokkos::create_mirror_view( tensor.m_j_coord );
261 typename value_array_type::HostMirror
262 host_value = Kokkos::create_mirror_view( tensor.m_value );
263 typename entry_array_type::HostMirror
264 host_num_k = Kokkos::create_mirror_view( tensor.m_num_k );
265 typename entry_array_type::HostMirror
266 host_num_j = Kokkos::create_mirror_view( tensor.m_num_j );
267 typename entry_array_type::HostMirror
268 host_k_row_map = Kokkos::create_mirror_view( tensor.m_k_row_map );
269 typename entry_array_type::HostMirror
270 host_j_row_map = Kokkos::create_mirror_view( tensor.m_j_row_map );
271
272 // Compute k row map
273 size_type sum = 0;
274 host_k_row_map(0) = 0;
275 for ( size_type i = 0 ; i < dimension ; ++i ) {
276 sum += k_coord_work[i];
277 host_k_row_map(i+1) = sum;
278 host_num_k(i) = 0;
279 }
280
281 // Compute j row map
282 sum = 0;
283 host_j_row_map(0) = 0;
284 for ( size_type i = 0 ; i < k_entry_count ; ++i ) {
285 sum += j_coord_work[i];
286 host_j_row_map(i+1) = sum;
287 host_num_j(i) = 0;
288 }
289
290 for ( size_type i = 0 ; i < dimension ; ++i ) {
291 k_coord_work[i] = host_k_row_map[i];
292 }
293 for ( size_type i = 0 ; i < k_entry_count ; ++i ) {
294 j_coord_work[i] = host_j_row_map[i];
295 }
296
297 for (typename Cijk_type::i_iterator i_it=Cijk.i_begin();
298 i_it!=Cijk.i_end(); ++i_it) {
299 OrdinalType i = index(i_it);
300 for (typename Cijk_type::ik_iterator k_it = Cijk.k_begin(i_it);
301 k_it != Cijk.k_end(i_it); ++k_it) {
302 OrdinalType k = index(k_it);
303 const size_type kEntry = k_coord_work[i];
304 ++k_coord_work[i];
305 host_k_coord(kEntry) = k ;
306 ++host_num_k(i);
307 for (typename Cijk_type::ikj_iterator j_it = Cijk.j_begin(k_it);
308 j_it != Cijk.j_end(k_it); ++j_it) {
309 OrdinalType j = index(j_it);
310 ValueType c = Stokhos::value(j_it);
311 if (j >= k) {
312 const size_type jEntry = j_coord_work[kEntry];
313 ++j_coord_work[kEntry];
314 host_value(jEntry) = (j != k) ? c : 0.5*c;
315 host_j_coord(jEntry) = j ;
316 ++host_num_j(kEntry);
317 ++tensor.m_nnz;
318 }
319 }
320 }
321 }
322
323 // Copy data to device if necessary
324 Kokkos::deep_copy( tensor.m_k_coord , host_k_coord );
325 Kokkos::deep_copy( tensor.m_j_coord , host_j_coord );
326 Kokkos::deep_copy( tensor.m_value , host_value );
327 Kokkos::deep_copy( tensor.m_num_k , host_num_k );
328 Kokkos::deep_copy( tensor.m_num_j , host_num_j );
329 Kokkos::deep_copy( tensor.m_k_row_map , host_k_row_map );
330 Kokkos::deep_copy( tensor.m_j_row_map , host_j_row_map );
331
332 tensor.m_flops = 5*tensor.m_nnz + dimension;
333
334 return tensor ;
335 }
336};
337
338template< class Device , typename OrdinalType , typename ValueType >
339FlatSparse3Tensor<ValueType, Device>
343 const Teuchos::ParameterList& params = Teuchos::ParameterList() )
344{
345 return FlatSparse3Tensor<ValueType, Device>::create( basis, Cijk, params );
346}
347
348template <typename ValueType, typename Device>
349class BlockMultiply< FlatSparse3Tensor< ValueType , Device > >
350{
351public:
352
353 typedef typename Device::size_type size_type ;
354 typedef FlatSparse3Tensor< ValueType , Device > tensor_type ;
355
356 template< typename MatrixValue , typename VectorValue >
357 KOKKOS_INLINE_FUNCTION
358 static void apply( const tensor_type & tensor ,
359 const MatrixValue * const a ,
360 const VectorValue * const x ,
361 VectorValue * const y )
362 {
363
364 const size_type nDim = tensor.dimension();
365
366 // Loop over i
367 for ( size_type i = 0; i < nDim; ++i) {
368 VectorValue ytmp = 0;
369
370 // Loop over k for this i
371 const size_type nk = tensor.num_k(i);
372 const size_type kBeg = tensor.k_begin(i);
373 const size_type kEnd = kBeg + nk;
374 for (size_type kEntry = kBeg; kEntry < kEnd; ++kEntry) {
375 const size_type k = tensor.k_coord(kEntry);
376 const MatrixValue ak = a[k];
377 const VectorValue xk = x[k];
378
379 // Loop over j for this i,k
380 const size_type nj = tensor.num_j(kEntry);
381 const size_type jBeg = tensor.j_begin(kEntry);
382 const size_type jEnd = jBeg + nj;
383 for (size_type jEntry = jBeg; jEntry < jEnd; ++jEntry) {
384 const size_type j = tensor.j_coord(jEntry);
385 ytmp += tensor.value(jEntry) * ( a[j] * xk + ak * x[j] );
386 }
387 }
388
389 y[i] += ytmp ;
390 }
391 }
392
393 KOKKOS_INLINE_FUNCTION
394 static size_type matrix_size( const tensor_type & tensor )
395 { return tensor.dimension(); }
396
397 KOKKOS_INLINE_FUNCTION
398 static size_type vector_size( const tensor_type & tensor )
399 { return tensor.dimension(); }
400};
401
402} /* namespace Stokhos */
403
404//----------------------------------------------------------------------------
405//----------------------------------------------------------------------------
406
407#endif /* #ifndef STOKHOS_FLAT_SPARSE_3_TENSOR_HPP */
static KOKKOS_INLINE_FUNCTION size_type vector_size(const tensor_type &tensor)
static KOKKOS_INLINE_FUNCTION size_type matrix_size(const tensor_type &tensor)
static KOKKOS_INLINE_FUNCTION void apply(const tensor_type &tensor, const MatrixValue *const a, const VectorValue *const x, VectorValue *const y)
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
KOKKOS_INLINE_FUNCTION size_type entry_count() const
Number of sparse entries.
Kokkos::View< value_type[], execution_space > value_array_type
KOKKOS_INLINE_FUNCTION const size_type & k_coord(const size_type kEntry) const
k coordinate for k entry 'kEntry'
FlatSparse3Tensor & operator=(const FlatSparse3Tensor &rhs)
KOKKOS_INLINE_FUNCTION size_type j_end(size_type kEntry) const
End j entries with a k entry 'kEntry'.
KOKKOS_INLINE_FUNCTION size_type j_begin(size_type kEntry) const
Begin j entries with a k entry 'kEntry'.
KOKKOS_INLINE_FUNCTION size_type num_j(size_type kEntry) const
Number of j entries with a k entry 'kEntry'.
KOKKOS_INLINE_FUNCTION size_type num_flops() const
Number flop's per multiply-add.
KOKKOS_INLINE_FUNCTION size_type num_k(size_type i) const
Number of k entries with a coordinate 'i'.
KOKKOS_INLINE_FUNCTION const value_type & value(const size_type jEntry) const
Value for j entry 'jEntry'.
KOKKOS_INLINE_FUNCTION size_type dimension() const
Dimension of the tensor.
KOKKOS_INLINE_FUNCTION const size_type & j_coord(const size_type jEntry) const
j coordinate for j entry 'jEntry'
static FlatSparse3Tensor create(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())
KOKKOS_INLINE_FUNCTION size_type k_begin(size_type i) const
Begin k entries with a coordinate 'i'.
Kokkos::View< size_type[], execution_space > row_map_array_type
FlatSparse3Tensor(const FlatSparse3Tensor &rhs)
execution_space::size_type size_type
Kokkos::View< size_type[], execution_space > coord_array_type
Kokkos::View< size_type[], execution_space > entry_array_type
KOKKOS_INLINE_FUNCTION size_type num_non_zeros() const
Number of non-zero's.
KOKKOS_INLINE_FUNCTION size_type k_end(size_type i) const
End k entries with a coordinate 'i'.
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.
ordinal_type num_k() const
Number of k entries in C(i,j,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.
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
Top-level namespace for Stokhos classes and functions.
FlatSparse3Tensor< ValueType, Device > create_flat_sparse_3_tensor(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())