Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
Kokkos_CrsMatrix_MP_Vector.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 KOKKOS_CRSMATRIX_MP_VECTOR_HPP
43 #define KOKKOS_CRSMATRIX_MP_VECTOR_HPP
44 
45 #include "Sacado_MP_Vector.hpp"
47 #include "Kokkos_Sparse.hpp"
48 #include "Kokkos_Blas1_MP_Vector.hpp" // for some utilities
49 
50 #include "Kokkos_Core.hpp"
51 #include "Stokhos_Multiply.hpp"
52 
53 #include "Teuchos_TestForException.hpp"
54 
55 //----------------------------------------------------------------------------
56 // Specializations of KokkosSparse::CrsMatrix for Sacado::MP::Vector scalar type
57 //----------------------------------------------------------------------------
58 
59 namespace Stokhos {
60 
61 namespace details {
62 
63 template <typename Matrix, typename InputVector, typename OutputVector,
64  typename Update = MultiplyAssign,
65  typename Enabled = void>
66 class MPMultiply {};
67 
68 // Kernel implementing y = A * x where
69 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
70 // x, y == Kokkos::View< Sacado::MP::Vector<...>*,...>,
71 // x and y are rank 1, any layout
72 // We spell everything out here to make sure the ranks and devices match.
73 //
74 // This implementation uses overloaded operators for MP::Vector.
75 template <typename MatrixDevice,
76  typename MatrixStorage,
77  typename MatrixOrdinal,
78  typename MatrixMemory,
79  typename MatrixSize,
80  typename InputStorage,
81  typename ... InputP,
82  typename OutputStorage,
83  typename ... OutputP,
84  typename Update>
85 class MPMultiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
86  MatrixOrdinal,
87  MatrixDevice,
88  MatrixMemory,
89  MatrixSize>,
90  Kokkos::View< Sacado::MP::Vector<InputStorage>*,
91  InputP... >,
92  Kokkos::View< Sacado::MP::Vector<OutputStorage>*,
93  OutputP... >,
94  Update
95 #ifdef KOKKOS_HAVE_CUDA
96  , typename std::enable_if<
97  !std::is_same<MatrixDevice,Kokkos::Cuda>::value >::type
98 #endif
99  >
100 {
101 
102 public:
103 
108 
109  typedef MatrixDevice execution_space;
110  typedef typename execution_space::size_type size_type;
111 
112  typedef KokkosSparse::CrsMatrix< MatrixValue,
113  MatrixOrdinal,
114  MatrixDevice,
115  MatrixMemory,
116  MatrixSize > matrix_type;
117  typedef Kokkos::View< InputVectorValue*,
118  InputP... > input_vector_type;
119  typedef Kokkos::View< OutputVectorValue*,
120  OutputP... > output_vector_type;
122 
127 
129  const input_vector_type & x,
130  const output_vector_type & y,
131  const update_type& update )
132  : m_A( A )
133  , m_x( x )
134  , m_y( y )
135  , m_update( update )
136  {}
137 
138  KOKKOS_INLINE_FUNCTION
139  void operator()( const size_type iRow ) const
140  {
141  // Compute mat-vec for this row
142  const size_type iEntryBegin = m_A.graph.row_map[iRow];
143  const size_type iEntryEnd = m_A.graph.row_map[iRow+1];
144  scalar_type sum = 0.0;
145  for (size_type iEntry = iEntryBegin; iEntry < iEntryEnd; ++iEntry) {
146  size_type iCol = m_A.graph.entries(iEntry);
147  sum += m_A.values(iEntry) * m_x(iCol);
148  }
149  m_update( m_y(iRow), sum );
150  } // operator()
151 
152  static void apply( const matrix_type & A,
153  const input_vector_type & x,
154  const output_vector_type & y,
155  const update_type & update )
156  {
157  const size_type row_count = A.graph.row_map.dimension_0()-1;
158  Kokkos::parallel_for( row_count, MPMultiply(A,x,y,update) );
159  }
160 };
161 
162 // Kernel implementing y = A * x where
163 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
164 // x, y == Kokkos::View< Sacado::MP::Vector<...>**,...>,
165 // x and y are rank 2, any layout
166 // We spell everything out here to make sure the ranks and devices match.
167 //
168 // This implementation uses overloaded operators for MP::Vector.
169 template <typename MatrixDevice,
170  typename MatrixStorage,
171  typename MatrixOrdinal,
172  typename MatrixMemory,
173  typename MatrixSize,
174  typename InputStorage,
175  typename ... InputP,
176  typename OutputStorage,
177  typename ... OutputP,
178  typename Update>
179 class MPMultiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
180  MatrixOrdinal,
181  MatrixDevice,
182  MatrixMemory,
183  MatrixSize >,
184  Kokkos::View< Sacado::MP::Vector<InputStorage>**,
185  InputP... >,
186  Kokkos::View< Sacado::MP::Vector<OutputStorage>**,
187  OutputP... >,
188  Update
189 #ifdef KOKKOS_HAVE_CUDA
190  , typename std::enable_if<
191  !std::is_same<MatrixDevice,Kokkos::Cuda>::value >::type
192 #endif
193  >
194 {
195 public:
200 
201  typedef MatrixDevice execution_space;
202  typedef typename execution_space::size_type size_type;
203 
204  typedef KokkosSparse::CrsMatrix< MatrixValue,
205  MatrixOrdinal,
206  MatrixDevice,
207  MatrixMemory,
208  MatrixSize > matrix_type;
209  typedef typename matrix_type::values_type matrix_values_type;
210  typedef Kokkos::View< InputVectorValue**,
211  InputP... > input_vector_type;
212  typedef Kokkos::View< OutputVectorValue**,
213  OutputP... > output_vector_type;
215 
220 
222  const input_vector_type & x,
223  const output_vector_type & y,
224  const update_type& update )
225  : m_A( A )
226  , m_x( x )
227  , m_y( y )
228  , m_update( update )
229  {}
230 
231  KOKKOS_INLINE_FUNCTION
232  void operator()( const size_type iRow ) const
233  {
235  // Loop over columns of x, y
236  const size_type num_col = m_y.dimension_1();
237  for (size_type col=0; col<num_col; ++col) {
238  // Compute mat-vec for this row
239  const size_type iEntryBegin = m_A.graph.row_map[iRow];
240  const size_type iEntryEnd = m_A.graph.row_map[iRow+1];
241  sum = 0.0;
242  for (size_type iEntry = iEntryBegin; iEntry < iEntryEnd; ++iEntry) {
243  size_type iCol = m_A.graph.entries(iEntry);
244  sum += m_A.values(iEntry) * m_x(iCol,col);
245  }
246  m_update( m_y(iRow,col), sum );
247  } // x, y column loop
248  } // operator()
249 
250 public:
251 
252  static void apply( const matrix_type & A,
253  const input_vector_type & x,
254  const output_vector_type & y,
255  const update_type & update )
256  {
257  const size_type row_count = A.graph.row_map.dimension_0()-1;
258  Kokkos::parallel_for( row_count, MPMultiply(A,x,y,update) );
259  }
260 };
261 
262 } // namespace details
263 
264 // Kernel implementing y = A * x where
265 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
266 // x, y == Kokkos::View< Sacado::MP::Vector<...>*,...>,
267 // x and y are rank 1, any layout
268 // We spell everything out here to make sure the ranks and devices match.
269 //
270 // This implementation uses overloaded operators for MP::Vector.
271 template <typename MatrixDevice,
272  typename MatrixStorage,
273  typename MatrixOrdinal,
274  typename MatrixMemory,
275  typename MatrixSize,
276  typename InputStorage,
277  typename ... InputP,
278  typename OutputStorage,
279  typename ... OutputP>
280 class Multiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
281  MatrixOrdinal,
282  MatrixDevice,
283  MatrixMemory,
284  MatrixSize >,
285  Kokkos::View< Sacado::MP::Vector<InputStorage>*,
286  InputP... >,
287  Kokkos::View< Sacado::MP::Vector<OutputStorage>*,
288  OutputP... >
289  >
290 {
291 public:
295 
296  typedef MatrixDevice execution_space;
297  typedef typename execution_space::size_type size_type;
298 
299  typedef KokkosSparse::CrsMatrix< MatrixValue,
300  MatrixOrdinal,
301  MatrixDevice,
302  MatrixMemory,
303  MatrixSize > matrix_type;
304  typedef typename matrix_type::values_type matrix_values_type;
305  typedef Kokkos::View< InputVectorValue*,
306  InputP... > input_vector_type;
307  typedef Kokkos::View< OutputVectorValue*,
308  OutputP... > output_vector_type;
309 
310 public:
311 
312  static void apply( const matrix_type & A,
313  const input_vector_type & x,
314  const output_vector_type & y )
315  {
317  multiply_type::apply(A,x,y,details::MultiplyAssign());
318  }
319 };
320 
321 // Kernel implementing y = A * x where
322 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
323 // x, y == Kokkos::View< Sacado::MP::Vector<...>**,...>,
324 // x and y are rank 2, any layout
325 // We spell everything out here to make sure the ranks and devices match.
326 //
327 // This implementation uses overloaded operators for MP::Vector.
328 template <typename MatrixDevice,
329  typename MatrixStorage,
330  typename MatrixOrdinal,
331  typename MatrixMemory,
332  typename MatrixSize,
333  typename InputStorage,
334  typename ... InputP,
335  typename OutputStorage,
336  typename ... OutputP>
337 class Multiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
338  MatrixOrdinal,
339  MatrixDevice,
340  MatrixMemory,
341  MatrixSize >,
342  Kokkos::View< Sacado::MP::Vector<InputStorage>**,
343  InputP... >,
344  Kokkos::View< Sacado::MP::Vector<OutputStorage>**,
345  OutputP... >
346  >
347 {
348 public:
352 
353  typedef MatrixDevice execution_space;
354  typedef typename execution_space::size_type size_type;
355 
356  typedef KokkosSparse::CrsMatrix< MatrixValue,
357  MatrixOrdinal,
358  MatrixDevice,
359  MatrixMemory,
360  MatrixSize > matrix_type;
361  typedef typename matrix_type::values_type matrix_values_type;
362  typedef Kokkos::View< InputVectorValue**,
363  InputP... > input_vector_type;
364  typedef Kokkos::View< OutputVectorValue**,
365  OutputP... > output_vector_type;
366 
367 public:
368 
369  static void apply( const matrix_type & A,
370  const input_vector_type & x,
371  const output_vector_type & y )
372  {
374  multiply_type::apply(A,x,y,details::MultiplyAssign());
375  }
376 };
377 
378 } // namespace Stokhos
379 
380 namespace KokkosSparse {
381 
382 template <typename AlphaType,
383  typename BetaType,
384  typename MatrixType,
385  typename InputType,
386  typename ... InputP,
387  typename OutputType,
388  typename ... OutputP>
389 typename std::enable_if<
390  Kokkos::is_view_mp_vector< Kokkos::View< InputType, InputP... > >::value &&
391  Kokkos::is_view_mp_vector< Kokkos::View< OutputType, OutputP... > >::value
392  >::type
394  const char mode[],
395  const AlphaType& a,
396  const MatrixType& A,
397  const Kokkos::View< InputType, InputP... >& x,
398  const BetaType& b,
399  const Kokkos::View< OutputType, OutputP... >& y,
400  const RANK_ONE)
401 {
402  typedef Kokkos::View< OutputType, OutputP... > OutputVectorType;
403  typedef Kokkos::View< InputType, InputP... > InputVectorType;
404  typedef typename InputVectorType::array_type::non_const_value_type value_type;
405 
406  if(mode[0]!='N') {
408  "Stokhos spmv not implemented for transposed or conjugated matrix-vector multiplies");
409  }
410 
411  if (!Sacado::is_constant(a) || !Sacado::is_constant(b)) {
413  "MV_Multiply not implemented for non-constant a or b");
414  }
415 
416  value_type aa = Sacado::Value<AlphaType>::eval(a);
417  value_type bb = Sacado::Value<BetaType>::eval(b);
418  if (bb == value_type(0)) {
419  if (aa == value_type(1)) {
420  // y = A*x
421  typedef Stokhos::details::MultiplyAssign UpdateType;
422  typedef Stokhos::details::MPMultiply<MatrixType,
423  InputVectorType,OutputVectorType,UpdateType> multiply_type;
424  multiply_type::apply( A, x, y, UpdateType() );
425  }
426  else {
427  // y = a*A*x
429  typedef Stokhos::details::MPMultiply<MatrixType,
430  InputVectorType,OutputVectorType,UpdateType> multiply_type;
431  multiply_type::apply( A, x, y, UpdateType(aa) );
432  }
433  }
434  else if (bb == value_type(1)) {
435  if (aa == value_type(1)) {
436  // y += A*x
437  typedef Stokhos::details::MultiplyUpdate UpdateType;
438  typedef Stokhos::details::MPMultiply<MatrixType,
439  InputVectorType,OutputVectorType,UpdateType> multiply_type;
440  multiply_type::apply( A, x, y, UpdateType() );
441  }
442  else {
443  // y += a*A*x
445  typedef Stokhos::details::MPMultiply<MatrixType,
446  InputVectorType,OutputVectorType,UpdateType> multiply_type;
447  multiply_type::apply( A, x, y, UpdateType(aa) );
448  }
449  }
450  else {
451  // y = a*A*x + b*y
453  typedef Stokhos::details::MPMultiply<MatrixType,
454  InputVectorType,OutputVectorType,UpdateType> multiply_type;
455  multiply_type::apply( A, x, y, UpdateType(aa,bb) );
456  }
457 }
458 
459 template <typename AlphaType,
460  typename BetaType,
461  typename MatrixType,
462  typename InputType,
463  typename ... InputP,
464  typename OutputType,
465  typename ... OutputP>
466 typename std::enable_if<
467  Kokkos::is_view_mp_vector< Kokkos::View< InputType, InputP... > >::value &&
468  Kokkos::is_view_mp_vector< Kokkos::View< OutputType, OutputP... > >::value
469  >::type
471  const char mode[],
472  const AlphaType& a,
473  const MatrixType& A,
474  const Kokkos::View< InputType, InputP... >& x,
475  const BetaType& b,
476  const Kokkos::View< OutputType, OutputP... >& y,
477  const RANK_TWO)
478 {
479  if(mode[0]!='N') {
481  "Stokhos spmv not implemented for transposed or conjugated matrix-vector multiplies");
482  }
483  if (y.dimension_1() == 1) {
484  auto y_1D = subview(y, Kokkos::ALL(), 0);
485  auto x_1D = subview(x, Kokkos::ALL(), 0);
486  spmv(mode, a, A, x_1D, b, y_1D, RANK_ONE());
487  }
488  else {
489  typedef Kokkos::View< OutputType, OutputP... > OutputVectorType;
490  typedef Kokkos::View< InputType, InputP... > InputVectorType;
491  typedef typename InputVectorType::array_type::non_const_value_type value_type;
492 
493  if (!Sacado::is_constant(a) || !Sacado::is_constant(b)) {
495  "Stokhos spmv not implemented for non-constant a or b");
496  }
497 
498  value_type aa = Sacado::Value<AlphaType>::eval(a);
499  value_type bb = Sacado::Value<BetaType>::eval(b);
500  if (bb == value_type(0)) {
501  if (aa == value_type(1)) {
502  // y = A*x
503  typedef Stokhos::details::MultiplyAssign UpdateType;
504  typedef Stokhos::details::MPMultiply<MatrixType,
505  InputVectorType,OutputVectorType,UpdateType> multiply_type;
506  multiply_type::apply( A, x, y, UpdateType() );
507  }
508  else {
509  // y = a*A*x
511  typedef Stokhos::details::MPMultiply<MatrixType,
512  InputVectorType,OutputVectorType,UpdateType> multiply_type;
513  multiply_type::apply( A, x, y, UpdateType(aa) );
514  }
515  }
516  else if (bb == value_type(1)) {
517  if (aa == value_type(1)) {
518  // y += A*x
519  typedef Stokhos::details::MultiplyUpdate UpdateType;
520  typedef Stokhos::details::MPMultiply<MatrixType,
521  InputVectorType,OutputVectorType,UpdateType> multiply_type;
522  multiply_type::apply( A, x, y, UpdateType() );
523  }
524  else {
525  // y += a*A*x
527  typedef Stokhos::details::MPMultiply<MatrixType,
528  InputVectorType,OutputVectorType,UpdateType> multiply_type;
529  multiply_type::apply( A, x, y, UpdateType(aa) );
530  }
531  }
532  else {
533  // y = a*A*x + b*y
535  typedef Stokhos::details::MPMultiply<MatrixType,
536  InputVectorType,OutputVectorType,UpdateType> multiply_type;
537  multiply_type::apply( A, x, y, UpdateType(aa,bb) );
538  }
539  }
540 }
541 
542 }
543 
544 #endif /* #ifndef KOKKOS_CRSMATRIX_MP_VECTOR_HPP */
KOKKOS_INLINE_FUNCTION void raise_error(const char *msg)
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< RD, RP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< XD, XP... > >::value >::type sum(const Kokkos::View< RD, RP... > &r, const Kokkos::View< XD, XP... > &x)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
Top-level namespace for Stokhos classes and functions.
KOKKOS_INLINE_FUNCTION bool is_constant(const T &x)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267
void update(const ValueType &alpha, VectorType &x, const ValueType &beta, const VectorType &y)
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< InputType, InputP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< OutputType, OutputP... > >::value >::type spmv(const char mode[], const AlphaType &a, const MatrixType &A, const Kokkos::View< InputType, InputP... > &x, const BetaType &b, const Kokkos::View< OutputType, OutputP... > &y, const RANK_ONE)