Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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"
44 
46 #include "Kokkos_Sparse.hpp"
52 
53 // For computing DeviceConfig
54 #include "Kokkos_Core.hpp"
55 
56 // Helper functions
57 template< typename IntType >
58 inline
59 IntType map_fem_graph_coord( const IntType& N,
60  const IntType& i,
61  const IntType& j,
62  const IntType& k )
63 {
64  return k + N * ( j + N * i );
65 }
66 
67 template < typename ordinal >
68 inline
69 ordinal generate_fem_graph( ordinal N,
70  std::vector< std::vector<ordinal> >& graph )
71 {
72  graph.resize( N * N * N, std::vector<ordinal>() );
73 
74  ordinal total = 0;
75 
76  for ( int i = 0; i < (int) N; ++i ) {
77  for ( int j = 0; j < (int) N; ++j ) {
78  for ( int k = 0; k < (int) N; ++k ) {
79 
80  const ordinal row = map_fem_graph_coord((int)N,i,j,k);
81 
82  graph[row].reserve(27);
83 
84  for ( int ii = -1; ii < 2; ++ii ) {
85  for ( int jj = -1; jj < 2; ++jj ) {
86  for ( int kk = -1; kk < 2; ++kk ) {
87  if ( 0 <= i + ii && i + ii < (int) N &&
88  0 <= j + jj && j + jj < (int) N &&
89  0 <= k + kk && k + kk < (int) N ) {
90  ordinal col = map_fem_graph_coord((int)N,i+ii,j+jj,k+kk);
91 
92  graph[row].push_back(col);
93  }
94  }}}
95  total += graph[row].size();
96  }}}
97 
98  return total;
99 }
100 
101 template <typename scalar, typename ordinal>
102 inline
103 scalar generate_matrix_coefficient( const ordinal nFEM,
104  const ordinal nStoch,
105  const ordinal iRowFEM,
106  const ordinal iColFEM,
107  const ordinal iStoch )
108 {
109  const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
110  ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
111 
112  const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
113 
114  return A_fem + A_stoch;
115  //return 1.0;
116 }
117 
118 template <typename scalar, typename ordinal>
119 inline
120 scalar generate_vector_coefficient( const ordinal nFEM,
121  const ordinal nStoch,
122  const ordinal iColFEM,
123  const ordinal iStoch )
124 {
125  const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
126  const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
127  return X_fem + X_stoch;
128  //return 1.0;
129 }
130 
131 template <typename kokkos_cijk_type, typename ordinal_type>
134 {
135  using Teuchos::RCP;
136  using Teuchos::rcp;
137  using Teuchos::Array;
138 
139  typedef typename kokkos_cijk_type::value_type value_type;
145 
146  // Create product basis
147  Array< RCP<const one_d_basis> > bases(stoch_dim);
148  for (ordinal_type i=0; i<stoch_dim; i++)
149  bases[i] = rcp(new legendre_basis(poly_ord, true));
150  RCP<const product_basis> basis = rcp(new product_basis(bases));
151 
152  // Triple product tensor
153  RCP<Cijk> cijk = basis->computeTripleProductTensor();
154 
155  // Kokkos triple product tensor
156  kokkos_cijk_type kokkos_cijk =
157  Stokhos::create_product_tensor<execution_space>(*basis, *cijk);
158 
159  return kokkos_cijk;
160 }
161 
162 // Reasonable tolerances for common precisions
163 template <typename Scalar> struct ScalarTol {};
164 template <> struct ScalarTol<float> { static float tol() { return 1e-4; } };
165 template <> struct ScalarTol<double> { static double tol() { return 1e-10; } };
166 
167 // Compare two rank-2 views for equality, to given precision
168 template <typename array_type, typename scalar_type>
169 bool compare_rank_2_views(const array_type& y,
170  const array_type& y_exp,
171  const scalar_type rel_tol,
172  const scalar_type abs_tol,
173  Teuchos::FancyOStream& out)
174 {
175  typedef typename array_type::size_type size_type;
176  typename array_type::HostMirror hy = Kokkos::create_mirror_view(y);
177  typename array_type::HostMirror hy_exp = Kokkos::create_mirror_view(y_exp);
178  Kokkos::deep_copy(hy, y);
179  Kokkos::deep_copy(hy_exp, y_exp);
180 
181  size_type num_rows = y.dimension_0();
182  size_type num_cols = y.dimension_1();
183  bool success = true;
184  for (size_type i=0; i<num_rows; ++i) {
185  for (size_type j=0; j<num_cols; ++j) {
186  scalar_type diff = std::abs( hy(i,j) - hy_exp(i,j) );
187  scalar_type tol = rel_tol*std::abs(hy_exp(i,j)) + abs_tol;
188  bool s = diff < tol;
189  out << "y_expected(" << i << "," << j << ") - "
190  << "y(" << i << "," << j << ") = " << hy_exp(i,j)
191  << " - " << hy(i,j) << " == "
192  << diff << " < " << tol << " : ";
193  if (s)
194  out << "passed";
195  else
196  out << "failed";
197  out << std::endl;
198  success = success && s;
199  }
200  }
201 
202  return success;
203 }
204 
205 template <typename vector_type, typename scalar_type>
206 bool compareRank1(const vector_type& y,
207  const vector_type& y_exp,
208  const scalar_type rel_tol,
209  const scalar_type abs_tol,
210  Teuchos::FancyOStream& out)
211 {
212  typedef typename vector_type::size_type size_type;
213  typename vector_type::HostMirror hy = Kokkos::create_mirror_view(y);
214  typename vector_type::HostMirror hy_exp = Kokkos::create_mirror_view(y_exp);
215  Kokkos::deep_copy(hy, y);
216  Kokkos::deep_copy(hy_exp, y_exp);
217 
218  size_type num_rows = y.dimension_0();
219  bool success = true;
220  for (size_type i=0; i<num_rows; ++i) {
221  for (size_type j=0; j<Kokkos::dimension_scalar(y); ++j) {
222  scalar_type diff = std::abs( hy(i).fastAccessCoeff(j) - hy_exp(i).fastAccessCoeff(j) );
223  scalar_type tol = rel_tol*std::abs(hy_exp(i).fastAccessCoeff(j)) + abs_tol;
224  bool s = diff < tol;
225  out << "y_expected(" << i << ").coeff(" << j << ") - "
226  << "y(" << i << ").coeff(" << j << ") = " << hy_exp(i).fastAccessCoeff(j)
227  << " - " << hy(i).fastAccessCoeff(j) << " == "
228  << diff << " < " << tol << " : ";
229  if (s)
230  out << "passed";
231  else
232  out << "failed";
233  out << std::endl;
234  success = success && s;
235  }
236  }
237  return success;
238 }
239 
240 template <typename vector_type, typename scalar_type>
241 bool compareRank2(const vector_type& y,
242  const vector_type& y_exp,
243  const scalar_type rel_tol,
244  const scalar_type abs_tol,
245  Teuchos::FancyOStream& out)
246 {
247  typedef typename vector_type::size_type size_type;
248  typename vector_type::HostMirror hy = Kokkos::create_mirror_view(y);
249  typename vector_type::HostMirror hy_exp = Kokkos::create_mirror_view(y_exp);
250  Kokkos::deep_copy(hy, y);
251  Kokkos::deep_copy(hy_exp, y_exp);
252 
253  size_type num_rows = y.dimension_0();
254  size_type num_cols = y.dimension_1();
255  bool success = true;
256 
257  for (size_type col = 0; col < num_cols; ++col){
258  for (size_type i=0; i<num_rows; ++i) {
259  for (size_type j=0; j<Kokkos::dimension_scalar(y); ++j) {
260  scalar_type diff = std::abs( hy(i,col).fastAccessCoeff(j) - hy_exp(i,col).fastAccessCoeff(j) );
261  scalar_type tol = rel_tol*std::abs(hy_exp(i,col).fastAccessCoeff(j)) + abs_tol;
262  bool s = diff < tol;
263  out << "y_expected(" << i << ").coeff(" << j << ") - "
264  << "y(" << i << ").coeff(" << j << ") = " << hy_exp(i,col).fastAccessCoeff(j)
265  << " - " << hy(i,col).fastAccessCoeff(j) << " == "
266  << diff << " < " << tol << " : ";
267  if (s)
268  out << "passed";
269  else
270  out << "failed";
271  out << std::endl;
272  success = success && s;
273  }
274  }
275  }
276 
277 
278  return success;
279 }
280 
281 
282 // Helper function to build a diagonal matrix
283 template <typename MatrixType, typename CijkType>
284 MatrixType
286  typename MatrixType::ordinal_type pce_size,
287  const CijkType& cijk) {
288  typedef typename MatrixType::ordinal_type ordinal_type;
289  typedef typename MatrixType::StaticCrsGraphType matrix_graph_type;
290  typedef typename MatrixType::values_type matrix_values_type;
291 
292  std::vector< std::vector<ordinal_type> > graph(nrow);
293  for (ordinal_type i=0; i<nrow; ++i)
294  graph[i] = std::vector<ordinal_type>(1, i);
295  ordinal_type graph_length = nrow;
296 
297  matrix_graph_type matrix_graph =
298  Kokkos::create_staticcrsgraph<matrix_graph_type>("graph", graph);
299  matrix_values_type matrix_values =
300  Kokkos::make_view<matrix_values_type>("values", cijk, graph_length, pce_size);
301 
302  MatrixType matrix("matrix", nrow, matrix_values, matrix_graph);
303  return matrix;
304 }
305 
306 //
307 // Tests
308 //
309 
310 // Kernel to set diagonal of a matrix to prescribed values
311 template <typename MatrixType>
314  typedef typename MatrixType::size_type size_type;
317 
318  const MatrixType m_matrix;
319  ReplaceDiagonalValuesKernel(const MatrixType matrix) : m_matrix(matrix) {};
320 
321  // Replace diagonal entry for row 'i' with a value
322  KOKKOS_INLINE_FUNCTION
323  void operator() (const size_type i) const {
324  const ordinal_type row = i;
325  const ordinal_type col = i;
326  value_type val = value_type(row);
327  m_matrix.replaceValues(row, &col, 1, &val, false, true);
328  }
329 
330  // Kernel launch
331  static void apply(const MatrixType matrix) {
332  const size_type nrow = matrix.numRows();
333  Kokkos::parallel_for( nrow, ReplaceDiagonalValuesKernel(matrix) );
334  }
335 
336  // Check the result is as expected
337  static bool check(const MatrixType matrix,
338  Teuchos::FancyOStream& out) {
339  typedef typename MatrixType::values_type matrix_values_type;
340  typename matrix_values_type::HostMirror host_matrix_values =
341  Kokkos::create_mirror_view(matrix.values);
342  Kokkos::deep_copy(host_matrix_values, matrix.values);
343  const ordinal_type nrow = matrix.numRows();
344  const ordinal_type pce_size = Kokkos::dimension_scalar(host_matrix_values);
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>
364  typedef typename MatrixType::size_type size_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  const ordinal_type pce_size = Kokkos::dimension_scalar(host_matrix_values);
395  bool success = true;
396  value_type val_expected(Kokkos::cijk(matrix.values));
397  for (ordinal_type row=0; row<nrow; ++row) {
398  val_expected = row;
399  bool s = compareVecs(host_matrix_values(row),
400  "matrix_values(row)",
401  val_expected,
402  "val_expected",
403  0.0, 0.0, out);
404  success = success && s;
405  }
406  return success;
407  }
408 };
409 
410 // Kernel to add values to the diagonal of a matrix where each thread
411 // adds to the same row (checks atomic really works)
412 template <typename MatrixType>
415  typedef typename MatrixType::size_type size_type;
418 
419  const MatrixType m_matrix;
420  AddDiagonalValuesAtomicKernel(const MatrixType matrix) : m_matrix(matrix) {};
421 
422  // Replace diagonal entry for row 'i' with a value
423  KOKKOS_INLINE_FUNCTION
424  void operator() (const size_type i) const {
425  const ordinal_type row = 0;
426  const ordinal_type col = 0;
428  m_matrix.sumIntoValues(row, &col, 1, &val, false, true);
429  }
430 
431  // Kernel launch
432  static void apply(const MatrixType matrix) {
433  const size_type nrow = matrix.numRows();
434  Kokkos::parallel_for( nrow, AddDiagonalValuesAtomicKernel(matrix) );
435  }
436 
437  // Check the result is as expected
438  static bool check(const MatrixType matrix,
439  Teuchos::FancyOStream& out) {
440  typedef typename MatrixType::values_type matrix_values_type;
441  typename matrix_values_type::HostMirror host_matrix_values =
442  Kokkos::create_mirror_view(matrix.values);
443  Kokkos::deep_copy(host_matrix_values, matrix.values);
444  const ordinal_type nrow = matrix.numRows();
445  const ordinal_type pce_size = Kokkos::dimension_scalar(host_matrix_values);
446  bool success = true;
447  value_type val_expected(Kokkos::cijk(matrix.values));
448  for (ordinal_type row=0; row<nrow; ++row) {
449  value_type val;
450  if (row == 0)
451  val_expected = nrow*(nrow-1)/2;
452  else
453  val_expected = 0.0;
454  bool s = compareVecs(host_matrix_values(row),
455  "matrix_values(row)",
456  val_expected,
457  "val_expected",
458  0.0, 0.0, out);
459  success = success && s;
460  }
461  return success;
462  }
463 };
464 
466  Kokkos_CrsMatrix_PCE, ReplaceValues, MatrixScalar )
467 {
468  typedef typename MatrixScalar::ordinal_type Ordinal;
469  typedef typename MatrixScalar::execution_space Device;
470  typedef typename MatrixScalar::cijk_type Cijk;
471  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
472 
473  // Build Cijk tensor
474  const Ordinal stoch_dim = 2;
475  const Ordinal poly_ord = 3;
476  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
477 
478  // Build diagonal matrix
479  const Ordinal nrow = 10;
480  const Ordinal pce_size = cijk.dimension();
481  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
482 
483  // Launch our kernel
485  kernel::apply(matrix);
486 
487  // Check the result
488  success = kernel::check(matrix, out);
489 }
490 
492  Kokkos_CrsMatrix_PCE, SumIntoValues, MatrixScalar )
493 {
494  typedef typename MatrixScalar::ordinal_type Ordinal;
495  typedef typename MatrixScalar::execution_space Device;
496  typedef typename MatrixScalar::cijk_type Cijk;
497  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
498 
499  // Build Cijk tensor
500  const Ordinal stoch_dim = 2;
501  const Ordinal poly_ord = 3;
502  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
503 
504  // Build diagonal matrix
505  const Ordinal nrow = 10;
506  const Ordinal pce_size = cijk.dimension();
507  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
508 
509  // Launch our kernel
510  typedef AddDiagonalValuesKernel<Matrix> kernel;
511  kernel::apply(matrix);
512 
513  // Check the result
514  success = kernel::check(matrix, out);
515 }
516 
518  Kokkos_CrsMatrix_PCE, SumIntoValuesAtomic, MatrixScalar )
519 {
520  typedef typename MatrixScalar::ordinal_type Ordinal;
521  typedef typename MatrixScalar::execution_space Device;
522  typedef typename MatrixScalar::cijk_type Cijk;
523  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
524 
525  // Build Cijk tensor
526  const Ordinal stoch_dim = 2;
527  const Ordinal poly_ord = 3;
528  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
529 
530  // Build diagonal matrix
531  const Ordinal nrow = 10;
532  const Ordinal pce_size = cijk.dimension();
533  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
534 
535  // Launch our kernel
537  kernel::apply(matrix);
538 
539  // Check the result
540  success = kernel::check(matrix, out);
541 }
542 
543 template <typename PCEType, typename Multiply>
544 bool test_embedded_pce(const typename PCEType::ordinal_type nGrid,
545  const typename PCEType::ordinal_type stoch_dim,
546  const typename PCEType::ordinal_type poly_ord,
547  Kokkos::DeviceConfig dev_config,
548  Multiply multiply_op,
549  Teuchos::FancyOStream& out)
550 {
551  typedef typename PCEType::ordinal_type ordinal_type;
552  typedef typename PCEType::value_type scalar_type;
553  typedef typename PCEType::storage_type storage_type;
554  typedef typename PCEType::cijk_type cijk_type;
556  typedef Kokkos::LayoutLeft Layout;
557  typedef Kokkos::View< PCEType*, Layout, execution_space > block_vector_type;
558  typedef KokkosSparse::CrsMatrix< PCEType, ordinal_type, execution_space > block_matrix_type;
559  typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
560  typedef typename block_matrix_type::values_type matrix_values_type;
561 
562  // Build Cijk tensor
563  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
564  const ordinal_type stoch_length = cijk.dimension();
565  // const ordinal_type align = 8;
566  // const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
567  const ordinal_type stoch_length_aligned = stoch_length;
568 
569  // Check pce_length == storage_type::static_size for static storage
570  TEUCHOS_TEST_FOR_EXCEPTION(
571  storage_type::is_static && storage_type::static_size != stoch_length,
572  std::logic_error,
573  "Static storage size must equal pce size");
574 
575  // Generate FEM graph:
576  const ordinal_type fem_length = nGrid * nGrid * nGrid;
577  std::vector< std::vector<ordinal_type> > fem_graph;
578  const ordinal_type fem_graph_length = generate_fem_graph( nGrid, fem_graph );
579 
580  //------------------------------
581  // Generate input/output multivectors -- Sacado dimension is always last,
582  // regardless of LayoutLeft/Right
583 
584  block_vector_type x =
585  Kokkos::make_view<block_vector_type>("x", cijk, fem_length, stoch_length_aligned);
586  block_vector_type y =
587  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, stoch_length_aligned);
588 
589  typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
590  typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
591 
592  // View the block vector as an array of the embedded intrinsic type.
593  typename block_vector_type::HostMirror::array_type hax = hx ;
594  typename block_vector_type::HostMirror::array_type hay = hy ;
595 
596  for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
597  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
598  hax(iRowStoch, iRowFEM) =
599  generate_vector_coefficient<scalar_type>(
600  fem_length, stoch_length, iRowFEM, iRowStoch );
601  hay(iRowStoch, iRowFEM) = 0.0;
602  }
603  }
604 
605  Kokkos::deep_copy( x, hx );
606  Kokkos::deep_copy( y, hy );
607 
608  //------------------------------
609  // Generate block matrix -- it is always LayoutRight (currently)
610 
611  matrix_graph_type matrix_graph =
612  Kokkos::create_staticcrsgraph<matrix_graph_type>(
613  std::string("test crs graph"), fem_graph);
614  matrix_values_type matrix_values =
615  Kokkos::make_view<matrix_values_type>(
616  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, stoch_length_aligned);
617  block_matrix_type matrix(
618  "block_matrix", fem_length, matrix_values, matrix_graph);
619  matrix.dev_config = dev_config;
620 
621  typename matrix_values_type::HostMirror hM =
622  Kokkos::create_mirror_view( matrix.values );
623 
624  typename matrix_values_type::HostMirror::array_type haM = hM ;
625 
626  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
627  const ordinal_type row_size = fem_graph[iRowFEM].size();
628  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
629  ++iRowEntryFEM, ++iEntryFEM) {
630  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
631 
632  for (ordinal_type k=0; k<stoch_length; ++k) {
633  haM(iEntryFEM, k) =
634  generate_matrix_coefficient<scalar_type>(
635  fem_length, stoch_length, iRowFEM, iColFEM, k);
636  }
637  }
638  }
639 
640  Kokkos::deep_copy( matrix.values, hM );
641 
642  //------------------------------
643  // multiply
644 
645  multiply_op( matrix, x, y );
646 
647  //------------------------------
648  // generate correct answer
649 
650  typedef typename block_vector_type::array_type array_type;
651  array_type ay_expected =
652  array_type("ay_expected", stoch_length_aligned, fem_length);
653  typename array_type::HostMirror hay_expected =
654  Kokkos::create_mirror_view(ay_expected);
655  typename cijk_type::HostMirror host_cijk =
657  Kokkos::deep_copy(host_cijk, cijk);
658  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
659  const ordinal_type row_size = fem_graph[iRowFEM].size();
660  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
661  ++iRowEntryFEM, ++iEntryFEM) {
662  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
663  for (ordinal_type i=0; i<stoch_length; ++i) {
664  const ordinal_type num_entry = host_cijk.num_entry(i);
665  const ordinal_type entry_beg = host_cijk.entry_begin(i);
666  const ordinal_type entry_end = entry_beg + num_entry;
667  scalar_type tmp = 0;
668  for (ordinal_type entry = entry_beg; entry < entry_end; ++entry) {
669  const ordinal_type j = host_cijk.coord(entry,0);
670  const ordinal_type k = host_cijk.coord(entry,1);
671  const scalar_type a_j =
672  generate_matrix_coefficient<scalar_type>(
673  fem_length, stoch_length, iRowFEM, iColFEM, j);
674  const scalar_type a_k =
675  generate_matrix_coefficient<scalar_type>(
676  fem_length, stoch_length, iRowFEM, iColFEM, k);
677  const scalar_type x_j =
678  generate_vector_coefficient<scalar_type>(
679  fem_length, stoch_length, iColFEM, j);
680  const scalar_type x_k =
681  generate_vector_coefficient<scalar_type>(
682  fem_length, stoch_length, iColFEM, k);
683  tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
684  }
685  hay_expected(i, iRowFEM) += tmp;
686  }
687  }
688  }
689  Kokkos::deep_copy( ay_expected, hay_expected );
690 
691  //------------------------------
692  // check
693 
694  typename block_vector_type::array_type ay = y;
697  bool success = compare_rank_2_views(ay, ay_expected, rel_tol, abs_tol, out);
698 
699  return success;
700 }
701 
703  Kokkos_CrsMatrix_PCE, Multiply, Scalar, MultiplyOp )
704 {
705  typedef typename Scalar::ordinal_type Ordinal;
706 
707  const Ordinal nGrid = 5;
708  const Ordinal stoch_dim = 2;
709  const Ordinal poly_ord = 3;
710  Kokkos::DeviceConfig dev_config;
711 
712  success = test_embedded_pce<Scalar>(
713  nGrid, stoch_dim, poly_ord, dev_config, MultiplyOp(), out);
714 }
715 
716 struct Kokkos_MV_Multiply_Op {
717  template <typename Matrix, typename InputVector, typename OutputVector>
718  void operator() (const Matrix& A,
719  const InputVector& x,
720  OutputVector& y) const {
721  KokkosSparse::spmv("N", typename Matrix::value_type(1.0) , A, x, typename Matrix::value_type(0.0), y);
722  }
723 };
724 
725 template <typename Tag>
726 struct Stokhos_MV_Multiply_Op {
727  Tag tag;
728  Stokhos_MV_Multiply_Op(const Tag& tg = Tag()) : tag(tg) {}
729 
730  template <typename Matrix, typename InputVector, typename OutputVector>
731  void operator() (const Matrix& A,
732  const InputVector& x,
733  OutputVector& y) const {
734  Stokhos::multiply(A, x, y, tag);
735  }
736 };
737 
739  Kokkos_CrsMatrix_PCE, MeanMultiplyRank1, Scalar )
740 {
741  typedef typename Scalar::ordinal_type Ordinal;
742 
743  const Ordinal nGrid = 5;
744  const Ordinal stoch_dim = 5;
745  const Ordinal poly_ord = 3;
746  Kokkos::DeviceConfig dev_config;
747 
748  typedef typename Scalar::ordinal_type ordinal_type;
749  typedef typename Scalar::value_type scalar_type;
750  typedef typename Scalar::storage_type storage_type;
751  typedef typename Scalar::cijk_type cijk_type;
753  typedef Kokkos::LayoutLeft Layout;
754  typedef Kokkos::View< Scalar*, Layout, execution_space > block_vector_type;
755  typedef KokkosSparse::CrsMatrix< Scalar, ordinal_type, execution_space > block_matrix_type;
756  typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
757  typedef typename block_matrix_type::values_type matrix_values_type;
758 
759  // Build Cijk tensor
760  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
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"), 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
883  success = compareRank1(y, y_expected, rel_tol, abs_tol, out);
884 }
885 
886 
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;
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  Kokkos::DeviceConfig dev_config;
906 
907  // Build Cijk tensor
908  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
909  const ordinal_type stoch_length = cijk.dimension();
910  const ordinal_type align = 8;
911  const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
912  const ordinal_type num_cols = 2;
913  // Generate FEM graph:
914  const ordinal_type fem_length = nGrid * nGrid * nGrid;
915  std::vector< std::vector<ordinal_type> > fem_graph;
916  const ordinal_type fem_graph_length = generate_fem_graph( nGrid, fem_graph );
917 
918  block_vector_type x =
919  Kokkos::make_view<block_vector_type>("x", cijk, fem_length, num_cols, stoch_length_aligned);
920  block_vector_type y =
921  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, num_cols, stoch_length_aligned);
922 
923  block_vector_type y_expected =
924  Kokkos::make_view<block_vector_type>("y_expected", cijk, fem_length, num_cols, stoch_length_aligned);
925 
926  typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
927  typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
928  typename block_vector_type::HostMirror hy_expected =
929  Kokkos::create_mirror_view( y_expected );
930 
931  for (ordinal_type i=0; i<num_cols; ++i){
932  for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
933  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
934  hx(iRowFEM,i).fastAccessCoeff(iRowStoch) =
935  generate_vector_coefficient<scalar_type>(
936  fem_length, stoch_length, iRowFEM, iRowStoch );
937  hy(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
938  hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
939  }
940  }
941  }
942  Kokkos::deep_copy( x, hx );
943  Kokkos::deep_copy( y, hy );
944  Kokkos::deep_copy( y_expected, hy_expected );
945 
946  //------------------------------
947  // Generate matrix with stochastic dimension 1
948  matrix_graph_type matrix_graph =
949  Kokkos::create_staticcrsgraph<matrix_graph_type>(
950  std::string("test crs graph"), fem_graph);
951  matrix_values_type matrix_values =
952  Kokkos::make_view<matrix_values_type>(
953  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, ordinal_type(1));
954  block_matrix_type matrix(
955  "block_matrix", fem_length, matrix_values, matrix_graph);
956  matrix.dev_config = dev_config;
957 
958  typename matrix_values_type::HostMirror hM =
959  Kokkos::create_mirror_view( matrix.values );
960 
961  typename matrix_values_type::HostMirror::array_type haM = hM ;
962 
963  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
964  const ordinal_type row_size = fem_graph[iRowFEM].size();
965  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
966  ++iRowEntryFEM, ++iEntryFEM) {
967  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
968 
969  haM(iEntryFEM, 0) =
970  generate_matrix_coefficient<scalar_type>(
971  fem_length, 1, iRowFEM, iColFEM, 0);
972  for (ordinal_type i=0; i<num_cols; ++i){
973  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
974  hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) +=
975  haM(iEntryFEM, 0) * hx(iColFEM,i).fastAccessCoeff(iRowStoch);
976  }
977  }
978  }
979  }
980 
981  Kokkos::deep_copy( matrix.values, hM );
982  Kokkos::deep_copy( y_expected, hy_expected );
983 
984  /*
985  //Generate same matrix with stochastic dim = Kokkos::dimension_scalar(x) (i.e. not = 1)
986  matrix_values_type full_matrix_values =
987  Kokkos::make_view<matrix_values_type>(
988  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, stoch_length_aligned);
989  block_matrix_type full_matrix(
990  "block_matrix", fem_length, full_matrix_values, matrix_graph);
991  matrix.dev_config = dev_config;
992 
993  typename matrix_values_type::HostMirror full_hM =
994  Kokkos::create_mirror_view( full_matrix.values );
995 
996  typename matrix_values_type::HostMirror::array_type full_haM = full_hM ;
997 
998  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
999  const ordinal_type row_size = fem_graph[iRowFEM].size();
1000  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
1001  ++iRowEntryFEM, ++iEntryFEM) {
1002  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
1003 
1004  for (ordinal_type k=0; k<stoch_length; ++k) {
1005  if (k == 0)
1006  full_haM(iEntryFEM, k) =
1007  generate_matrix_coefficient<scalar_type>(
1008  fem_length, 1, iRowFEM, iColFEM, k);
1009  else
1010  full_haM(iEntryFEM, k) = 0.0;
1011  }
1012  }
1013  }
1014 
1015  Kokkos::deep_copy( full_matrix.values, full_hM );
1016  */
1017 
1018  //------------------------------
1019  // multiply
1020 
1021  KokkosSparse::spmv("N", Scalar(1.0) , matrix, x, Scalar(0.0), y);
1022 
1023  //------------------------------
1024  // multiply with full matrix
1025 
1026  //Kokkos::MV_Multiply( y_expected, full_matrix, x );
1027 
1028  //------------------------------
1029  // check
1030 
1033  success = compareRank2(y, y_expected, rel_tol, abs_tol, out);
1034 }
1035 
1036 
1038 
1039 #define CRSMATRIX_UQ_PCE_TESTS_MATRIXSCALAR( SCALAR ) \
1040  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1041  Kokkos_CrsMatrix_PCE, ReplaceValues, SCALAR ) \
1042  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1043  Kokkos_CrsMatrix_PCE, SumIntoValues, SCALAR ) \
1044  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1045  Kokkos_CrsMatrix_PCE, SumIntoValuesAtomic, SCALAR )
1046 #define CRSMATRIX_UQ_PCE_MEAN_MULTIPLY_TESTS( SCALAR ) \
1047  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1048  Kokkos_CrsMatrix_PCE, MeanMultiplyRank1, SCALAR ) \
1049  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1050  Kokkos_CrsMatrix_PCE, MeanMultiplyRank2, SCALAR )
1051 #define CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR_OP( SCALAR, OP ) \
1052  TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( \
1053  Kokkos_CrsMatrix_PCE, Multiply, SCALAR, OP )
1054 
1055 #define CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR( SCALAR ) \
1056  CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR_OP( SCALAR, KokkosMultiply )
1057 
1058 #define CRSMATRIX_UQ_PCE_TESTS_STORAGE( STORAGE ) \
1059  typedef Sacado::UQ::PCE<STORAGE> UQ_PCE_ ## STORAGE; \
1060  CRSMATRIX_UQ_PCE_TESTS_MATRIXSCALAR( UQ_PCE_ ## STORAGE ) \
1061  CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR( UQ_PCE_ ## STORAGE ) \
1062  CRSMATRIX_UQ_PCE_MEAN_MULTIPLY_TESTS( UQ_PCE_ ## STORAGE )
1063 
1064 #define CRSMATRIX_UQ_PCE_TESTS_ORDINAL_SCALAR_DEVICE( ORDINAL, SCALAR, DEVICE ) \
1065  typedef Stokhos::DynamicStorage<ORDINAL,SCALAR,DEVICE> DS; \
1066  CRSMATRIX_UQ_PCE_TESTS_STORAGE( DS )
1067 
1068 #define CRSMATRIX_UQ_PCE_TESTS_DEVICE( DEVICE ) \
1069  CRSMATRIX_UQ_PCE_TESTS_ORDINAL_SCALAR_DEVICE( int, double, DEVICE )
Stokhos::StandardStorage< int, double > storage_type
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)
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
void operator()(const Matrix &A, const InputVector &x, OutputVector &y) const
expr val()
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
Kokkos::DefaultExecutionSpace execution_space
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
static void apply(const MatrixType matrix)
MatrixType buildDiagonalMatrix(typename MatrixType::ordinal_type nrow, typename MatrixType::ordinal_type pce_size, const CijkType &cijk)
bool compareRank2(const vector_type &y, const vector_type &y_exp, const scalar_type rel_tol, const scalar_type abs_tol, Teuchos::FancyOStream &out)
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(Kokkos_CrsMatrix_PCE, ReplaceValues, MatrixScalar)
Kokkos_MV_Multiply_Op KokkosMultiply
void multiply(const CrsMatrix< MatrixValue, Device, Layout > &A, const InputMultiVectorType &x, OutputMultiVectorType &y, const std::vector< OrdinalType > &col_indices, SingleColumnMultivectorMultiply)
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
static void apply(const MatrixType matrix)
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< View< T, P... > >::value, unsigned >::type dimension_scalar(const View< T, P... > &view)
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)
bool test_embedded_pce(const typename PCEType::ordinal_type nGrid, const typename PCEType::ordinal_type stoch_dim, const typename PCEType::ordinal_type poly_ord, Kokkos::DeviceConfig dev_config, Multiply multiply_op, Teuchos::FancyOStream &out)
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr expr expr expr j
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
kokkos_cijk_type build_cijk(ordinal_type stoch_dim, ordinal_type poly_ord)
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)
Legendre polynomial basis.
TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(Kokkos_CrsMatrix_PCE, Multiply, Scalar, MultiplyOp)
Abstract base class for 1-D orthogonal polynomials.
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
pce_type Scalar
bool compareRank1(const vector_type &y, const vector_type &y_exp, const scalar_type rel_tol, const scalar_type abs_tol, Teuchos::FancyOStream &out)
static void apply(const MatrixType matrix)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267
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)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
void operator()(const Matrix &A, const InputVector &x, OutputVector &y) const