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