Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Stokhos_KokkosCrsMatrixUQPCEUnitTest.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 #include "Teuchos_UnitTestHelpers.hpp"
43 #include "Stokhos_UnitTestHelpers.hpp"
44 
45 #include "Stokhos_Sacado_Kokkos_UQ_PCE.hpp"
46 #include "KokkosSparse_CrsMatrix.hpp"
47 #include "KokkosSparse_spmv.hpp"
48 #include "Kokkos_CrsMatrix_UQ_PCE.hpp"
49 #include "Kokkos_CrsMatrix_UQ_PCE_Cuda.hpp"
50 #include "Stokhos_LegendreBasis.hpp"
51 #include "Stokhos_CompletePolynomialBasis.hpp"
52 #include "Stokhos_Sparse3Tensor.hpp"
53 
54 // For computing DeviceConfig
55 #include "Kokkos_Core.hpp"
56 
57 // Helper functions
58 template< typename IntType >
59 inline
60 IntType map_fem_graph_coord( const IntType& N,
61  const IntType& i,
62  const IntType& j,
63  const IntType& k )
64 {
65  return k + N * ( j + N * i );
66 }
67 
68 template < typename ordinal >
69 inline
70 ordinal generate_fem_graph( ordinal N,
71  std::vector< std::vector<ordinal> >& graph )
72 {
73  graph.resize( N * N * N, std::vector<ordinal>() );
74 
75  ordinal total = 0;
76 
77  for ( int i = 0; i < (int) N; ++i ) {
78  for ( int j = 0; j < (int) N; ++j ) {
79  for ( int k = 0; k < (int) N; ++k ) {
80 
81  const ordinal row = map_fem_graph_coord((int)N,i,j,k);
82 
83  graph[row].reserve(27);
84 
85  for ( int ii = -1; ii < 2; ++ii ) {
86  for ( int jj = -1; jj < 2; ++jj ) {
87  for ( int kk = -1; kk < 2; ++kk ) {
88  if ( 0 <= i + ii && i + ii < (int) N &&
89  0 <= j + jj && j + jj < (int) N &&
90  0 <= k + kk && k + kk < (int) N ) {
91  ordinal col = map_fem_graph_coord((int)N,i+ii,j+jj,k+kk);
92 
93  graph[row].push_back(col);
94  }
95  }}}
96  total += graph[row].size();
97  }}}
98 
99  return total;
100 }
101 
102 template <typename scalar, typename ordinal>
103 inline
104 scalar generate_matrix_coefficient( const ordinal nFEM,
105  const ordinal nStoch,
106  const ordinal iRowFEM,
107  const ordinal iColFEM,
108  const ordinal iStoch )
109 {
110  const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
111  ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
112 
113  const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
114 
115  return A_fem + A_stoch;
116  //return 1.0;
117 }
118 
119 template <typename scalar, typename ordinal>
120 inline
121 scalar generate_vector_coefficient( const ordinal nFEM,
122  const ordinal nStoch,
123  const ordinal iColFEM,
124  const ordinal iStoch )
125 {
126  const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
127  const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
128  return X_fem + X_stoch;
129  //return 1.0;
130 }
131 
132 template <typename kokkos_cijk_type, typename ordinal_type>
133 kokkos_cijk_type build_cijk(ordinal_type stoch_dim,
134  ordinal_type poly_ord)
135 {
136  using Teuchos::RCP;
137  using Teuchos::rcp;
138  using Teuchos::Array;
139 
140  typedef typename kokkos_cijk_type::value_type value_type;
141  typedef typename kokkos_cijk_type::execution_space execution_space;
142  typedef Stokhos::OneDOrthogPolyBasis<ordinal_type,value_type> one_d_basis;
143  typedef Stokhos::LegendreBasis<ordinal_type,value_type> legendre_basis;
144  typedef Stokhos::CompletePolynomialBasis<ordinal_type,value_type> product_basis;
145  typedef Stokhos::Sparse3Tensor<ordinal_type,value_type> Cijk;
146 
147  // Create product basis
148  Array< RCP<const one_d_basis> > bases(stoch_dim);
149  for (ordinal_type i=0; i<stoch_dim; i++)
150  bases[i] = rcp(new legendre_basis(poly_ord, true));
151  RCP<const product_basis> basis = rcp(new product_basis(bases));
152 
153  // Triple product tensor
154  RCP<Cijk> cijk = basis->computeTripleProductTensor();
155 
156  // Kokkos triple product tensor
157  kokkos_cijk_type kokkos_cijk =
158  Stokhos::create_product_tensor<execution_space>(*basis, *cijk);
159 
160  return kokkos_cijk;
161 }
162 
163 // Reasonable tolerances for common precisions
164 template <typename Scalar> struct ScalarTol {};
165 template <> struct ScalarTol<float> { static float tol() { return 1e-4; } };
166 template <> struct ScalarTol<double> { static double tol() { return 1e-10; } };
167 
168 // Compare two rank-2 views for equality, to given precision
169 template <typename array_type, typename scalar_type>
170 bool compare_rank_2_views(const array_type& y,
171  const array_type& y_exp,
172  const scalar_type rel_tol,
173  const scalar_type abs_tol,
174  Teuchos::FancyOStream& out)
175 {
176  typedef typename array_type::size_type size_type;
177  typename array_type::HostMirror hy = Kokkos::create_mirror_view(y);
178  typename array_type::HostMirror hy_exp = Kokkos::create_mirror_view(y_exp);
179  Kokkos::deep_copy(hy, y);
180  Kokkos::deep_copy(hy_exp, y_exp);
181 
182  size_type num_rows = y.extent(0);
183  size_type num_cols = y.extent(1);
184  bool success = true;
185  for (size_type i=0; i<num_rows; ++i) {
186  for (size_type j=0; j<num_cols; ++j) {
187  scalar_type diff = std::abs( hy(i,j) - hy_exp(i,j) );
188  scalar_type tol = rel_tol*std::abs(hy_exp(i,j)) + abs_tol;
189  bool s = diff < tol;
190  out << "y_expected(" << i << "," << j << ") - "
191  << "y(" << i << "," << j << ") = " << hy_exp(i,j)
192  << " - " << hy(i,j) << " == "
193  << diff << " < " << tol << " : ";
194  if (s)
195  out << "passed";
196  else
197  out << "failed";
198  out << std::endl;
199  success = success && s;
200  }
201  }
202 
203  return success;
204 }
205 
206 template <typename vector_type, typename scalar_type>
207 bool compareRank1(const vector_type& y,
208  const vector_type& y_exp,
209  const scalar_type rel_tol,
210  const scalar_type abs_tol,
211  Teuchos::FancyOStream& out)
212 {
213  typedef typename vector_type::size_type size_type;
214  typename vector_type::HostMirror hy = Kokkos::create_mirror_view(y);
215  typename vector_type::HostMirror hy_exp = Kokkos::create_mirror_view(y_exp);
216  Kokkos::deep_copy(hy, y);
217  Kokkos::deep_copy(hy_exp, y_exp);
218 
219  size_type num_rows = y.extent(0);
220  bool success = true;
221  for (size_type i=0; i<num_rows; ++i) {
222  for (size_type j=0; j<Kokkos::dimension_scalar(y); ++j) {
223  scalar_type diff = std::abs( hy(i).fastAccessCoeff(j) - hy_exp(i).fastAccessCoeff(j) );
224  scalar_type tol = rel_tol*std::abs(hy_exp(i).fastAccessCoeff(j)) + abs_tol;
225  bool s = diff < tol;
226  out << "y_expected(" << i << ").coeff(" << j << ") - "
227  << "y(" << i << ").coeff(" << j << ") = " << hy_exp(i).fastAccessCoeff(j)
228  << " - " << hy(i).fastAccessCoeff(j) << " == "
229  << diff << " < " << tol << " : ";
230  if (s)
231  out << "passed";
232  else
233  out << "failed";
234  out << std::endl;
235  success = success && s;
236  }
237  }
238  return success;
239 }
240 
241 template <typename vector_type, typename scalar_type>
242 bool compareRank2(const vector_type& y,
243  const vector_type& y_exp,
244  const scalar_type rel_tol,
245  const scalar_type abs_tol,
246  Teuchos::FancyOStream& out)
247 {
248  typedef typename vector_type::size_type size_type;
249  typename vector_type::HostMirror hy = Kokkos::create_mirror_view(y);
250  typename vector_type::HostMirror hy_exp = Kokkos::create_mirror_view(y_exp);
251  Kokkos::deep_copy(hy, y);
252  Kokkos::deep_copy(hy_exp, y_exp);
253 
254  size_type num_rows = y.extent(0);
255  size_type num_cols = y.extent(1);
256  bool success = true;
257 
258  for (size_type col = 0; col < num_cols; ++col){
259  for (size_type i=0; i<num_rows; ++i) {
260  for (size_type j=0; j<Kokkos::dimension_scalar(y); ++j) {
261  scalar_type diff = std::abs( hy(i,col).fastAccessCoeff(j) - hy_exp(i,col).fastAccessCoeff(j) );
262  scalar_type tol = rel_tol*std::abs(hy_exp(i,col).fastAccessCoeff(j)) + abs_tol;
263  bool s = diff < tol;
264  out << "y_expected(" << i << ").coeff(" << j << ") - "
265  << "y(" << i << ").coeff(" << j << ") = " << hy_exp(i,col).fastAccessCoeff(j)
266  << " - " << hy(i,col).fastAccessCoeff(j) << " == "
267  << diff << " < " << tol << " : ";
268  if (s)
269  out << "passed";
270  else
271  out << "failed";
272  out << std::endl;
273  success = success && s;
274  }
275  }
276  }
277 
278 
279  return success;
280 }
281 
282 
283 // Helper function to build a diagonal matrix
284 template <typename MatrixType, typename CijkType>
285 MatrixType
286 buildDiagonalMatrix(typename MatrixType::ordinal_type nrow,
287  typename MatrixType::ordinal_type pce_size,
288  const CijkType& cijk) {
289  typedef typename MatrixType::ordinal_type ordinal_type;
290  typedef typename MatrixType::StaticCrsGraphType matrix_graph_type;
291  typedef typename MatrixType::values_type matrix_values_type;
292 
293  std::vector< std::vector<ordinal_type> > graph(nrow);
294  for (ordinal_type i=0; i<nrow; ++i)
295  graph[i] = std::vector<ordinal_type>(1, i);
296  ordinal_type graph_length = nrow;
297 
298  matrix_graph_type matrix_graph =
299  Kokkos::create_staticcrsgraph<matrix_graph_type>("graph", graph);
300  matrix_values_type matrix_values =
301  Kokkos::make_view<matrix_values_type>("values", cijk, graph_length, pce_size);
302 
303  MatrixType matrix("matrix", nrow, matrix_values, matrix_graph);
304  return matrix;
305 }
306 
307 //
308 // Tests
309 //
310 
311 // Kernel to set diagonal of a matrix to prescribed values
312 template <typename MatrixType>
313 struct ReplaceDiagonalValuesKernel {
314  typedef typename MatrixType::execution_space execution_space;
315  typedef typename MatrixType::size_type size_type;
316  typedef typename MatrixType::value_type value_type;
317  typedef typename MatrixType::ordinal_type ordinal_type;
318 
319  const MatrixType m_matrix;
320  ReplaceDiagonalValuesKernel(const MatrixType matrix) : m_matrix(matrix) {};
321 
322  // Replace diagonal entry for row 'i' with a value
323  KOKKOS_INLINE_FUNCTION
324  void operator() (const size_type i) const {
325  const ordinal_type row = i;
326  const ordinal_type col = i;
327  value_type val = value_type(row);
328  m_matrix.replaceValues(row, &col, 1, &val, false, true);
329  }
330 
331  // Kernel launch
332  static void apply(const MatrixType matrix) {
333  const size_type nrow = matrix.numRows();
334  Kokkos::parallel_for( nrow, ReplaceDiagonalValuesKernel(matrix) );
335  }
336 
337  // Check the result is as expected
338  static bool check(const MatrixType matrix,
339  Teuchos::FancyOStream& out) {
340  typedef typename MatrixType::values_type matrix_values_type;
341  typename matrix_values_type::HostMirror host_matrix_values =
342  Kokkos::create_mirror_view(matrix.values);
343  Kokkos::deep_copy(host_matrix_values, matrix.values);
344  const ordinal_type nrow = matrix.numRows();
345  bool success = true;
346  value_type val_expected(Kokkos::cijk(matrix.values));
347  for (ordinal_type row=0; row<nrow; ++row) {
348  val_expected = row;
349  bool s = compareVecs(host_matrix_values(row),
350  "matrix_values(row)",
351  val_expected,
352  "val_expected",
353  0.0, 0.0, out);
354  success = success && s;
355  }
356  return success;
357  }
358 };
359 
360 // Kernel to add values to the diagonal of a matrix
361 template <typename MatrixType>
362 struct AddDiagonalValuesKernel {
363  typedef typename MatrixType::execution_space execution_space;
364  typedef typename MatrixType::size_type size_type;
365  typedef typename MatrixType::value_type value_type;
366  typedef typename MatrixType::ordinal_type ordinal_type;
367 
368  const MatrixType m_matrix;
369  AddDiagonalValuesKernel(const MatrixType matrix) : m_matrix(matrix) {};
370 
371  // Replace diagonal entry for row 'i' with a value
372  KOKKOS_INLINE_FUNCTION
373  void operator() (const size_type i) const {
374  const ordinal_type row = i;
375  const ordinal_type col = i;
376  value_type val = value_type(row);
377  m_matrix.sumIntoValues(row, &col, 1, &val, false, true);
378  }
379 
380  // Kernel launch
381  static void apply(const MatrixType matrix) {
382  const size_type nrow = matrix.numRows();
383  Kokkos::parallel_for( nrow, AddDiagonalValuesKernel(matrix) );
384  }
385 
386  // Check the result is as expected
387  static bool check(const MatrixType matrix,
388  Teuchos::FancyOStream& out) {
389  typedef typename MatrixType::values_type matrix_values_type;
390  typename matrix_values_type::HostMirror host_matrix_values =
391  Kokkos::create_mirror_view(matrix.values);
392  Kokkos::deep_copy(host_matrix_values, matrix.values);
393  const ordinal_type nrow = matrix.numRows();
394  bool success = true;
395  value_type val_expected(Kokkos::cijk(matrix.values));
396  for (ordinal_type row=0; row<nrow; ++row) {
397  val_expected = row;
398  bool s = compareVecs(host_matrix_values(row),
399  "matrix_values(row)",
400  val_expected,
401  "val_expected",
402  0.0, 0.0, out);
403  success = success && s;
404  }
405  return success;
406  }
407 };
408 
409 // Kernel to add values to the diagonal of a matrix where each thread
410 // adds to the same row (checks atomic really works)
411 template <typename MatrixType>
412 struct AddDiagonalValuesAtomicKernel {
413  typedef typename MatrixType::execution_space execution_space;
414  typedef typename MatrixType::size_type size_type;
415  typedef typename MatrixType::value_type value_type;
416  typedef typename MatrixType::ordinal_type ordinal_type;
417 
418  const MatrixType m_matrix;
419  AddDiagonalValuesAtomicKernel(const MatrixType matrix) : m_matrix(matrix) {};
420 
421  // Replace diagonal entry for row 'i' with a value
422  KOKKOS_INLINE_FUNCTION
423  void operator() (const size_type i) const {
424  const ordinal_type row = 0;
425  const ordinal_type col = 0;
426  value_type val = value_type(i);
427  m_matrix.sumIntoValues(row, &col, 1, &val, false, true);
428  }
429 
430  // Kernel launch
431  static void apply(const MatrixType matrix) {
432  const size_type nrow = matrix.numRows();
433  Kokkos::parallel_for( nrow, AddDiagonalValuesAtomicKernel(matrix) );
434  }
435 
436  // Check the result is as expected
437  static bool check(const MatrixType matrix,
438  Teuchos::FancyOStream& out) {
439  typedef typename MatrixType::values_type matrix_values_type;
440  typename matrix_values_type::HostMirror host_matrix_values =
441  Kokkos::create_mirror_view(matrix.values);
442  Kokkos::deep_copy(host_matrix_values, matrix.values);
443  const ordinal_type nrow = matrix.numRows();
444  bool success = true;
445  value_type val_expected(Kokkos::cijk(matrix.values));
446  for (ordinal_type row=0; row<nrow; ++row) {
447  value_type val;
448  if (row == 0)
449  val_expected = nrow*(nrow-1)/2;
450  else
451  val_expected = 0.0;
452  bool s = compareVecs(host_matrix_values(row),
453  "matrix_values(row)",
454  val_expected,
455  "val_expected",
456  0.0, 0.0, out);
457  success = success && s;
458  }
459  return success;
460  }
461 };
462 
463 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(
464  Kokkos_CrsMatrix_PCE, ReplaceValues, MatrixScalar )
465 {
466  typedef typename MatrixScalar::ordinal_type Ordinal;
467  typedef typename MatrixScalar::execution_space Device;
468  typedef typename MatrixScalar::cijk_type Cijk;
469  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
470 
471  // Build Cijk tensor
472  const Ordinal stoch_dim = 2;
473  const Ordinal poly_ord = 3;
474  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
475 
476  // Build diagonal matrix
477  const Ordinal nrow = 10;
478  const Ordinal pce_size = cijk.dimension();
479  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
480 
481  // Launch our kernel
482  typedef ReplaceDiagonalValuesKernel<Matrix> kernel;
483  kernel::apply(matrix);
484 
485  // Check the result
486  success = kernel::check(matrix, out);
487 }
488 
489 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(
490  Kokkos_CrsMatrix_PCE, SumIntoValues, MatrixScalar )
491 {
492  typedef typename MatrixScalar::ordinal_type Ordinal;
493  typedef typename MatrixScalar::execution_space Device;
494  typedef typename MatrixScalar::cijk_type Cijk;
495  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
496 
497  // Build Cijk tensor
498  const Ordinal stoch_dim = 2;
499  const Ordinal poly_ord = 3;
500  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
501 
502  // Build diagonal matrix
503  const Ordinal nrow = 10;
504  const Ordinal pce_size = cijk.dimension();
505  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
506 
507  // Launch our kernel
508  typedef AddDiagonalValuesKernel<Matrix> kernel;
509  kernel::apply(matrix);
510 
511  // Check the result
512  success = kernel::check(matrix, out);
513 }
514 
515 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(
516  Kokkos_CrsMatrix_PCE, SumIntoValuesAtomic, MatrixScalar )
517 {
518  typedef typename MatrixScalar::ordinal_type Ordinal;
519  typedef typename MatrixScalar::execution_space Device;
520  typedef typename MatrixScalar::cijk_type Cijk;
521  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
522 
523  // Build Cijk tensor
524  const Ordinal stoch_dim = 2;
525  const Ordinal poly_ord = 3;
526  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
527 
528  // Build diagonal matrix
529  const Ordinal nrow = 10;
530  const Ordinal pce_size = cijk.dimension();
531  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
532 
533  // Launch our kernel
534  typedef AddDiagonalValuesAtomicKernel<Matrix> kernel;
535  kernel::apply(matrix);
536 
537  // Check the result
538  success = kernel::check(matrix, out);
539 }
540 
541 template <typename PCEType, typename Multiply>
542 bool test_embedded_pce(const typename PCEType::ordinal_type nGrid,
543  const typename PCEType::ordinal_type stoch_dim,
544  const typename PCEType::ordinal_type poly_ord,
545  KokkosSparse::DeviceConfig dev_config,
546  Multiply multiply_op,
547  Teuchos::FancyOStream& out)
548 {
549  typedef typename PCEType::ordinal_type ordinal_type;
550  typedef typename PCEType::value_type scalar_type;
551  typedef typename PCEType::storage_type storage_type;
552  typedef typename PCEType::cijk_type cijk_type;
553  typedef typename storage_type::execution_space execution_space;
554  typedef Kokkos::LayoutLeft Layout;
555  typedef Kokkos::View< PCEType*, Layout, execution_space > block_vector_type;
556  typedef KokkosSparse::CrsMatrix< PCEType, ordinal_type, execution_space > block_matrix_type;
557  typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
558  typedef typename block_matrix_type::values_type matrix_values_type;
559 
560  // Build Cijk tensor
561  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
562  const ordinal_type stoch_length = cijk.dimension();
563  // const ordinal_type align = 8;
564  // const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
565  const ordinal_type stoch_length_aligned = stoch_length;
566 
567  // Check pce_length == storage_type::static_size for static storage
568  TEUCHOS_TEST_FOR_EXCEPTION(
569  storage_type::is_static && storage_type::static_size != stoch_length,
570  std::logic_error,
571  "Static storage size must equal pce size");
572 
573  // Generate FEM graph:
574  const ordinal_type fem_length = nGrid * nGrid * nGrid;
575  std::vector< std::vector<ordinal_type> > fem_graph;
576  const ordinal_type fem_graph_length = generate_fem_graph( nGrid, fem_graph );
577 
578  //------------------------------
579  // Generate input/output multivectors -- Sacado dimension is always last,
580  // regardless of LayoutLeft/Right
581 
582  block_vector_type x =
583  Kokkos::make_view<block_vector_type>("x", cijk, fem_length, stoch_length_aligned);
584  block_vector_type y =
585  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, stoch_length_aligned);
586 
587  typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
588  typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
589 
590  // View the block vector as an array of the embedded intrinsic type.
591  typename block_vector_type::HostMirror::array_type hax = hx ;
592  typename block_vector_type::HostMirror::array_type hay = hy ;
593 
594  for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
595  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
596  hax(iRowStoch, iRowFEM) =
597  generate_vector_coefficient<scalar_type>(
598  fem_length, stoch_length, iRowFEM, iRowStoch );
599  hay(iRowStoch, iRowFEM) = 0.0;
600  }
601  }
602 
603  Kokkos::deep_copy( x, hx );
604  Kokkos::deep_copy( y, hy );
605 
606  //------------------------------
607  // Generate block matrix -- it is always LayoutRight (currently)
608 
609  matrix_graph_type matrix_graph =
610  Kokkos::create_staticcrsgraph<matrix_graph_type>(
611  std::string("test crs graph"), fem_graph);
612  matrix_values_type matrix_values =
613  Kokkos::make_view<matrix_values_type>(
614  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, stoch_length_aligned);
615  block_matrix_type matrix(
616  "block_matrix", fem_length, matrix_values, matrix_graph);
617  matrix.dev_config = dev_config;
618 
619  typename matrix_values_type::HostMirror hM =
620  Kokkos::create_mirror_view( matrix.values );
621 
622  typename matrix_values_type::HostMirror::array_type haM = hM ;
623 
624  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
625  const ordinal_type row_size = fem_graph[iRowFEM].size();
626  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
627  ++iRowEntryFEM, ++iEntryFEM) {
628  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
629 
630  for (ordinal_type k=0; k<stoch_length; ++k) {
631  haM(iEntryFEM, k) =
632  generate_matrix_coefficient<scalar_type>(
633  fem_length, stoch_length, iRowFEM, iColFEM, k);
634  }
635  }
636  }
637 
638  Kokkos::deep_copy( matrix.values, hM );
639 
640  //------------------------------
641  // multiply
642 
643  multiply_op( matrix, x, y );
644 
645  //------------------------------
646  // generate correct answer
647 
648  typedef typename block_vector_type::array_type array_type;
649  array_type ay_expected =
650  array_type("ay_expected", stoch_length_aligned, fem_length);
651  typename array_type::HostMirror hay_expected =
652  Kokkos::create_mirror_view(ay_expected);
653  typename cijk_type::HostMirror host_cijk =
654  Kokkos::create_mirror_view(cijk);
655  Kokkos::deep_copy(host_cijk, cijk);
656  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
657  const ordinal_type row_size = fem_graph[iRowFEM].size();
658  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
659  ++iRowEntryFEM, ++iEntryFEM) {
660  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
661  for (ordinal_type i=0; i<stoch_length; ++i) {
662  const ordinal_type num_entry = host_cijk.num_entry(i);
663  const ordinal_type entry_beg = host_cijk.entry_begin(i);
664  const ordinal_type entry_end = entry_beg + num_entry;
665  scalar_type tmp = 0;
666  for (ordinal_type entry = entry_beg; entry < entry_end; ++entry) {
667  const ordinal_type j = host_cijk.coord(entry,0);
668  const ordinal_type k = host_cijk.coord(entry,1);
669  const scalar_type a_j =
670  generate_matrix_coefficient<scalar_type>(
671  fem_length, stoch_length, iRowFEM, iColFEM, j);
672  const scalar_type a_k =
673  generate_matrix_coefficient<scalar_type>(
674  fem_length, stoch_length, iRowFEM, iColFEM, k);
675  const scalar_type x_j =
676  generate_vector_coefficient<scalar_type>(
677  fem_length, stoch_length, iColFEM, j);
678  const scalar_type x_k =
679  generate_vector_coefficient<scalar_type>(
680  fem_length, stoch_length, iColFEM, k);
681  tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
682  }
683  hay_expected(i, iRowFEM) += tmp;
684  }
685  }
686  }
687  Kokkos::deep_copy( ay_expected, hay_expected );
688 
689  //------------------------------
690  // check
691 
692  typename block_vector_type::array_type ay = y;
693  scalar_type rel_tol = ScalarTol<scalar_type>::tol();
694  scalar_type abs_tol = ScalarTol<scalar_type>::tol();
695  bool success = compare_rank_2_views(ay, ay_expected, rel_tol, abs_tol, out);
696 
697  return success;
698 }
699 
700 TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(
701  Kokkos_CrsMatrix_PCE, Multiply, Scalar, MultiplyOp )
702 {
703  typedef typename Scalar::ordinal_type Ordinal;
704 
705  const Ordinal nGrid = 5;
706  const Ordinal stoch_dim = 2;
707  const Ordinal poly_ord = 3;
708  KokkosSparse::DeviceConfig dev_config;
709 
710  success = test_embedded_pce<Scalar>(
711  nGrid, stoch_dim, poly_ord, dev_config, MultiplyOp(), out);
712 }
713 
714 struct Kokkos_MV_Multiply_Op {
715  template <typename Matrix, typename InputVector, typename OutputVector>
716  void operator() (const Matrix& A,
717  const InputVector& x,
718  OutputVector& y) const {
719  KokkosSparse::spmv("N", typename Matrix::value_type(1.0) , A, x, typename Matrix::value_type(0.0), y);
720  }
721 };
722 
723 template <typename Tag>
724 struct Stokhos_MV_Multiply_Op {
725  Tag tag;
726  Stokhos_MV_Multiply_Op(const Tag& tg = Tag()) : tag(tg) {}
727 
728  template <typename Matrix, typename InputVector, typename OutputVector>
729  void operator() (const Matrix& A,
730  const InputVector& x,
731  OutputVector& y) const {
732  Stokhos::multiply(A, x, y, tag);
733  }
734 };
735 
736 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(
737  Kokkos_CrsMatrix_PCE, MeanMultiplyRank1, Scalar )
738 {
739  typedef typename Scalar::ordinal_type Ordinal;
740 
741  const Ordinal nGrid = 5;
742  const Ordinal stoch_dim = 5;
743  const Ordinal poly_ord = 3;
744  KokkosSparse::DeviceConfig dev_config;
745 
746  typedef typename Scalar::ordinal_type ordinal_type;
747  typedef typename Scalar::value_type scalar_type;
748  typedef typename Scalar::storage_type storage_type;
749  typedef typename Scalar::cijk_type cijk_type;
750  typedef typename storage_type::execution_space execution_space;
751  typedef Kokkos::LayoutLeft Layout;
752  typedef Kokkos::View< Scalar*, Layout, execution_space > block_vector_type;
753  typedef KokkosSparse::CrsMatrix< Scalar, ordinal_type, execution_space > block_matrix_type;
754  typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
755  typedef typename block_matrix_type::values_type matrix_values_type;
756 
757  // Build Cijk tensor
758  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
759  cijk_type mean_cijk =
760  Stokhos::create_mean_based_product_tensor<execution_space, ordinal_type, scalar_type>();
761  const ordinal_type stoch_length = cijk.dimension();
762  const ordinal_type align = 8;
763  const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
764 
765  // Generate FEM graph:
766  const ordinal_type fem_length = nGrid * nGrid * nGrid;
767  std::vector< std::vector<ordinal_type> > fem_graph;
768  const ordinal_type fem_graph_length = generate_fem_graph( nGrid, fem_graph );
769 
770  block_vector_type x =
771  Kokkos::make_view<block_vector_type>("x", cijk, fem_length, stoch_length_aligned);
772  block_vector_type y =
773  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, stoch_length_aligned);
774 
775  block_vector_type y_expected =
776  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, stoch_length_aligned);
777 
778  typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
779  typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
780  typename block_vector_type::HostMirror hy_expected =
781  Kokkos::create_mirror_view( y_expected );
782 
783  // View the block vector as an array of the embedded intrinsic type.
784  typename block_vector_type::HostMirror::array_type hax = hx ;
785  typename block_vector_type::HostMirror::array_type hay = hy ;
786  typename block_vector_type::HostMirror::array_type hay_expected =
787  hy_expected ;
788 
789  for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
790  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
791  hax(iRowStoch,iRowFEM) =
792  generate_vector_coefficient<scalar_type>(
793  fem_length, stoch_length, iRowFEM, iRowStoch );
794  hay(iRowStoch,iRowFEM) = 0.0;
795  hay_expected(iRowStoch,iRowFEM) = 0.0;
796  }
797  }
798  Kokkos::deep_copy( x, hx );
799  Kokkos::deep_copy( y, hy );
800  Kokkos::deep_copy( y_expected, hy_expected );
801 
802  //------------------------------
803  // Generate block matrix -- it is always LayoutRight (currently)
804  matrix_graph_type matrix_graph =
805  Kokkos::create_staticcrsgraph<matrix_graph_type>(
806  std::string("test crs graph"), fem_graph);
807  matrix_values_type matrix_values =
808  Kokkos::make_view<matrix_values_type>(
809  Kokkos::ViewAllocateWithoutInitializing("matrix"), mean_cijk, fem_graph_length, ordinal_type(1)); //instead of stoch_length
810  block_matrix_type matrix(
811  "block_matrix", fem_length, matrix_values, matrix_graph);
812  matrix.dev_config = dev_config;
813 
814  typename matrix_values_type::HostMirror hM =
815  Kokkos::create_mirror_view( matrix.values );
816  typename matrix_values_type::HostMirror::array_type haM = hM ;
817 
818  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
819  const ordinal_type row_size = fem_graph[iRowFEM].size();
820  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
821  ++iRowEntryFEM, ++iEntryFEM) {
822  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
823 
824  haM(iEntryFEM, 0) =
825  generate_matrix_coefficient<scalar_type>(
826  fem_length, 1, iRowFEM, iColFEM, 0);
827  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
828  hay_expected(iRowStoch,iRowFEM) +=
829  haM(iEntryFEM, 0) * hax(iRowStoch,iColFEM);
830  }
831  }
832  }
833  Kokkos::deep_copy( matrix.values, hM );
834  Kokkos::deep_copy( y_expected, hy_expected );
835 
836  /*
837  //Generate same matrix with stochastic dim = Kokkos::dimension_scalar(x) (i.e. not = 1)
838  matrix_values_type full_matrix_values =
839  Kokkos::make_view<matrix_values_type>(
840  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, stoch_length_aligned);
841  block_matrix_type full_matrix(
842  "block_matrix", fem_length, full_matrix_values, matrix_graph);
843  matrix.dev_config = dev_config;
844 
845  typename matrix_values_type::HostMirror full_hM =
846  Kokkos::create_mirror_view( full_matrix.values );
847  typename matrix_values_type::HostMirror::array_type full_haM = full_hM ;
848 
849  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
850  const ordinal_type row_size = fem_graph[iRowFEM].size();
851  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
852  ++iRowEntryFEM, ++iEntryFEM) {
853  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
854 
855  for (ordinal_type k=0; k<stoch_length; ++k) {
856  if (k == 0)
857  full_haM(iEntryFEM, k) =
858  generate_matrix_coefficient<scalar_type>(
859  fem_length, 1, iRowFEM, iColFEM, k);
860  else
861  full_haM(iEntryFEM, k) = 0.0;
862  }
863  }
864  }
865 
866  Kokkos::deep_copy( full_matrix.values, full_hM );
867  */
868 
869  //------------------------------
870  // multiply
871 
872  KokkosSparse::spmv("N", Scalar(1.0) , matrix, x, Scalar(0.0), y);
873 
874  //------------------------------
875  // multiply with same matrix but with sacado_size = x.sacado_size
876 
877  //Kokkos::MV_Multiply( y_expected, full_matrix, x );
878 
879  //------------------------------
880  // check
881  scalar_type rel_tol = ScalarTol<scalar_type>::tol();
882  scalar_type abs_tol = ScalarTol<scalar_type>::tol();
883  success = compareRank1(y, y_expected, rel_tol, abs_tol, out);
884 }
885 
886 
887 TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(
888  Kokkos_CrsMatrix_PCE, MeanMultiplyRank2, Scalar )
889 {
890  typedef typename Scalar::ordinal_type ordinal_type;
891  typedef typename Scalar::value_type scalar_type;
892  typedef typename Scalar::storage_type storage_type;
893  typedef typename Scalar::cijk_type cijk_type;
894  typedef typename storage_type::execution_space execution_space;
895  typedef Kokkos::LayoutLeft Layout;
896  typedef Kokkos::View< Scalar**, Layout, execution_space > block_vector_type;
897  typedef KokkosSparse::CrsMatrix< Scalar, ordinal_type, execution_space > block_matrix_type;
898  typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
899  typedef typename block_matrix_type::values_type matrix_values_type;
900 
901 
902  const ordinal_type nGrid = 5;
903  const ordinal_type stoch_dim = 2;
904  const ordinal_type poly_ord = 3;
905  KokkosSparse::DeviceConfig dev_config;
906 
907  // Build Cijk tensor
908  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
909  cijk_type mean_cijk =
910  Stokhos::create_mean_based_product_tensor<execution_space, ordinal_type, scalar_type>();
911  const ordinal_type stoch_length = cijk.dimension();
912  const ordinal_type align = 8;
913  const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
914  const ordinal_type num_cols = 2;
915  // Generate FEM graph:
916  const ordinal_type fem_length = nGrid * nGrid * nGrid;
917  std::vector< std::vector<ordinal_type> > fem_graph;
918  const ordinal_type fem_graph_length = generate_fem_graph( nGrid, fem_graph );
919 
920  block_vector_type x =
921  Kokkos::make_view<block_vector_type>("x", cijk, fem_length, num_cols, stoch_length_aligned);
922  block_vector_type y =
923  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, num_cols, stoch_length_aligned);
924 
925  block_vector_type y_expected =
926  Kokkos::make_view<block_vector_type>("y_expected", cijk, fem_length, num_cols, stoch_length_aligned);
927 
928  typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
929  typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
930  typename block_vector_type::HostMirror hy_expected =
931  Kokkos::create_mirror_view( y_expected );
932 
933  for (ordinal_type i=0; i<num_cols; ++i){
934  for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
935  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
936  hx(iRowFEM,i).fastAccessCoeff(iRowStoch) =
937  generate_vector_coefficient<scalar_type>(
938  fem_length, stoch_length, iRowFEM, iRowStoch );
939  hy(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
940  hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
941  }
942  }
943  }
944  Kokkos::deep_copy( x, hx );
945  Kokkos::deep_copy( y, hy );
946  Kokkos::deep_copy( y_expected, hy_expected );
947 
948  //------------------------------
949  // Generate matrix with stochastic dimension 1
950  matrix_graph_type matrix_graph =
951  Kokkos::create_staticcrsgraph<matrix_graph_type>(
952  std::string("test crs graph"), fem_graph);
953  matrix_values_type matrix_values =
954  Kokkos::make_view<matrix_values_type>(
955  Kokkos::ViewAllocateWithoutInitializing("matrix"), mean_cijk, fem_graph_length, ordinal_type(1));
956  block_matrix_type matrix(
957  "block_matrix", fem_length, matrix_values, matrix_graph);
958  matrix.dev_config = dev_config;
959 
960  typename matrix_values_type::HostMirror hM =
961  Kokkos::create_mirror_view( matrix.values );
962 
963  typename matrix_values_type::HostMirror::array_type haM = hM ;
964 
965  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
966  const ordinal_type row_size = fem_graph[iRowFEM].size();
967  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
968  ++iRowEntryFEM, ++iEntryFEM) {
969  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
970 
971  haM(iEntryFEM, 0) =
972  generate_matrix_coefficient<scalar_type>(
973  fem_length, 1, iRowFEM, iColFEM, 0);
974  for (ordinal_type i=0; i<num_cols; ++i){
975  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
976  hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) +=
977  haM(iEntryFEM, 0) * hx(iColFEM,i).fastAccessCoeff(iRowStoch);
978  }
979  }
980  }
981  }
982 
983  Kokkos::deep_copy( matrix.values, hM );
984  Kokkos::deep_copy( y_expected, hy_expected );
985 
986  /*
987  //Generate same matrix with stochastic dim = Kokkos::dimension_scalar(x) (i.e. not = 1)
988  matrix_values_type full_matrix_values =
989  Kokkos::make_view<matrix_values_type>(
990  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, stoch_length_aligned);
991  block_matrix_type full_matrix(
992  "block_matrix", fem_length, full_matrix_values, matrix_graph);
993  matrix.dev_config = dev_config;
994 
995  typename matrix_values_type::HostMirror full_hM =
996  Kokkos::create_mirror_view( full_matrix.values );
997 
998  typename matrix_values_type::HostMirror::array_type full_haM = full_hM ;
999 
1000  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
1001  const ordinal_type row_size = fem_graph[iRowFEM].size();
1002  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
1003  ++iRowEntryFEM, ++iEntryFEM) {
1004  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
1005 
1006  for (ordinal_type k=0; k<stoch_length; ++k) {
1007  if (k == 0)
1008  full_haM(iEntryFEM, k) =
1009  generate_matrix_coefficient<scalar_type>(
1010  fem_length, 1, iRowFEM, iColFEM, k);
1011  else
1012  full_haM(iEntryFEM, k) = 0.0;
1013  }
1014  }
1015  }
1016 
1017  Kokkos::deep_copy( full_matrix.values, full_hM );
1018  */
1019 
1020  //------------------------------
1021  // multiply
1022 
1023  KokkosSparse::spmv("N", Scalar(1.0) , matrix, x, Scalar(0.0), y);
1024 
1025  //------------------------------
1026  // multiply with full matrix
1027 
1028  //Kokkos::MV_Multiply( y_expected, full_matrix, x );
1029 
1030  //------------------------------
1031  // check
1032 
1033  scalar_type rel_tol = ScalarTol<scalar_type>::tol();
1034  scalar_type abs_tol = ScalarTol<scalar_type>::tol();
1035  success = compareRank2(y, y_expected, rel_tol, abs_tol, out);
1036 }
1037 
1038 
1039 typedef Kokkos_MV_Multiply_Op KokkosMultiply;
1040 
1041 #define CRSMATRIX_UQ_PCE_TESTS_MATRIXSCALAR( SCALAR ) \
1042  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1043  Kokkos_CrsMatrix_PCE, ReplaceValues, SCALAR ) \
1044  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1045  Kokkos_CrsMatrix_PCE, SumIntoValues, SCALAR ) \
1046  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1047  Kokkos_CrsMatrix_PCE, SumIntoValuesAtomic, SCALAR )
1048 #define CRSMATRIX_UQ_PCE_MEAN_MULTIPLY_TESTS( SCALAR ) \
1049  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1050  Kokkos_CrsMatrix_PCE, MeanMultiplyRank1, SCALAR ) \
1051  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1052  Kokkos_CrsMatrix_PCE, MeanMultiplyRank2, SCALAR )
1053 #define CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR_OP( SCALAR, OP ) \
1054  TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( \
1055  Kokkos_CrsMatrix_PCE, Multiply, SCALAR, OP )
1056 
1057 #define CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR( SCALAR ) \
1058  CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR_OP( SCALAR, KokkosMultiply )
1059 
1060 #define CRSMATRIX_UQ_PCE_TESTS_STORAGE( STORAGE ) \
1061  typedef Sacado::UQ::PCE<STORAGE> UQ_PCE_ ## STORAGE; \
1062  CRSMATRIX_UQ_PCE_TESTS_MATRIXSCALAR( UQ_PCE_ ## STORAGE ) \
1063  CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR( UQ_PCE_ ## STORAGE ) \
1064  CRSMATRIX_UQ_PCE_MEAN_MULTIPLY_TESTS( UQ_PCE_ ## STORAGE )
1065 
1066 #define CRSMATRIX_UQ_PCE_TESTS_ORDINAL_SCALAR_DEVICE( ORDINAL, SCALAR, DEVICE ) \
1067  typedef Stokhos::DynamicStorage<ORDINAL,SCALAR,DEVICE> DS; \
1068  CRSMATRIX_UQ_PCE_TESTS_STORAGE( DS )
1069 
1070 #define CRSMATRIX_UQ_PCE_TESTS_DEVICE( DEVICE ) \
1071  CRSMATRIX_UQ_PCE_TESTS_ORDINAL_SCALAR_DEVICE( int, double, DEVICE )
AddDiagonalValuesAtomicKernel::operator()
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:308
Stokhos_Sparse3Tensor.hpp
ReplaceDiagonalValuesKernel::operator()
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:211
storage_type
Stokhos::StandardStorage< int, double > storage_type
Definition: Stokhos_SacadoETPCEUnitTest.cpp:50
Teuchos_UnitTestHelpers.hpp
Stokhos::compareVecs
bool compareVecs(const VectorType1 &a1, const std::string &a1_name, const VectorType2 &a2, const std::string &a2_name, const ValueType &rel_tol, const ValueType &abs_tol, Teuchos::FancyOStream &out)
Definition: Stokhos_UnitTestHelpers.hpp:123
stoch_dim
const int stoch_dim
Definition: Stokhos_TpetraCrsMatrixUQPCEUnitTest.hpp:179
poly_ord
const int poly_ord
Definition: Stokhos_TpetraCrsMatrixUQPCEUnitTest.hpp:180
Stokhos::Sparse3Tensor
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Definition: Stokhos_Sparse3Tensor.hpp:56
test_embedded_pce
bool test_embedded_pce(const typename PCEType::ordinal_type nGrid, const typename PCEType::ordinal_type stoch_dim, const typename PCEType::ordinal_type poly_ord, KokkosSparse::DeviceConfig dev_config, Multiply multiply_op, Teuchos::FancyOStream &out)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:542
generate_fem_graph
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
Definition: TestEpetra.cpp:77
AddDiagonalValuesKernel::check
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:387
DynamicVecTest::execution_space
Kokkos::DefaultExecutionSpace execution_space
Definition: Stokhos_SacadoMPVectorCommTests.cpp:682
Kokkos_MV_Multiply_Op
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:584
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
buildDiagonalMatrix
MatrixType buildDiagonalMatrix(typename MatrixType::ordinal_type nrow, typename MatrixType::ordinal_type mp_vector_size)
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:174
AddDiagonalValuesAtomicKernel::apply
static void apply(const MatrixType matrix)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:431
KokkosKernelsUnitTest::generate_vector_coefficient
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
Definition: Stokhos_KokkosArrayKernelsUnitTest.hpp:152
Stokhos_MV_Multiply_Op
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:594
Kokkos::dimension_scalar
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< View< T, P...> >::value, unsigned >::type dimension_scalar(const View< T, P...> &view)
Definition: KokkosExp_View_UQ_PCE_Contiguous.hpp:208
compareRank2
bool compareRank2(const vector_type &y, const vector_type &y_exp, const scalar_type rel_tol, const scalar_type abs_tol, Teuchos::FancyOStream &out)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:242
map_fem_graph_coord
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
Definition: TestEpetra.cpp:67
j
j
Definition: Sacado_Fad_Exp_MP_Vector.hpp:527
Scalar
pce_type Scalar
Definition: linear2d_diffusion_scalar_types.hpp:49
Stokhos::multiply
void multiply(const CrsMatrix< MatrixValue, Device, Layout > &A, const InputMultiVectorType &x, OutputMultiVectorType &y, const std::vector< OrdinalType > &col_indices, SingleColumnMultivectorMultiply)
Definition: Stokhos_CrsMatrix.hpp:621
TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL
TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(Kokkos_SG_SpMv, CrsProductTensorCijk, Scalar, Device)
Definition: Stokhos_KokkosArrayKernelsUnitTest_Host.hpp:43
ReplaceDiagonalValuesKernel::check
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:338
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
KokkosMultiply
Kokkos_MV_Multiply_Op KokkosMultiply
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:606
Kokkos_MV_Multiply_Op::operator()
void operator()(const Matrix &A, const InputVector &x, OutputVector &y) const
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:586
AddDiagonalValuesKernel::apply
static void apply(const MatrixType matrix)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:381
Teuchos::basic_FancyOStream
AddDiagonalValuesAtomicKernel::check
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:437
Stokhos_MV_Multiply_Op::Stokhos_MV_Multiply_Op
Stokhos_MV_Multiply_Op(const Tag &tg=Tag())
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:726
Kokkos::deep_copy
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
Definition: Stokhos_CrsMatrix.hpp:688
Sacado::UQ::abs
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
Definition: Sacado_UQ_PCE_Imp.hpp:1181
fastAccessCoeff
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j)-expr2.val(j)
Stokhos::CompletePolynomialBasis
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Definition: Stokhos_CompletePolynomialBasis.hpp:72
build_cijk
kokkos_cijk_type build_cijk(ordinal_type stoch_dim, ordinal_type poly_ord)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:133
Stokhos_UnitTestHelpers.hpp
Stokhos_LegendreBasis.hpp
Kokkos::cijk
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Definition: KokkosExp_View_UQ_PCE_Contiguous.hpp:217
Stokhos::LegendreBasis
Legendre polynomial basis.
Definition: Stokhos_LegendreBasis.hpp:67
Kokkos_CrsMatrix_UQ_PCE.hpp
KokkosKernelsUnitTest::generate_matrix_coefficient
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
Definition: Stokhos_KokkosArrayKernelsUnitTest.hpp:135
val
expr val()
AddDiagonalValuesKernel::operator()
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:259
Stokhos::OneDOrthogPolyBasis
Abstract base class for 1-D orthogonal polynomials.
Definition: Stokhos_OneDOrthogPolyBasis.hpp:81
TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL
TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(Kokkos_CrsMatrix_MP, ReplaceValues, MatrixScalar)
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:350
compareRank1
bool compareRank1(const vector_type &y, const vector_type &y_exp, const scalar_type rel_tol, const scalar_type abs_tol, Teuchos::FancyOStream &out)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:207
Stokhos_MV_Multiply_Op::operator()
void operator()(const Matrix &A, const InputVector &x, OutputVector &y) const
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:599
Stokhos_CompletePolynomialBasis.hpp
compare_rank_2_views
bool compare_rank_2_views(const array_type &y, const array_type &y_exp, const scalar_type rel_tol, const scalar_type abs_tol, Teuchos::FancyOStream &out)
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:135
Kokkos_CrsMatrix_UQ_PCE_Cuda.hpp
Stokhos_Sacado_Kokkos_UQ_PCE.hpp
ReplaceDiagonalValuesKernel::apply
static void apply(const MatrixType matrix)
Definition: Stokhos_KokkosCrsMatrixUQPCEUnitTest.hpp:332
double
double
Definition: tpetra_mat_vec.cpp:247
Stokhos_MV_Multiply_Op::tag
Tag tag
Definition: Stokhos_KokkosCrsMatrixMPVectorUnitTest.hpp:595
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
Kokkos::create_mirror_view
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
Definition: Stokhos_CrsMatrix.hpp:677
Generated by   doxygen 1.8.5