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"
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 
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< 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< 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<MatrixDevice,Kokkos::Cuda>::value >::type
124 #endif
125  >
126 {
127 
128 public:
129 
134 
135  typedef MatrixDevice execution_space;
136  typedef typename execution_space::size_type size_type;
137 
138  typedef KokkosSparse::CrsMatrix< 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;
148 
153 
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< 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<MatrixDevice,Kokkos::Cuda>::value >::type
218 #endif
219  >
220 {
221 public:
226 
227  typedef MatrixDevice execution_space;
228  typedef typename execution_space::size_type size_type;
229 
230  typedef KokkosSparse::CrsMatrix< 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;
241 
246 
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  {
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 } // namespace details
289 
290 // Kernel implementing y = A * x where
291 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
292 // x, y == Kokkos::View< Sacado::MP::Vector<...>*,...>,
293 // x and y are rank 1, any layout
294 // We spell everything out here to make sure the ranks and devices match.
295 //
296 // This implementation uses overloaded operators for MP::Vector.
297 template <typename MatrixDevice,
298  typename MatrixStorage,
299  typename MatrixOrdinal,
300  typename MatrixMemory,
301  typename MatrixSize,
302  typename InputStorage,
303  typename ... InputP,
304  typename OutputStorage,
305  typename ... OutputP>
306 class Multiply< 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  >
316 {
317 public:
321 
322  typedef MatrixDevice execution_space;
323  typedef typename execution_space::size_type size_type;
324 
325  typedef KokkosSparse::CrsMatrix< MatrixValue,
326  MatrixOrdinal,
327  MatrixDevice,
328  MatrixMemory,
329  MatrixSize > matrix_type;
330  typedef typename matrix_type::values_type matrix_values_type;
331  typedef Kokkos::View< const InputVectorValue*,
332  InputP... > input_vector_type;
333  typedef Kokkos::View< OutputVectorValue*,
334  OutputP... > output_vector_type;
335 
336 public:
337 
338  static void apply( const matrix_type & A,
339  const input_vector_type & x,
340  const output_vector_type & y )
341  {
343  multiply_type::apply(A,x,y,details::MultiplyAssign());
344  }
345 };
346 
347 // Kernel implementing y = A * x where
348 // A == KokkosSparse::CrsMatrix< Sacado::MP::Vector<...>,...>,
349 // x, y == Kokkos::View< Sacado::MP::Vector<...>**,...>,
350 // x and y are rank 2, any layout
351 // We spell everything out here to make sure the ranks and devices match.
352 //
353 // This implementation uses overloaded operators for MP::Vector.
354 template <typename MatrixDevice,
355  typename MatrixStorage,
356  typename MatrixOrdinal,
357  typename MatrixMemory,
358  typename MatrixSize,
359  typename InputStorage,
360  typename ... InputP,
361  typename OutputStorage,
362  typename ... OutputP>
363 class Multiply< KokkosSparse::CrsMatrix< Sacado::MP::Vector<MatrixStorage>,
364  MatrixOrdinal,
365  MatrixDevice,
366  MatrixMemory,
367  MatrixSize >,
368  Kokkos::View< const Sacado::MP::Vector<InputStorage>**,
369  InputP... >,
370  Kokkos::View< Sacado::MP::Vector<OutputStorage>**,
371  OutputP... >
372  >
373 {
374 public:
378 
379  typedef MatrixDevice execution_space;
380  typedef typename execution_space::size_type size_type;
381 
382  typedef KokkosSparse::CrsMatrix< MatrixValue,
383  MatrixOrdinal,
384  MatrixDevice,
385  MatrixMemory,
386  MatrixSize > matrix_type;
387  typedef typename matrix_type::values_type matrix_values_type;
388  typedef Kokkos::View< const InputVectorValue**,
389  InputP... > input_vector_type;
390  typedef Kokkos::View< OutputVectorValue**,
391  OutputP... > output_vector_type;
392 
393 public:
394 
395  static void apply( const matrix_type & A,
396  const input_vector_type & x,
397  const output_vector_type & y )
398  {
400  multiply_type::apply(A,x,y,details::MultiplyAssign());
401  }
402 };
403 
404 } // namespace Stokhos
405 
406 namespace KokkosSparse {
407 
408 template <typename AlphaType,
409  typename BetaType,
410  typename MatrixType,
411  typename InputType,
412  typename ... InputP,
413  typename OutputType,
414  typename ... OutputP>
415 typename std::enable_if<
416  Kokkos::is_view_mp_vector< Kokkos::View< InputType, InputP... > >::value &&
417  Kokkos::is_view_mp_vector< Kokkos::View< OutputType, OutputP... > >::value
418  >::type
420  const char mode[],
421  const AlphaType& a,
422  const MatrixType& A,
423  const Kokkos::View< InputType, InputP... >& x,
424  const BetaType& b,
425  const Kokkos::View< OutputType, OutputP... >& y,
426  const RANK_ONE)
427 {
428  typedef Kokkos::View< OutputType, OutputP... > OutputVectorType;
429  typedef Kokkos::View< InputType, InputP... > InputVectorType;
430  using input_vector_type = const_type_t<InputVectorType>;
431  typedef typename InputVectorType::array_type::non_const_value_type value_type;
432 
433  if(mode[0]!='N') {
435  "Stokhos spmv not implemented for transposed or conjugated matrix-vector multiplies");
436  }
437 
438  if (!Sacado::is_constant(a) || !Sacado::is_constant(b)) {
440  "MV_Multiply not implemented for non-constant a or b");
441  }
442 
443  value_type aa = Sacado::Value<AlphaType>::eval(a);
444  value_type bb = Sacado::Value<BetaType>::eval(b);
445  if (bb == value_type(0)) {
446  if (aa == value_type(1)) {
447  // y = A*x
448  typedef Stokhos::details::MultiplyAssign UpdateType;
449  typedef Stokhos::details::MPMultiply<MatrixType,
450  input_vector_type, OutputVectorType,
451  UpdateType> multiply_type;
452  multiply_type::apply( A, x, y, UpdateType() );
453  }
454  else {
455  // y = a*A*x
457  typedef Stokhos::details::MPMultiply<MatrixType,
458  input_vector_type, OutputVectorType,
459  UpdateType> multiply_type;
460  multiply_type::apply( A, x, y, UpdateType(aa) );
461  }
462  }
463  else if (bb == value_type(1)) {
464  if (aa == value_type(1)) {
465  // y += A*x
466  typedef Stokhos::details::MultiplyUpdate UpdateType;
467  typedef Stokhos::details::MPMultiply<MatrixType,
468  input_vector_type, OutputVectorType,
469  UpdateType> multiply_type;
470  multiply_type::apply( A, x, y, UpdateType() );
471  }
472  else {
473  // y += a*A*x
475  typedef Stokhos::details::MPMultiply<MatrixType,
476  input_vector_type, OutputVectorType,
477  UpdateType> multiply_type;
478  multiply_type::apply( A, x, y, UpdateType(aa) );
479  }
480  }
481  else {
482  // y = a*A*x + b*y
484  typedef Stokhos::details::MPMultiply<MatrixType,
485  input_vector_type, OutputVectorType,
486  UpdateType> multiply_type;
487  multiply_type::apply( A, x, y, UpdateType(aa,bb) );
488  }
489 }
490 
491 template <typename AlphaType,
492  typename BetaType,
493  typename MatrixType,
494  typename InputType,
495  typename ... InputP,
496  typename OutputType,
497  typename ... OutputP>
498 typename std::enable_if<
499  Kokkos::is_view_mp_vector< Kokkos::View< InputType, InputP... > >::value &&
500  Kokkos::is_view_mp_vector< Kokkos::View< OutputType, OutputP... > >::value
501  >::type
503  const char mode[],
504  const AlphaType& a,
505  const MatrixType& A,
506  const Kokkos::View< InputType, InputP... >& x,
507  const BetaType& b,
508  const Kokkos::View< OutputType, OutputP... >& y,
509  const RANK_TWO)
510 {
511  if(mode[0]!='N') {
513  "Stokhos spmv not implemented for transposed or conjugated matrix-vector multiplies");
514  }
515  if (y.extent(1) == 1) {
516  auto y_1D = subview(y, Kokkos::ALL(), 0);
517  auto x_1D = subview(x, Kokkos::ALL(), 0);
518  spmv(mode, a, A, x_1D, b, y_1D, RANK_ONE());
519  }
520  else {
521  typedef Kokkos::View< OutputType, OutputP... > OutputVectorType;
522  typedef Kokkos::View< InputType, InputP... > InputVectorType;
523  using input_vector_type = const_type_t<InputVectorType>;
524  typedef typename InputVectorType::array_type::non_const_value_type value_type;
525 
526  if (!Sacado::is_constant(a) || !Sacado::is_constant(b)) {
528  "Stokhos spmv not implemented for non-constant a or b");
529  }
530 
531  value_type aa = Sacado::Value<AlphaType>::eval(a);
532  value_type bb = Sacado::Value<BetaType>::eval(b);
533  if (bb == value_type(0)) {
534  if (aa == value_type(1)) {
535  // y = A*x
536  typedef Stokhos::details::MultiplyAssign UpdateType;
537  typedef Stokhos::details::MPMultiply<MatrixType,
538  input_vector_type, OutputVectorType,
539  UpdateType> multiply_type;
540  multiply_type::apply( A, x, y, UpdateType() );
541  }
542  else {
543  // y = a*A*x
545  typedef Stokhos::details::MPMultiply<MatrixType,
546  input_vector_type, OutputVectorType,
547  UpdateType> multiply_type;
548  multiply_type::apply( A, x, y, UpdateType(aa) );
549  }
550  }
551  else if (bb == value_type(1)) {
552  if (aa == value_type(1)) {
553  // y += A*x
554  typedef Stokhos::details::MultiplyUpdate UpdateType;
555  typedef Stokhos::details::MPMultiply<MatrixType,
556  input_vector_type, OutputVectorType,
557  UpdateType> multiply_type;
558  multiply_type::apply( A, x, y, UpdateType() );
559  }
560  else {
561  // y += a*A*x
563  typedef Stokhos::details::MPMultiply<MatrixType,
564  input_vector_type, OutputVectorType,
565  UpdateType> multiply_type;
566  multiply_type::apply( A, x, y, UpdateType(aa) );
567  }
568  }
569  else {
570  // y = a*A*x + b*y
572  typedef Stokhos::details::MPMultiply<MatrixType,
573  input_vector_type, OutputVectorType,
574  UpdateType> multiply_type;
575  multiply_type::apply( A, x, y, UpdateType(aa,bb) );
576  }
577  }
578 }
579 
580 }
581 
582 #endif /* #ifndef KOKKOS_CRSMATRIX_MP_VECTOR_HPP */
KOKKOS_INLINE_FUNCTION void raise_error(const char *msg)
KOKKOS_INLINE_FUNCTION bool is_constant(const T &x)
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)
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)