Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Stokhos_LinearSparse3Tensor.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_LINEAR_SPARSE_3_TENSOR_HPP
43#define STOKHOS_LINEAR_SPARSE_3_TENSOR_HPP
44
45#include "Kokkos_Core.hpp"
46
47#include "Stokhos_Multiply.hpp"
50#include "Teuchos_ParameterList.hpp"
51#include "Stokhos_TinyVec.hpp"
52
53//----------------------------------------------------------------------------
54//----------------------------------------------------------------------------
55
56namespace Stokhos {
57
61template< typename ValueType , class ExecutionSpace , int BlockSize >
63public:
64
65 typedef ExecutionSpace execution_space ;
66 typedef typename execution_space::size_type size_type ;
67 typedef ValueType value_type ;
68
69 static const int block_size = BlockSize;
70
71private:
72
73 typedef Kokkos::View< value_type[], execution_space > value_array_type ;
74
81
82public:
83
84 inline
86
87 inline
89 m_value() ,
90 m_dim() ,
92 m_nnz(0) ,
93 m_flops(0) ,
94 m_symmetric(false) {}
95
96 inline
98 m_value( rhs.m_value ) ,
99 m_dim( rhs.m_dim ),
101 m_nnz( rhs.m_nnz ) ,
102 m_flops( rhs.m_flops ) ,
103 m_symmetric( rhs.m_symmetric ) {}
104
105 inline
107 {
108 m_value = rhs.m_value ;
109 m_dim = rhs.m_dim ;
110 m_aligned_dim = rhs.m_aligned_dim;
111 m_nnz = rhs.m_nnz;
112 m_flops = rhs.m_flops;
113 m_symmetric = rhs.m_symmetric;
114 return *this ;
115 }
116
118 KOKKOS_INLINE_FUNCTION
119 size_type dimension() const { return m_dim ; }
120
122 KOKKOS_INLINE_FUNCTION
124
126 KOKKOS_INLINE_FUNCTION
128 { return m_value.extent(0); }
129
131 KOKKOS_INLINE_FUNCTION
132 bool symmetric() const
133 { return m_symmetric; }
134
136 KOKKOS_INLINE_FUNCTION
137 const value_type & value( const size_type entry ) const
138 { return m_value( entry ); }
139
141 KOKKOS_INLINE_FUNCTION
143 { return m_nnz; }
144
146 KOKKOS_INLINE_FUNCTION
148 { return m_flops; }
149
150 template <typename OrdinalType>
154 const Teuchos::ParameterList& params)
155 {
156 const bool symmetric = params.get<bool>("Symmetric");
157
158 // Allocate tensor data -- currently assuming isotropic
159 const size_type dim = basis.size();
160 LinearSparse3Tensor tensor ;
161 tensor.m_dim = dim;
162 tensor.m_aligned_dim = dim;
163 if (tensor.m_aligned_dim % block_size)
165 tensor.m_symmetric = symmetric;
166 tensor.m_nnz = symmetric ? 2 : 3 ;
167 tensor.m_value = value_array_type( "value" , tensor.m_nnz );
168
169 // Create mirror, is a view if is host memory
170 typename value_array_type::HostMirror
171 host_value = Kokkos::create_mirror_view( tensor.m_value );
172
173 // Get Cijk values
174 Teuchos::Array< Teuchos::RCP<const Stokhos::OneDOrthogPolyBasis<OrdinalType,ValueType> > > bases = basis.getCoordinateBases();
175 Teuchos::RCP< Stokhos::Dense3Tensor<OrdinalType,ValueType> > cijk =
176 bases[0]->computeTripleProductTensor();
177 // For non-isotropic, need to take products of these over basis components
178 host_value(0) = (*cijk)(0,0,0);
179 host_value(1) = (*cijk)(0,1,1);
180 if (!symmetric)
181 host_value(2) = (*cijk)(1,1,1);
182
183 // Copy data to device if necessary
184 Kokkos::deep_copy( tensor.m_value , host_value );
185
186 tensor.m_flops = 8*dim;
187 if (!symmetric)
188 tensor.m_flops += 2*dim ;
189
190 return tensor ;
191 }
192};
193
194template< class Device , typename OrdinalType , typename ValueType , int BlockSize >
195LinearSparse3Tensor<ValueType, Device,BlockSize>
199 const Teuchos::ParameterList& params)
200{
202 basis, Cijk, params );
203}
204
205template < typename ValueType, typename Device, int BlockSize >
206class BlockMultiply< LinearSparse3Tensor< ValueType , Device , BlockSize > >
207{
208public:
209
210 typedef typename Device::size_type size_type ;
211 typedef LinearSparse3Tensor< ValueType , Device , BlockSize > tensor_type ;
212
213 template< typename MatrixValue , typename VectorValue >
214 KOKKOS_INLINE_FUNCTION
215 static void apply( const tensor_type & tensor ,
216 const MatrixValue * const a ,
217 const VectorValue * const x ,
218 VectorValue * const y )
219 {
220 const size_type block_size = tensor_type::block_size;
221 typedef TinyVec<ValueType,block_size,true> TV;
222 const size_type dim = tensor.dimension();
223
224 const ValueType c0 = tensor.value(0);
225 const ValueType c1 = tensor.value(1);
226 const ValueType a0 = a[0];
227 const ValueType x0 = x[0];
228
229 if (block_size > 1) {
230
231 TV vc0(c0), vc1(c1), va0(a0), vx0(x0), vy0;
232 TV ai, ai2, xi, yi;
233
234 const MatrixValue *aa = a;
235 const VectorValue *xx = x;
236 VectorValue *yy = y;
237 vy0.zero();
238
239 const size_type nBlock = dim / block_size;
240 const size_type iEnd = nBlock * block_size;
241
242 if (tensor.symmetric()) {
243
244 size_type i=0;
245 for ( ; i < iEnd; i+=block_size,aa+=block_size,xx+=block_size,yy+=block_size) {
246 ai.aligned_load(aa);
247 ai2 = ai;
248 xi.aligned_load(xx);
249 yi.aligned_load(yy);
250
251 // y[i] += c1*(a0*xi + ai*x0);
252 ai.times_equal(vx0);
253 ai2.times_equal(xi);
254 xi.times_equal(va0);
255 xi.plus_equal(ai);
256 xi.times_equal(vc1);
257 yi.plus_equal(xi);
258 yi.aligned_scatter(yy);
259
260 // y0 += c1*ai*xi;
261 ai2.times_equal(vc1);
262 vy0.plus_equal(ai2);
263 }
264 ValueType y0 = vy0.sum();
265
266 // Do remaining entries with a scalar loop
267 for ( ; i < dim; ++i) {
268 const ValueType ai = *aa++;
269 const ValueType xi = *xx++;
270 *yy++ += c1*(a0*xi + ai*x0);
271 y0 += c1*ai*xi;
272 }
273 y[0] += y0 + (c0-3.0*c1)*a0*x0;
274 }
275 else {
276
277 const ValueType c2 = tensor.value(2);
278 TV vc2(c2);
279 size_type i=0;
280 for ( ; i < iEnd; i+=block_size,aa+=block_size,xx+=block_size,yy+=block_size) {
281 ai.aligned_load(aa);
282 ai2 = ai;
283 xi.aligned_load(xx);
284 yi.aligned_load(yy);
285
286 // y[i] += c1*(a0*xi + ai*x0) + c2*aixi;
287 ai.times_equal(vx0);
288 ai2.times_equal(xi);
289 xi.times_equal(va0);
290 xi.plus_equal(ai);
291 xi.times_equal(vc1);
292 yi.plus_equal(xi);
293 ai = ai2;
294 ai.times_equal(vc2);
295 yi.plus_equal(ai);
296 yi.aligned_scatter(yy);
297
298 // y0 += c1*aixi;
299 ai2.times_equal(vc1);
300 vy0.plus_equal(ai2);
301 }
302 ValueType y0 = vy0.sum();
303
304 // Do remaining entries with a scalar loop
305 for ( ; i < dim; ++i) {
306 const ValueType ai = *aa++;
307 const ValueType xi = *xx++;
308 const ValueType aixi = ai*xi;
309 *yy++ += c1*(a0*xi + ai*x0) + c2*aixi;
310 y0 += c1*aixi;
311 }
312 y[0] += y0 + (c0-3.0*c1-c2)*a0*x0;
313
314 }
315
316 }
317
318 else {
319
320 ValueType y0 = c0*a0*x0;
321
322 if (tensor.symmetric()) {
323
324 for ( size_type i = 1; i < dim; ++i) {
325 const ValueType ai = a[i];
326 const ValueType xi = x[i];
327 y[i] += c1*(a0*xi + ai*x0);
328 y0 += c1*ai*xi;
329 }
330 y[0] += y0;
331
332 }
333 else {
334
335 const ValueType c2 = tensor.value(2);
336 for ( size_type i = 1; i < dim; ++i) {
337 const ValueType ai = a[i];
338 const ValueType xi = x[i];
339 const ValueType aixi = ai*xi;
340 y[i] += c1*(a0*xi + ai*x0) + c2*aixi;
341 y0 += c1*aixi;
342 }
343 y[0] += y0;
344
345 }
346
347 }
348
349 }
350
351 KOKKOS_INLINE_FUNCTION
352 static size_type matrix_size( const tensor_type & tensor )
353 { return tensor.dimension(); }
354
355 KOKKOS_INLINE_FUNCTION
356 static size_type vector_size( const tensor_type & tensor )
357 { return tensor.dimension(); }
358};
359
360} /* namespace Stokhos */
361
362//----------------------------------------------------------------------------
363//----------------------------------------------------------------------------
364
365#endif /* #ifndef STOKHOS_LINEAR_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.
LinearSparse3Tensor & operator=(const LinearSparse3Tensor &rhs)
static LinearSparse3Tensor create(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params)
KOKKOS_INLINE_FUNCTION const value_type & value(const size_type entry) const
Value for entry 'entry'.
KOKKOS_INLINE_FUNCTION size_type dimension() const
Dimension of the tensor.
KOKKOS_INLINE_FUNCTION size_type entry_count() const
Number of sparse entries.
LinearSparse3Tensor(const LinearSparse3Tensor &rhs)
KOKKOS_INLINE_FUNCTION size_type aligned_dimension() const
Dimension of the tensor.
Kokkos::View< value_type[], execution_space > value_array_type
KOKKOS_INLINE_FUNCTION bool symmetric() const
Is tensor built from symmetric PDFs.
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.
virtual ordinal_type size() const =0
Return total size of basis.
Abstract base class for multivariate orthogonal polynomials generated from tensor products of univari...
virtual Teuchos::Array< Teuchos::RCP< const OneDOrthogPolyBasis< ordinal_type, value_type > > > getCoordinateBases() const =0
Return array of coordinate bases.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
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.
LinearSparse3Tensor< ValueType, Device, BlockSize > create_linear_sparse_3_tensor(const Stokhos::ProductBasis< OrdinalType, ValueType > &basis, const Stokhos::Sparse3Tensor< OrdinalType, ValueType > &Cijk, const Teuchos::ParameterList &params)