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