Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Kokkos_CrsMatrix_MP_Vector.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef KOKKOS_CRSMATRIX_MP_VECTOR_HPP
11 #define KOKKOS_CRSMATRIX_MP_VECTOR_HPP
12 
13 #include <type_traits>
14 
15 #include "Sacado_MP_Vector.hpp"
16 #include "Kokkos_View_MP_Vector.hpp"
17 #include "KokkosSparse_CrsMatrix.hpp"
18 #include "Kokkos_ArithTraits_MP_Vector.hpp"
19 #include "KokkosSparse_spmv.hpp"
20 #include "Kokkos_Blas1_MP_Vector.hpp" // for some utilities
21 
22 #include "Kokkos_Core.hpp"
23 #include "Stokhos_Multiply.hpp"
24 
25 #include "Teuchos_TestForException.hpp"
26 
27 //----------------------------------------------------------------------------
28 // Specializations of KokkosSparse::CrsMatrix for Sacado::MP::Vector scalar type
29 //----------------------------------------------------------------------------
30 
31 namespace { // (anonymous)
32 
33 template<class T, class = std::void_t<> >
34 struct const_type_impl {
35  using type = T;
36 };
37 
38 template<class T>
39 struct const_type_impl<T,
40  std::void_t<typename T::const_type> > {
41  using type = typename T::const_type;
42 };
43 
44 template<class T>
45 using const_type_t = typename const_type_impl<T>::type;
46 
47 } // namespace (anonymous)
48 
49 namespace Stokhos {
50 
51 namespace details {
52 
53 template <typename Matrix, typename InputVector, typename OutputVector,
54  typename Update = MultiplyAssign,
55  typename Enabled = void>
56 class MPMultiply {};
57 
58 // Kernel implementing y = A * x where
59 // A == KokkosSparse::CrsMatrix< const Sacado::MP::Vector<...>,...>,
60 // x, y == Kokkos::View< Sacado::MP::Vector<...>*,...>,
61 // x and y are rank 1, any layout
62 // We spell everything out here to make sure the ranks and devices match.
63 //
64 // This implementation uses overloaded operators for MP::Vector.
65 template <typename MatrixDevice,
66  typename MatrixStorage,
67  typename MatrixOrdinal,
68  typename MatrixMemory,
69  typename MatrixSize,
70  typename InputStorage,
71  typename ... InputP,
72  typename OutputStorage,
73  typename ... OutputP,
74  typename Update>
75 class MPMultiply< KokkosSparse::CrsMatrix< const Sacado::MP::Vector<MatrixStorage>,
76  MatrixOrdinal,
77  MatrixDevice,
78  MatrixMemory,
79  MatrixSize>,
80  Kokkos::View< const Sacado::MP::Vector<InputStorage>*,
81  InputP... >,
82  Kokkos::View< Sacado::MP::Vector<OutputStorage>*,
83  OutputP... >,
84  Update
85 #ifdef KOKKOS_ENABLE_CUDA
86  , typename std::enable_if<
87  !std::is_same<typename MatrixDevice::execution_space,Kokkos::Cuda>::value >::type
88 #endif
89  >
90 {
91 
92 public:
93 
94  typedef Sacado::MP::Vector<MatrixStorage> MatrixValue;
95  typedef Sacado::MP::Vector<InputStorage> InputVectorValue;
96  typedef Sacado::MP::Vector<OutputStorage> OutputVectorValue;
97  typedef OutputVectorValue scalar_type;
98 
99  typedef typename MatrixDevice::execution_space execution_space;
100  typedef typename execution_space::size_type size_type;
101 
102  typedef KokkosSparse::CrsMatrix< const MatrixValue,
103  MatrixOrdinal,
104  MatrixDevice,
105  MatrixMemory,
106  MatrixSize > matrix_type;
107  typedef Kokkos::View< const InputVectorValue*,
108  InputP... > input_vector_type;
109  typedef Kokkos::View< OutputVectorValue*,
110  OutputP... > output_vector_type;
111  typedef Update update_type;
112 
113  const matrix_type m_A;
114  const input_vector_type m_x;
115  const output_vector_type m_y;
116  const update_type m_update;
117 
118  MPMultiply( const matrix_type & A,
119  const input_vector_type & x,
120  const output_vector_type & y,
121  const update_type& update )
122  : m_A( A )
123  , m_x( x )
124  , m_y( y )
125  , m_update( update )
126  {}
127 
128  KOKKOS_INLINE_FUNCTION
129  void operator()( const size_type iRow ) const
130  {
131  // Compute mat-vec for this row
132  const size_type iEntryBegin = m_A.graph.row_map[iRow];
133  const size_type iEntryEnd = m_A.graph.row_map[iRow+1];
134  scalar_type sum = 0.0;
135  for (size_type iEntry = iEntryBegin; iEntry < iEntryEnd; ++iEntry) {
136  size_type iCol = m_A.graph.entries(iEntry);
137  sum += m_A.values(iEntry) * m_x(iCol);
138  }
139  m_update( m_y(iRow), sum );
140  } // operator()
141 
142  static void apply( const matrix_type & A,
143  const input_vector_type & x,
144  const output_vector_type & y,
145  const update_type & update )
146  {
147  const size_type row_count = A.graph.row_map.extent(0)-1;
148  Kokkos::parallel_for( row_count, MPMultiply(A,x,y,update) );
149  }
150 };
151 
152 // Kernel implementing y = A * x where
153 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
154 // x, y == Kokkos::View< Sacado::MP::Vector<...>**,...>,
155 // x and y are rank 2, any layout
156 // We spell everything out here to make sure the ranks and devices match.
157 //
158 // This implementation uses overloaded operators for MP::Vector.
159 template <typename MatrixDevice,
160  typename MatrixStorage,
161  typename MatrixOrdinal,
162  typename MatrixMemory,
163  typename MatrixSize,
164  typename InputStorage,
165  typename ... InputP,
166  typename OutputStorage,
167  typename ... OutputP,
168  typename Update>
169 class MPMultiply< KokkosSparse::CrsMatrix< const Sacado::MP::Vector<MatrixStorage>,
170  MatrixOrdinal,
171  MatrixDevice,
172  MatrixMemory,
173  MatrixSize >,
174  Kokkos::View< const Sacado::MP::Vector<InputStorage>**,
175  InputP... >,
176  Kokkos::View< Sacado::MP::Vector<OutputStorage>**,
177  OutputP... >,
178  Update
179 #ifdef KOKKOS_ENABLE_CUDA
180  , typename std::enable_if<
181  !std::is_same<typename MatrixDevice::execution_space,Kokkos::Cuda>::value >::type
182 #endif
183  >
184 {
185 public:
186  typedef Sacado::MP::Vector<MatrixStorage> MatrixValue;
187  typedef Sacado::MP::Vector<InputStorage> InputVectorValue;
188  typedef Sacado::MP::Vector<OutputStorage> OutputVectorValue;
189  typedef OutputVectorValue scalar_type;
190 
191  typedef typename MatrixDevice::execution_space execution_space;
192  typedef typename execution_space::size_type size_type;
193 
194  typedef KokkosSparse::CrsMatrix< const MatrixValue,
195  MatrixOrdinal,
196  MatrixDevice,
197  MatrixMemory,
198  MatrixSize > matrix_type;
199  typedef typename matrix_type::values_type matrix_values_type;
200  typedef Kokkos::View< const InputVectorValue**,
201  InputP... > input_vector_type;
202  typedef Kokkos::View< OutputVectorValue**,
203  OutputP... > output_vector_type;
204  typedef Update update_type;
205 
206  const matrix_type m_A;
207  const input_vector_type m_x;
208  const output_vector_type m_y;
209  const update_type m_update;
210 
211  MPMultiply( const matrix_type & A,
212  const input_vector_type & x,
213  const output_vector_type & y,
214  const update_type& update )
215  : m_A( A )
216  , m_x( x )
217  , m_y( y )
218  , m_update( update )
219  {}
220 
221  KOKKOS_INLINE_FUNCTION
222  void operator()( const size_type iRow ) const
223  {
224  scalar_type sum;
225  // Loop over columns of x, y
226  const size_type num_col = m_y.extent(1);
227  for (size_type col=0; col<num_col; ++col) {
228  // Compute mat-vec for this row
229  const size_type iEntryBegin = m_A.graph.row_map[iRow];
230  const size_type iEntryEnd = m_A.graph.row_map[iRow+1];
231  sum = 0.0;
232  for (size_type iEntry = iEntryBegin; iEntry < iEntryEnd; ++iEntry) {
233  size_type iCol = m_A.graph.entries(iEntry);
234  sum += m_A.values(iEntry) * m_x(iCol,col);
235  }
236  m_update( m_y(iRow,col), sum );
237  } // x, y column loop
238  } // operator()
239 
240 public:
241 
242  static void apply( const matrix_type & A,
243  const input_vector_type & x,
244  const output_vector_type & y,
245  const update_type & update )
246  {
247  const size_type row_count = A.graph.row_map.extent(0)-1;
248  Kokkos::parallel_for( row_count, MPMultiply(A,x,y,update) );
249  }
250 };
251 
252 // Kernel implementing y = A * x where
253 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
254 // x, y == Kokkos::View< Sacado::MP::Vector<...>*,...>,
255 // x and y are rank 1, any layout
256 // We spell everything out here to make sure the ranks and devices match.
257 //
258 // This implementation uses overloaded operators for MP::Vector.
259 template <typename MatrixDevice,
260  typename MatrixStorage,
261  typename MatrixOrdinal,
262  typename MatrixMemory,
263  typename MatrixSize,
264  typename InputStorage,
265  typename ... InputP,
266  typename OutputStorage,
267  typename ... OutputP,
268  typename Update>
269 class MPMultiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
270  MatrixOrdinal,
271  MatrixDevice,
272  MatrixMemory,
273  MatrixSize>,
274  Kokkos::View< const Sacado::MP::Vector<InputStorage>*,
275  InputP... >,
276  Kokkos::View< Sacado::MP::Vector<OutputStorage>*,
277  OutputP... >,
278  Update
279  >
280 {
281 
282 public:
283 
284  typedef Sacado::MP::Vector<MatrixStorage> MatrixValue;
285  typedef Sacado::MP::Vector<InputStorage> InputVectorValue;
286  typedef Sacado::MP::Vector<OutputStorage> OutputVectorValue;
287  typedef OutputVectorValue scalar_type;
288 
289  typedef typename MatrixDevice::execution_space execution_space;
290  typedef typename execution_space::size_type size_type;
291 
292  typedef KokkosSparse::CrsMatrix< MatrixValue,
293  MatrixOrdinal,
294  MatrixDevice,
295  MatrixMemory,
296  MatrixSize > matrix_type;
297  typedef typename matrix_type::const_type const_matrix_type;
298 
299  typedef Kokkos::View< const InputVectorValue*,
300  InputP... > input_vector_type;
301  typedef Kokkos::View< OutputVectorValue*,
302  OutputP... > output_vector_type;
303  typedef Update update_type;
304 
305  static void apply( const matrix_type & A,
306  const input_vector_type & x,
307  const output_vector_type & y,
308  const update_type & update )
309  {
310  const_matrix_type cA = A;
311  MPMultiply< const_matrix_type, input_vector_type, output_vector_type, update_type >::apply(
312  cA, x, y, update);
313  }
314 };
315 
316 // Kernel implementing y = A * x where
317 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
318 // x, y == Kokkos::View< Sacado::MP::Vector<...>**,...>,
319 // x and y are rank 2, any layout
320 // We spell everything out here to make sure the ranks and devices match.
321 //
322 // This implementation uses overloaded operators for MP::Vector.
323 template <typename MatrixDevice,
324  typename MatrixStorage,
325  typename MatrixOrdinal,
326  typename MatrixMemory,
327  typename MatrixSize,
328  typename InputStorage,
329  typename ... InputP,
330  typename OutputStorage,
331  typename ... OutputP,
332  typename Update>
333 class MPMultiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
334  MatrixOrdinal,
335  MatrixDevice,
336  MatrixMemory,
337  MatrixSize>,
338  Kokkos::View< const Sacado::MP::Vector<InputStorage>**,
339  InputP... >,
340  Kokkos::View< Sacado::MP::Vector<OutputStorage>**,
341  OutputP... >,
342  Update
343  >
344 {
345 
346 public:
347 
348  typedef Sacado::MP::Vector<MatrixStorage> MatrixValue;
349  typedef Sacado::MP::Vector<InputStorage> InputVectorValue;
350  typedef Sacado::MP::Vector<OutputStorage> OutputVectorValue;
351  typedef OutputVectorValue scalar_type;
352 
353  typedef typename MatrixDevice::execution_space 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::const_type const_matrix_type;
362 
363  typedef Kokkos::View< const InputVectorValue**,
364  InputP... > input_vector_type;
365  typedef Kokkos::View< OutputVectorValue**,
366  OutputP... > output_vector_type;
367  typedef Update update_type;
368 
369  static void apply( const matrix_type & A,
370  const input_vector_type & x,
371  const output_vector_type & y,
372  const update_type & update )
373  {
374  const_matrix_type cA = A;
375  MPMultiply< const_matrix_type, input_vector_type, output_vector_type, update_type >::apply(
376  cA, x, y, update);
377  }
378 };
379 
380 } // namespace details
381 
382 // Kernel implementing y = A * x where
383 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
384 // x, y == Kokkos::View< Sacado::MP::Vector<...>*,...>,
385 // x and y are rank 1, any layout
386 // We spell everything out here to make sure the ranks and devices match.
387 //
388 // This implementation uses overloaded operators for MP::Vector.
389 template <typename MatrixDevice,
390  typename MatrixStorage,
391  typename MatrixOrdinal,
392  typename MatrixMemory,
393  typename MatrixSize,
394  typename InputStorage,
395  typename ... InputP,
396  typename OutputStorage,
397  typename ... OutputP>
398 class Multiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
399  MatrixOrdinal,
400  MatrixDevice,
401  MatrixMemory,
402  MatrixSize >,
403  Kokkos::View< const Sacado::MP::Vector<InputStorage>*,
404  InputP... >,
405  Kokkos::View< Sacado::MP::Vector<OutputStorage>*,
406  OutputP... >
407  >
408 {
409 public:
410  typedef Sacado::MP::Vector<MatrixStorage> MatrixValue;
411  typedef Sacado::MP::Vector<InputStorage> InputVectorValue;
412  typedef Sacado::MP::Vector<OutputStorage> OutputVectorValue;
413 
414  typedef typename MatrixDevice::execution_space execution_space;
415  typedef typename execution_space::size_type size_type;
416 
417  typedef KokkosSparse::CrsMatrix< MatrixValue,
418  MatrixOrdinal,
419  MatrixDevice,
420  MatrixMemory,
421  MatrixSize > matrix_type;
422  typedef typename matrix_type::values_type matrix_values_type;
423  typedef Kokkos::View< const InputVectorValue*,
424  InputP... > input_vector_type;
425  typedef Kokkos::View< OutputVectorValue*,
426  OutputP... > output_vector_type;
427 
428 public:
429 
430  static void apply( const matrix_type & A,
431  const input_vector_type & x,
432  const output_vector_type & y )
433  {
434  typedef details::MPMultiply<matrix_type,input_vector_type,output_vector_type> multiply_type;
435  multiply_type::apply(A,x,y,details::MultiplyAssign());
436  }
437 };
438 
439 // Kernel implementing y = A * x where
440 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
441 // x, y == Kokkos::View< Sacado::MP::Vector<...>**,...>,
442 // x and y are rank 2, any layout
443 // We spell everything out here to make sure the ranks and devices match.
444 //
445 // This implementation uses overloaded operators for MP::Vector.
446 template <typename MatrixDevice,
447  typename MatrixStorage,
448  typename MatrixOrdinal,
449  typename MatrixMemory,
450  typename MatrixSize,
451  typename InputStorage,
452  typename ... InputP,
453  typename OutputStorage,
454  typename ... OutputP>
455 class Multiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
456  MatrixOrdinal,
457  MatrixDevice,
458  MatrixMemory,
459  MatrixSize >,
460  Kokkos::View< const Sacado::MP::Vector<InputStorage>**,
461  InputP... >,
462  Kokkos::View< Sacado::MP::Vector<OutputStorage>**,
463  OutputP... >
464  >
465 {
466 public:
467  typedef Sacado::MP::Vector<MatrixStorage> MatrixValue;
468  typedef Sacado::MP::Vector<InputStorage> InputVectorValue;
469  typedef Sacado::MP::Vector<OutputStorage> OutputVectorValue;
470 
471  typedef typename MatrixDevice::execution_space execution_space;
472  typedef typename execution_space::size_type size_type;
473 
474  typedef KokkosSparse::CrsMatrix< MatrixValue,
475  MatrixOrdinal,
476  MatrixDevice,
477  MatrixMemory,
478  MatrixSize > matrix_type;
479  typedef typename matrix_type::values_type matrix_values_type;
480  typedef Kokkos::View< const InputVectorValue**,
481  InputP... > input_vector_type;
482  typedef Kokkos::View< OutputVectorValue**,
483  OutputP... > output_vector_type;
484 
485 public:
486 
487  static void apply( const matrix_type & A,
488  const input_vector_type & x,
489  const output_vector_type & y )
490  {
491  typedef details::MPMultiply<matrix_type,input_vector_type,output_vector_type> multiply_type;
492  multiply_type::apply(A,x,y,details::MultiplyAssign());
493  }
494 };
495 
496 } // namespace Stokhos
497 
498 namespace KokkosSparse {
499 
500 template <
501 #if KOKKOSKERNELS_VERSION >= 40199
502  typename ExecutionSpace,
503 #endif
504 #if KOKKOSKERNELS_VERSION >= 40299
505  typename Handle,
506 #endif
507  typename AlphaType,
508  typename BetaType,
509  typename MatrixType,
510  typename InputType,
511  typename ... InputP,
512  typename OutputType,
513  typename ... OutputP>
514 typename std::enable_if<
515  Kokkos::is_view_mp_vector< Kokkos::View< InputType, InputP... > >::value &&
516  Kokkos::is_view_mp_vector< Kokkos::View< OutputType, OutputP... > >::value
517 #if KOKKOSKERNELS_VERSION >= 40299
518  && KokkosSparse::is_crs_matrix_v<MatrixType>
519  && (Kokkos::View< OutputType, OutputP... >::rank() == 1)
520 #endif
521  >::type
522 spmv(
523 #if KOKKOSKERNELS_VERSION >= 40199
524  const ExecutionSpace& space,
525 #endif
526 #if KOKKOSKERNELS_VERSION < 40299
527  KokkosKernels::Experimental::Controls,
528 #else
529  Handle* handle,
530 #endif
531  const char mode[],
532  const AlphaType& a,
533  const MatrixType& A,
534  const Kokkos::View< InputType, InputP... >& x,
535  const BetaType& b,
536  const Kokkos::View< OutputType, OutputP... >& y
537 #if KOKKOSKERNELS_VERSION < 40299
538  , const RANK_ONE
539 #endif
540 )
541 {
542  typedef Kokkos::View< OutputType, OutputP... > OutputVectorType;
543  typedef Kokkos::View< InputType, InputP... > InputVectorType;
544  using input_vector_type = const_type_t<InputVectorType>;
545  typedef typename InputVectorType::array_type::non_const_value_type value_type;
546 
547 #if KOKKOSKERNELS_VERSION >= 40199
548  if(space != ExecutionSpace()) {
549  Kokkos::Impl::raise_error(
550  "Stokhos spmv not implemented for non-default execution space instance");
551  }
552 #endif
553  if(mode[0]!='N') {
554  Kokkos::Impl::raise_error(
555  "Stokhos spmv not implemented for transposed or conjugated matrix-vector multiplies");
556  }
557 
558  if (!Sacado::is_constant(a) || !Sacado::is_constant(b)) {
559  Kokkos::Impl::raise_error(
560  "MV_Multiply not implemented for non-constant a or b");
561  }
562 
563  value_type aa = Sacado::Value<AlphaType>::eval(a);
564  value_type bb = Sacado::Value<BetaType>::eval(b);
565  if (bb == value_type(0)) {
566  if (aa == value_type(1)) {
567  // y = A*x
568  typedef Stokhos::details::MultiplyAssign UpdateType;
569  typedef Stokhos::details::MPMultiply<MatrixType,
570  input_vector_type, OutputVectorType,
571  UpdateType> multiply_type;
572  multiply_type::apply( A, x, y, UpdateType() );
573  }
574  else {
575  // y = a*A*x
576  typedef Stokhos::details::MultiplyScaledAssign<value_type> UpdateType;
577  typedef Stokhos::details::MPMultiply<MatrixType,
578  input_vector_type, OutputVectorType,
579  UpdateType> multiply_type;
580  multiply_type::apply( A, x, y, UpdateType(aa) );
581  }
582  }
583  else if (bb == value_type(1)) {
584  if (aa == value_type(1)) {
585  // y += A*x
586  typedef Stokhos::details::MultiplyUpdate UpdateType;
587  typedef Stokhos::details::MPMultiply<MatrixType,
588  input_vector_type, OutputVectorType,
589  UpdateType> multiply_type;
590  multiply_type::apply( A, x, y, UpdateType() );
591  }
592  else {
593  // y += a*A*x
594  typedef Stokhos::details::MultiplyScaledUpdate<value_type> UpdateType;
595  typedef Stokhos::details::MPMultiply<MatrixType,
596  input_vector_type, OutputVectorType,
597  UpdateType> multiply_type;
598  multiply_type::apply( A, x, y, UpdateType(aa) );
599  }
600  }
601  else {
602  // y = a*A*x + b*y
603  typedef Stokhos::details::MultiplyScaledUpdate2<value_type> UpdateType;
604  typedef Stokhos::details::MPMultiply<MatrixType,
605  input_vector_type, OutputVectorType,
606  UpdateType> multiply_type;
607  multiply_type::apply( A, x, y, UpdateType(aa,bb) );
608  }
609 }
610 
611 template <
612 #if KOKKOSKERNELS_VERSION >= 40199
613  typename ExecutionSpace,
614 #endif
615 #if KOKKOSKERNELS_VERSION >= 40299
616  typename Handle,
617 #endif
618  typename AlphaType,
619  typename BetaType,
620  typename MatrixType,
621  typename InputType,
622  typename ... InputP,
623  typename OutputType,
624  typename ... OutputP>
625 typename std::enable_if<
626  Kokkos::is_view_mp_vector< Kokkos::View< InputType, InputP... > >::value &&
627  Kokkos::is_view_mp_vector< Kokkos::View< OutputType, OutputP... > >::value
628 #if KOKKOSKERNELS_VERSION >= 40299
629  && KokkosSparse::is_crs_matrix_v<MatrixType>
630  && (Kokkos::View< OutputType, OutputP... >::rank() == 2)
631 #endif
632  >::type
633 spmv(
634 #if KOKKOSKERNELS_VERSION >= 40199
635  const ExecutionSpace& space,
636 #endif
637 #if KOKKOSKERNELS_VERSION < 40299
638  KokkosKernels::Experimental::Controls,
639 #else
640  Handle* handle,
641 #endif
642  const char mode[],
643  const AlphaType& a,
644  const MatrixType& A,
645  const Kokkos::View< InputType, InputP... >& x,
646  const BetaType& b,
647  const Kokkos::View< OutputType, OutputP... >& y
648 #if KOKKOSKERNELS_VERSION < 40299
649  , const RANK_TWO
650 #endif
651  )
652 {
653 #if KOKKOSKERNELS_VERSION >= 40199
654  if(space != ExecutionSpace()) {
655  Kokkos::Impl::raise_error(
656  "Stokhos spmv not implemented for non-default execution space instance");
657  }
658 #endif
659  if(mode[0]!='N') {
660  Kokkos::Impl::raise_error(
661  "Stokhos spmv not implemented for transposed or conjugated matrix-vector multiplies");
662  }
663  if (y.extent(1) == 1) {
664  auto y_1D = subview(y, Kokkos::ALL(), 0);
665  auto x_1D = subview(x, Kokkos::ALL(), 0);
666 #if KOKKOSKERNELS_VERSION >= 40299
667  spmv(space, handle, mode, a, A, x_1D, b, y_1D);
668 #elif (KOKKOSKERNELS_VERSION < 40299) && (KOKKOSKERNELS_VERSION >= 40199)
669  spmv(space, KokkosKernels::Experimental::Controls(), mode, a, A, x_1D, b, y_1D, RANK_ONE());
670 #else
671  spmv(KokkosKernels::Experimental::Controls(), mode, a, A, x_1D, b, y_1D, RANK_ONE());
672 #endif
673  }
674  else {
675  typedef Kokkos::View< OutputType, OutputP... > OutputVectorType;
676  typedef Kokkos::View< InputType, InputP... > InputVectorType;
677  using input_vector_type = const_type_t<InputVectorType>;
678  typedef typename InputVectorType::array_type::non_const_value_type value_type;
679 
680  if (!Sacado::is_constant(a) || !Sacado::is_constant(b)) {
681  Kokkos::Impl::raise_error(
682  "Stokhos spmv not implemented for non-constant a or b");
683  }
684 
685  value_type aa = Sacado::Value<AlphaType>::eval(a);
686  value_type bb = Sacado::Value<BetaType>::eval(b);
687  if (bb == value_type(0)) {
688  if (aa == value_type(1)) {
689  // y = A*x
690  typedef Stokhos::details::MultiplyAssign UpdateType;
691  typedef Stokhos::details::MPMultiply<MatrixType,
692  input_vector_type, OutputVectorType,
693  UpdateType> multiply_type;
694  multiply_type::apply( A, x, y, UpdateType() );
695  }
696  else {
697  // y = a*A*x
698  typedef Stokhos::details::MultiplyScaledAssign<value_type> UpdateType;
699  typedef Stokhos::details::MPMultiply<MatrixType,
700  input_vector_type, OutputVectorType,
701  UpdateType> multiply_type;
702  multiply_type::apply( A, x, y, UpdateType(aa) );
703  }
704  }
705  else if (bb == value_type(1)) {
706  if (aa == value_type(1)) {
707  // y += A*x
708  typedef Stokhos::details::MultiplyUpdate UpdateType;
709  typedef Stokhos::details::MPMultiply<MatrixType,
710  input_vector_type, OutputVectorType,
711  UpdateType> multiply_type;
712  multiply_type::apply( A, x, y, UpdateType() );
713  }
714  else {
715  // y += a*A*x
716  typedef Stokhos::details::MultiplyScaledUpdate<value_type> UpdateType;
717  typedef Stokhos::details::MPMultiply<MatrixType,
718  input_vector_type, OutputVectorType,
719  UpdateType> multiply_type;
720  multiply_type::apply( A, x, y, UpdateType(aa) );
721  }
722  }
723  else {
724  // y = a*A*x + b*y
725  typedef Stokhos::details::MultiplyScaledUpdate2<value_type> UpdateType;
726  typedef Stokhos::details::MPMultiply<MatrixType,
727  input_vector_type, OutputVectorType,
728  UpdateType> multiply_type;
729  multiply_type::apply( A, x, y, UpdateType(aa,bb) );
730  }
731  }
732 }
733 
734 }
735 
736 #endif /* #ifndef KOKKOS_CRSMATRIX_MP_VECTOR_HPP */
Stokhos_Multiply.hpp
Kokkos_Blas1_MP_Vector.hpp
DynamicVecTest::execution_space
Kokkos::DefaultExecutionSpace execution_space
Definition: Stokhos_SacadoMPVectorCommTests.cpp:651
Kokkos::Impl::raise_error
KOKKOS_INLINE_FUNCTION void raise_error(const char *msg)
Definition: Kokkos_View_Utils.hpp:32
Kokkos_View_MP_Vector.hpp
Sacado::is_constant
KOKKOS_INLINE_FUNCTION bool is_constant(const T &x)
Definition: Stokhos_Is_Constant.hpp:21
Kokkos_ArithTraits_MP_Vector.hpp
KokkosSparse::spmv
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(KokkosKernels::Experimental::Controls, 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)
Definition: Kokkos_CrsMatrix_UQ_PCE.hpp:1460
Teuchos_TestForException.hpp
Sacado_MP_Vector.hpp
Kokkos::is_view_mp_vector
Definition: Kokkos_View_MP_Vector_Utils.hpp:72
KokkosBlas::sum
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)
Definition: Kokkos_Blas1_UQ_PCE.hpp:358
Stokhos::update
void update(const ValueType &alpha, VectorType &x, const ValueType &beta, const VectorType &y)
Definition: Stokhos_Update.hpp:55
Generated on Fri Nov 22 2024 09:18:08 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5