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