Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_FlatSparse3Tensor_kji.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_KJI_HPP
43#define STOKHOS_FLAT_SPARSE_3_TENSOR_KJI_HPP
44
45#include "Kokkos_Core.hpp"
46
47#include "Stokhos_Multiply.hpp"
50#include "Teuchos_ParameterList.hpp"
51
52//----------------------------------------------------------------------------
53//----------------------------------------------------------------------------
54
55namespace Stokhos {
56
60template< typename ValueType , class ExecutionSpace >
62public:
63
64 typedef ExecutionSpace execution_space ;
65 typedef typename execution_space::size_type size_type ;
66 typedef ValueType value_type ;
67
68private:
69
70 typedef Kokkos::View< size_type[] , execution_space > coord_array_type ;
71 typedef Kokkos::View< value_type[], execution_space > value_array_type ;
72 typedef Kokkos::View< size_type[], execution_space > entry_array_type ;
73 typedef Kokkos::View< size_type[], execution_space > row_map_array_type ;
74
85
86
87public:
88
89 inline
91
92 inline
94 m_j_coord() ,
95 m_i_coord() ,
96 m_value() ,
97 m_num_j() ,
98 m_num_i() ,
99 m_j_row_map() ,
100 m_i_row_map() ,
101 m_nnz(0) ,
102 m_dim(0) ,
103 m_flops(0) {}
104
105 inline
107 m_j_coord( rhs.m_j_coord ) ,
108 m_i_coord( rhs.m_i_coord ) ,
109 m_value( rhs.m_value ) ,
110 m_num_j( rhs.m_num_j ) ,
111 m_num_i( rhs.m_num_i ) ,
112 m_j_row_map( rhs.m_j_row_map ) ,
113 m_i_row_map( rhs.m_i_row_map ) ,
114 m_nnz( rhs.m_nnz ) ,
115 m_dim( rhs.m_dim ) ,
116 m_flops( rhs.m_flops ) {}
117
118 inline
120 {
121 m_j_coord = rhs.m_j_coord ;
122 m_i_coord = rhs.m_i_coord ;
123 m_value = rhs.m_value ;
124 m_num_j = rhs.m_num_j ;
125 m_num_i = rhs.m_num_i ;
128 m_nnz = rhs.m_nnz;
129 m_dim = rhs.m_dim ;
130 m_flops = rhs.m_flops ;
131 return *this ;
132 }
133
135 KOKKOS_INLINE_FUNCTION
136 size_type dimension() const { return m_dim ; }
137
139 KOKKOS_INLINE_FUNCTION
140 size_type num_k() const { return m_j_row_map.extent(0) - 1 ; }
141
143 KOKKOS_INLINE_FUNCTION
145 { return m_i_coord.extent(0); }
146
148 KOKKOS_INLINE_FUNCTION
150 { return m_j_row_map[k]; }
151
153 KOKKOS_INLINE_FUNCTION
155 { return m_j_row_map[k] + m_num_j(k); }
156
158 KOKKOS_INLINE_FUNCTION
160 { return m_num_j(k); }
161
163 KOKKOS_INLINE_FUNCTION
164 const size_type& j_coord( const size_type jEntry ) const
165 { return m_j_coord( jEntry ); }
166
168 KOKKOS_INLINE_FUNCTION
170 { return m_i_row_map[jEntry]; }
171
173 KOKKOS_INLINE_FUNCTION
174 size_type i_end( size_type jEntry ) const
175 { return m_i_row_map[jEntry] + m_num_i(jEntry); }
176
178 KOKKOS_INLINE_FUNCTION
179 size_type num_i( size_type jEntry ) const
180 { return m_num_i(jEntry); }
181
183 KOKKOS_INLINE_FUNCTION
184 const size_type& i_coord( const size_type iEntry ) const
185 { return m_i_coord( iEntry ); }
186
188 KOKKOS_INLINE_FUNCTION
189 const value_type & value( const size_type iEntry ) const
190 { return m_value( iEntry ); }
191
193 KOKKOS_INLINE_FUNCTION
195 { return m_nnz; }
196
198 KOKKOS_INLINE_FUNCTION
200 { return m_flops; }
201
202 template <typename OrdinalType>
206 const Teuchos::ParameterList& params = Teuchos::ParameterList())
207 {
209
210 // Compute number of j's for each k
211 const size_type dimension = basis.size();
212 const size_type nk = Cijk.num_k();
213 std::vector< size_t > j_coord_work( nk , (size_t) 0 );
214 size_type j_entry_count = 0 ;
215 for (typename Cijk_type::k_iterator k_it=Cijk.k_begin();
216 k_it!=Cijk.k_end(); ++k_it) {
217 OrdinalType k = index(k_it);
218 for (typename Cijk_type::kj_iterator j_it = Cijk.j_begin(k_it);
219 j_it != Cijk.j_end(k_it); ++j_it) {
220 OrdinalType j = index(j_it);
221 if (j >= k) {
222 ++j_coord_work[k];
223 ++j_entry_count;
224 }
225 }
226 }
227
228 // Compute number of i's for each k and j
229 std::vector< size_t > i_coord_work( j_entry_count , (size_t) 0 );
230 size_type i_entry_count = 0 ;
231 size_type j_entry = 0 ;
232 for (typename Cijk_type::k_iterator k_it=Cijk.k_begin();
233 k_it!=Cijk.k_end(); ++k_it) {
234 OrdinalType k = index(k_it);
235 for (typename Cijk_type::kj_iterator j_it = Cijk.j_begin(k_it);
236 j_it != Cijk.j_end(k_it); ++j_it) {
237 OrdinalType j = index(j_it);
238 if (j >= k) {
239 for (typename Cijk_type::kji_iterator i_it = Cijk.i_begin(j_it);
240 i_it != Cijk.i_end(j_it); ++i_it) {
241 ++i_coord_work[j_entry];
242 ++i_entry_count;
243 }
244 ++j_entry;
245 }
246 }
247 }
248
249 /*
250 // Pad each row to have size divisible by alignment size
251 enum { Align = std::is_same<ExecutionSpace,Kokkos::Cuda>::value ? 32 : 2 };
252 for ( size_type i = 0 ; i < dimension ; ++i ) {
253 const size_t rem = coord_work[i] % Align;
254 if (rem > 0) {
255 const size_t pad = Align - rem;
256 coord_work[i] += pad;
257 entry_count += pad;
258 }
259 }
260 */
261
262 // Allocate tensor data
263 FlatSparse3Tensor_kji tensor ;
264 tensor.m_dim = dimension;
265 tensor.m_j_coord = coord_array_type( "j_coord" , j_entry_count );
266 tensor.m_i_coord = coord_array_type( "i_coord" , i_entry_count );
267 tensor.m_value = value_array_type( "value" , i_entry_count );
268 tensor.m_num_j = entry_array_type( "num_j" , nk );
269 tensor.m_num_i = entry_array_type( "num_i" , j_entry_count );
270 tensor.m_j_row_map = row_map_array_type( "j_row_map" , nk+1 );
271 tensor.m_i_row_map = row_map_array_type( "i_row_map" , j_entry_count+1 );
272 tensor.m_flops = 3*j_entry_count + 2*i_entry_count;
273
274 // Create mirror, is a view if is host memory
275 typename coord_array_type::HostMirror
276 host_j_coord = Kokkos::create_mirror_view( tensor.m_j_coord );
277 typename coord_array_type::HostMirror
278 host_i_coord = Kokkos::create_mirror_view( tensor.m_i_coord );
279 typename value_array_type::HostMirror
280 host_value = Kokkos::create_mirror_view( tensor.m_value );
281 typename entry_array_type::HostMirror
282 host_num_j = Kokkos::create_mirror_view( tensor.m_num_j );
283 typename entry_array_type::HostMirror
284 host_num_i = Kokkos::create_mirror_view( tensor.m_num_i );
285 typename entry_array_type::HostMirror
286 host_j_row_map = Kokkos::create_mirror_view( tensor.m_j_row_map );
287 typename entry_array_type::HostMirror
288 host_i_row_map = Kokkos::create_mirror_view( tensor.m_i_row_map );
289
290 // Compute j row map
291 size_type sum = 0;
292 host_j_row_map(0) = 0;
293 for ( size_type k = 0 ; k < nk ; ++k ) {
294 sum += j_coord_work[k];
295 host_j_row_map(k+1) = sum;
296 host_num_j(k) = 0;
297 }
298
299 // Compute i row map
300 sum = 0;
301 host_i_row_map(0) = 0;
302 for ( size_type j = 0 ; j < j_entry_count ; ++j ) {
303 sum += i_coord_work[j];
304 host_i_row_map(j+1) = sum;
305 host_num_i(j) = 0;
306 }
307
308 for ( size_type k = 0 ; k < nk ; ++k ) {
309 j_coord_work[k] = host_j_row_map[k];
310 }
311 for ( size_type j = 0 ; j < j_entry_count ; ++j ) {
312 i_coord_work[j] = host_i_row_map[j];
313 }
314
315 for (typename Cijk_type::k_iterator k_it=Cijk.k_begin();
316 k_it!=Cijk.k_end(); ++k_it) {
317 OrdinalType k = index(k_it);
318 for (typename Cijk_type::kj_iterator j_it = Cijk.j_begin(k_it);
319 j_it != Cijk.j_end(k_it); ++j_it) {
320 OrdinalType j = index(j_it);
321 if (j >= k) {
322 const size_type jEntry = j_coord_work[k];
323 ++j_coord_work[k];
324 host_j_coord(jEntry) = j ;
325 ++host_num_j(k);
326 for (typename Cijk_type::kji_iterator i_it = Cijk.i_begin(j_it);
327 i_it != Cijk.i_end(j_it); ++i_it) {
328 OrdinalType i = index(i_it);
329 ValueType c = Stokhos::value(i_it);
330 const size_type iEntry = i_coord_work[jEntry];
331 ++i_coord_work[jEntry];
332 host_value(iEntry) = (j != k) ? c : 0.5*c;
333 host_i_coord(iEntry) = i ;
334 ++host_num_i(jEntry);
335 ++tensor.m_nnz;
336 }
337 }
338 }
339 }
340
341 // Copy data to device if necessary
342 Kokkos::deep_copy( tensor.m_j_coord , host_j_coord );
343 Kokkos::deep_copy( tensor.m_i_coord , host_i_coord );
344 Kokkos::deep_copy( tensor.m_value , host_value );
345 Kokkos::deep_copy( tensor.m_num_j , host_num_j );
346 Kokkos::deep_copy( tensor.m_num_i , host_num_i );
347 Kokkos::deep_copy( tensor.m_j_row_map , host_j_row_map );
348 Kokkos::deep_copy( tensor.m_i_row_map , host_i_row_map );
349
350 return tensor ;
351 }
352};
353
354template< class Device , typename OrdinalType , typename ValueType >
355FlatSparse3Tensor_kji<ValueType, Device>
359 const Teuchos::ParameterList& params = Teuchos::ParameterList())
360{
362 basis, Cijk, params );
363}
364
365template < typename ValueType, typename Device >
366class BlockMultiply< FlatSparse3Tensor_kji< ValueType , Device > >
367{
368public:
369
370 typedef typename Device::size_type size_type ;
372
373 template< typename MatrixValue , typename VectorValue >
374 KOKKOS_INLINE_FUNCTION
375 static void apply( const tensor_type & tensor ,
376 const MatrixValue * const a ,
377 const VectorValue * const x ,
378 VectorValue * const y )
379 {
380 const size_type nk = tensor.num_k();
381
382 // Loop over k
383 for ( size_type k = 0; k < nk; ++k) {
384 const MatrixValue ak = a[k];
385 const VectorValue xk = x[k];
386
387 // Loop over j for this k
388 const size_type nj = tensor.num_j(k);
389 const size_type jBeg = tensor.j_begin(k);
390 const size_type jEnd = jBeg + nj;
391 for (size_type jEntry = jBeg; jEntry < jEnd; ++jEntry) {
392 const size_type j = tensor.j_coord(jEntry);
393 VectorValue tmp = a[j] * xk + ak * x[j];
394
395 // Loop over i for this k,j
396 const size_type ni = tensor.num_i(jEntry);
397 const size_type iBeg = tensor.i_begin(jEntry);
398 const size_type iEnd = iBeg + ni;
399 for (size_type iEntry = iBeg; iEntry < iEnd; ++iEntry) {
400 const size_type i = tensor.i_coord(iEntry);
401 y[i] += tensor.value(iEntry) * tmp;
402 }
403 }
404 }
405 }
406
407 KOKKOS_INLINE_FUNCTION
408 static size_type matrix_size( const tensor_type & tensor )
409 { return tensor.dimension(); }
410
411 KOKKOS_INLINE_FUNCTION
412 static size_type vector_size( const tensor_type & tensor )
413 { return tensor.dimension(); }
414};
415
416} /* namespace Stokhos */
417
418//----------------------------------------------------------------------------
419//----------------------------------------------------------------------------
420
421#endif /* #ifndef STOKHOS_FLAT_SPARSE_3_TENSOR_KJI_HPP */
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)
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_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'
Kokkos::View< value_type[], execution_space > value_array_type
KOKKOS_INLINE_FUNCTION const value_type & value(const size_type iEntry) const
Value for i entry 'iEntry'.
KOKKOS_INLINE_FUNCTION size_type j_begin(size_type k) const
Begin j entries with a coordinate 'k'.
KOKKOS_INLINE_FUNCTION size_type num_k() const
Number of k entries.
KOKKOS_INLINE_FUNCTION size_type num_non_zeros() const
Number of non-zero's.
KOKKOS_INLINE_FUNCTION size_type num_flops() const
Number flop's per multiply-add.
FlatSparse3Tensor_kji & operator=(const FlatSparse3Tensor_kji &rhs)
KOKKOS_INLINE_FUNCTION size_type num_i(size_type jEntry) const
Number of i entries with a j entry 'jEntry'.
Kokkos::View< size_type[], execution_space > entry_array_type
KOKKOS_INLINE_FUNCTION size_type i_begin(size_type jEntry) const
Begin i entries with a j entry 'jEntry'.
static FlatSparse3Tensor_kji create(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())
KOKKOS_INLINE_FUNCTION const size_type & i_coord(const size_type iEntry) const
i coordinate for i entry 'iEntry'
Kokkos::View< size_type[], execution_space > coord_array_type
KOKKOS_INLINE_FUNCTION size_type i_end(size_type jEntry) const
End i entries with a j entry 'jEntry'.
Kokkos::View< size_type[], execution_space > row_map_array_type
KOKKOS_INLINE_FUNCTION size_type j_end(size_type k) const
End j entries with a coordinate 'k'.
FlatSparse3Tensor_kji(const FlatSparse3Tensor_kji &rhs)
KOKKOS_INLINE_FUNCTION size_type num_j(size_type k) const
Number of j entries with a coordinate 'k'.
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.
Top-level namespace for Stokhos classes and functions.
FlatSparse3Tensor_kji< ValueType, Device > create_flat_sparse_3_tensor_kji(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params=Teuchos::ParameterList())