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 "KokkosSparse_CrsMatrix.hpp"
47 #include "KokkosSparse_spmv.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>
135 {
136  using Teuchos::RCP;
137  using Teuchos::rcp;
138  using Teuchos::Array;
139 
140  typedef typename kokkos_cijk_type::value_type value_type;
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
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>
315  typedef typename MatrixType::size_type size_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  const ordinal_type pce_size = Kokkos::dimension_scalar(host_matrix_values);
346  bool success = true;
347  value_type val_expected(Kokkos::cijk(matrix.values));
348  for (ordinal_type row=0; row<nrow; ++row) {
349  val_expected = row;
350  bool s = compareVecs(host_matrix_values(row),
351  "matrix_values(row)",
352  val_expected,
353  "val_expected",
354  0.0, 0.0, out);
355  success = success && s;
356  }
357  return success;
358  }
359 };
360 
361 // Kernel to add values to the diagonal of a matrix
362 template <typename MatrixType>
365  typedef typename MatrixType::size_type size_type;
368 
369  const MatrixType m_matrix;
370  AddDiagonalValuesKernel(const MatrixType matrix) : m_matrix(matrix) {};
371 
372  // Replace diagonal entry for row 'i' with a value
373  KOKKOS_INLINE_FUNCTION
374  void operator() (const size_type i) const {
375  const ordinal_type row = i;
376  const ordinal_type col = i;
377  value_type val = value_type(row);
378  m_matrix.sumIntoValues(row, &col, 1, &val, false, true);
379  }
380 
381  // Kernel launch
382  static void apply(const MatrixType matrix) {
383  const size_type nrow = matrix.numRows();
384  Kokkos::parallel_for( nrow, AddDiagonalValuesKernel(matrix) );
385  }
386 
387  // Check the result is as expected
388  static bool check(const MatrixType matrix,
389  Teuchos::FancyOStream& out) {
390  typedef typename MatrixType::values_type matrix_values_type;
391  typename matrix_values_type::HostMirror host_matrix_values =
392  Kokkos::create_mirror_view(matrix.values);
393  Kokkos::deep_copy(host_matrix_values, matrix.values);
394  const ordinal_type nrow = matrix.numRows();
395  const ordinal_type pce_size = Kokkos::dimension_scalar(host_matrix_values);
396  bool success = true;
397  value_type val_expected(Kokkos::cijk(matrix.values));
398  for (ordinal_type row=0; row<nrow; ++row) {
399  val_expected = row;
400  bool s = compareVecs(host_matrix_values(row),
401  "matrix_values(row)",
402  val_expected,
403  "val_expected",
404  0.0, 0.0, out);
405  success = success && s;
406  }
407  return success;
408  }
409 };
410 
411 // Kernel to add values to the diagonal of a matrix where each thread
412 // adds to the same row (checks atomic really works)
413 template <typename MatrixType>
416  typedef typename MatrixType::size_type size_type;
419 
420  const MatrixType m_matrix;
421  AddDiagonalValuesAtomicKernel(const MatrixType matrix) : m_matrix(matrix) {};
422 
423  // Replace diagonal entry for row 'i' with a value
424  KOKKOS_INLINE_FUNCTION
425  void operator() (const size_type i) const {
426  const ordinal_type row = 0;
427  const ordinal_type col = 0;
429  m_matrix.sumIntoValues(row, &col, 1, &val, false, true);
430  }
431 
432  // Kernel launch
433  static void apply(const MatrixType matrix) {
434  const size_type nrow = matrix.numRows();
435  Kokkos::parallel_for( nrow, AddDiagonalValuesAtomicKernel(matrix) );
436  }
437 
438  // Check the result is as expected
439  static bool check(const MatrixType matrix,
440  Teuchos::FancyOStream& out) {
441  typedef typename MatrixType::values_type matrix_values_type;
442  typename matrix_values_type::HostMirror host_matrix_values =
443  Kokkos::create_mirror_view(matrix.values);
444  Kokkos::deep_copy(host_matrix_values, matrix.values);
445  const ordinal_type nrow = matrix.numRows();
446  const ordinal_type pce_size = Kokkos::dimension_scalar(host_matrix_values);
447  bool success = true;
448  value_type val_expected(Kokkos::cijk(matrix.values));
449  for (ordinal_type row=0; row<nrow; ++row) {
450  value_type val;
451  if (row == 0)
452  val_expected = nrow*(nrow-1)/2;
453  else
454  val_expected = 0.0;
455  bool s = compareVecs(host_matrix_values(row),
456  "matrix_values(row)",
457  val_expected,
458  "val_expected",
459  0.0, 0.0, out);
460  success = success && s;
461  }
462  return success;
463  }
464 };
465 
467  Kokkos_CrsMatrix_PCE, ReplaceValues, MatrixScalar )
468 {
469  typedef typename MatrixScalar::ordinal_type Ordinal;
470  typedef typename MatrixScalar::execution_space Device;
471  typedef typename MatrixScalar::cijk_type Cijk;
472  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
473 
474  // Build Cijk tensor
475  const Ordinal stoch_dim = 2;
476  const Ordinal poly_ord = 3;
477  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
478 
479  // Build diagonal matrix
480  const Ordinal nrow = 10;
481  const Ordinal pce_size = cijk.dimension();
482  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
483 
484  // Launch our kernel
486  kernel::apply(matrix);
487 
488  // Check the result
489  success = kernel::check(matrix, out);
490 }
491 
493  Kokkos_CrsMatrix_PCE, SumIntoValues, MatrixScalar )
494 {
495  typedef typename MatrixScalar::ordinal_type Ordinal;
496  typedef typename MatrixScalar::execution_space Device;
497  typedef typename MatrixScalar::cijk_type Cijk;
498  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
499 
500  // Build Cijk tensor
501  const Ordinal stoch_dim = 2;
502  const Ordinal poly_ord = 3;
503  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
504 
505  // Build diagonal matrix
506  const Ordinal nrow = 10;
507  const Ordinal pce_size = cijk.dimension();
508  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
509 
510  // Launch our kernel
511  typedef AddDiagonalValuesKernel<Matrix> kernel;
512  kernel::apply(matrix);
513 
514  // Check the result
515  success = kernel::check(matrix, out);
516 }
517 
519  Kokkos_CrsMatrix_PCE, SumIntoValuesAtomic, MatrixScalar )
520 {
521  typedef typename MatrixScalar::ordinal_type Ordinal;
522  typedef typename MatrixScalar::execution_space Device;
523  typedef typename MatrixScalar::cijk_type Cijk;
524  typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
525 
526  // Build Cijk tensor
527  const Ordinal stoch_dim = 2;
528  const Ordinal poly_ord = 3;
529  Cijk cijk = build_cijk<Cijk>(stoch_dim, poly_ord);
530 
531  // Build diagonal matrix
532  const Ordinal nrow = 10;
533  const Ordinal pce_size = cijk.dimension();
534  Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
535 
536  // Launch our kernel
538  kernel::apply(matrix);
539 
540  // Check the result
541  success = kernel::check(matrix, out);
542 }
543 
544 template <typename PCEType, typename Multiply>
545 bool test_embedded_pce(const typename PCEType::ordinal_type nGrid,
546  const typename PCEType::ordinal_type stoch_dim,
547  const typename PCEType::ordinal_type poly_ord,
548  KokkosSparse::DeviceConfig dev_config,
549  Multiply multiply_op,
550  Teuchos::FancyOStream& out)
551 {
552  typedef typename PCEType::ordinal_type ordinal_type;
553  typedef typename PCEType::value_type scalar_type;
554  typedef typename PCEType::storage_type storage_type;
555  typedef typename PCEType::cijk_type cijk_type;
557  typedef Kokkos::LayoutLeft Layout;
558  typedef Kokkos::View< PCEType*, Layout, execution_space > block_vector_type;
559  typedef KokkosSparse::CrsMatrix< PCEType, ordinal_type, execution_space > block_matrix_type;
560  typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
561  typedef typename block_matrix_type::values_type matrix_values_type;
562 
563  // Build Cijk tensor
564  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
565  const ordinal_type stoch_length = cijk.dimension();
566  // const ordinal_type align = 8;
567  // const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
568  const ordinal_type stoch_length_aligned = stoch_length;
569 
570  // Check pce_length == storage_type::static_size for static storage
571  TEUCHOS_TEST_FOR_EXCEPTION(
572  storage_type::is_static && storage_type::static_size != stoch_length,
573  std::logic_error,
574  "Static storage size must equal pce size");
575 
576  // Generate FEM graph:
577  const ordinal_type fem_length = nGrid * nGrid * nGrid;
578  std::vector< std::vector<ordinal_type> > fem_graph;
579  const ordinal_type fem_graph_length = generate_fem_graph( nGrid, fem_graph );
580 
581  //------------------------------
582  // Generate input/output multivectors -- Sacado dimension is always last,
583  // regardless of LayoutLeft/Right
584 
585  block_vector_type x =
586  Kokkos::make_view<block_vector_type>("x", cijk, fem_length, stoch_length_aligned);
587  block_vector_type y =
588  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, stoch_length_aligned);
589 
590  typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
591  typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
592 
593  // View the block vector as an array of the embedded intrinsic type.
594  typename block_vector_type::HostMirror::array_type hax = hx ;
595  typename block_vector_type::HostMirror::array_type hay = hy ;
596 
597  for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
598  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
599  hax(iRowStoch, iRowFEM) =
600  generate_vector_coefficient<scalar_type>(
601  fem_length, stoch_length, iRowFEM, iRowStoch );
602  hay(iRowStoch, iRowFEM) = 0.0;
603  }
604  }
605 
606  Kokkos::deep_copy( x, hx );
607  Kokkos::deep_copy( y, hy );
608 
609  //------------------------------
610  // Generate block matrix -- it is always LayoutRight (currently)
611 
612  matrix_graph_type matrix_graph =
613  Kokkos::create_staticcrsgraph<matrix_graph_type>(
614  std::string("test crs graph"), fem_graph);
615  matrix_values_type matrix_values =
616  Kokkos::make_view<matrix_values_type>(
617  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, stoch_length_aligned);
618  block_matrix_type matrix(
619  "block_matrix", fem_length, matrix_values, matrix_graph);
620  matrix.dev_config = dev_config;
621 
622  typename matrix_values_type::HostMirror hM =
623  Kokkos::create_mirror_view( matrix.values );
624 
625  typename matrix_values_type::HostMirror::array_type haM = hM ;
626 
627  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
628  const ordinal_type row_size = fem_graph[iRowFEM].size();
629  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
630  ++iRowEntryFEM, ++iEntryFEM) {
631  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
632 
633  for (ordinal_type k=0; k<stoch_length; ++k) {
634  haM(iEntryFEM, k) =
635  generate_matrix_coefficient<scalar_type>(
636  fem_length, stoch_length, iRowFEM, iColFEM, k);
637  }
638  }
639  }
640 
641  Kokkos::deep_copy( matrix.values, hM );
642 
643  //------------------------------
644  // multiply
645 
646  multiply_op( matrix, x, y );
647 
648  //------------------------------
649  // generate correct answer
650 
651  typedef typename block_vector_type::array_type array_type;
652  array_type ay_expected =
653  array_type("ay_expected", stoch_length_aligned, fem_length);
654  typename array_type::HostMirror hay_expected =
655  Kokkos::create_mirror_view(ay_expected);
656  typename cijk_type::HostMirror host_cijk =
658  Kokkos::deep_copy(host_cijk, cijk);
659  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
660  const ordinal_type row_size = fem_graph[iRowFEM].size();
661  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
662  ++iRowEntryFEM, ++iEntryFEM) {
663  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
664  for (ordinal_type i=0; i<stoch_length; ++i) {
665  const ordinal_type num_entry = host_cijk.num_entry(i);
666  const ordinal_type entry_beg = host_cijk.entry_begin(i);
667  const ordinal_type entry_end = entry_beg + num_entry;
668  scalar_type tmp = 0;
669  for (ordinal_type entry = entry_beg; entry < entry_end; ++entry) {
670  const ordinal_type j = host_cijk.coord(entry,0);
671  const ordinal_type k = host_cijk.coord(entry,1);
672  const scalar_type a_j =
673  generate_matrix_coefficient<scalar_type>(
674  fem_length, stoch_length, iRowFEM, iColFEM, j);
675  const scalar_type a_k =
676  generate_matrix_coefficient<scalar_type>(
677  fem_length, stoch_length, iRowFEM, iColFEM, k);
678  const scalar_type x_j =
679  generate_vector_coefficient<scalar_type>(
680  fem_length, stoch_length, iColFEM, j);
681  const scalar_type x_k =
682  generate_vector_coefficient<scalar_type>(
683  fem_length, stoch_length, iColFEM, k);
684  tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
685  }
686  hay_expected(i, iRowFEM) += tmp;
687  }
688  }
689  }
690  Kokkos::deep_copy( ay_expected, hay_expected );
691 
692  //------------------------------
693  // check
694 
695  typename block_vector_type::array_type ay = y;
698  bool success = compare_rank_2_views(ay, ay_expected, rel_tol, abs_tol, out);
699 
700  return success;
701 }
702 
704  Kokkos_CrsMatrix_PCE, Multiply, Scalar, MultiplyOp )
705 {
706  typedef typename Scalar::ordinal_type Ordinal;
707 
708  const Ordinal nGrid = 5;
709  const Ordinal stoch_dim = 2;
710  const Ordinal poly_ord = 3;
711  KokkosSparse::DeviceConfig dev_config;
712 
713  success = test_embedded_pce<Scalar>(
714  nGrid, stoch_dim, poly_ord, dev_config, MultiplyOp(), out);
715 }
716 
717 struct Kokkos_MV_Multiply_Op {
718  template <typename Matrix, typename InputVector, typename OutputVector>
719  void operator() (const Matrix& A,
720  const InputVector& x,
721  OutputVector& y) const {
722  KokkosSparse::spmv("N", typename Matrix::value_type(1.0) , A, x, typename Matrix::value_type(0.0), y);
723  }
724 };
725 
726 template <typename Tag>
727 struct Stokhos_MV_Multiply_Op {
728  Tag tag;
729  Stokhos_MV_Multiply_Op(const Tag& tg = Tag()) : tag(tg) {}
730 
731  template <typename Matrix, typename InputVector, typename OutputVector>
732  void operator() (const Matrix& A,
733  const InputVector& x,
734  OutputVector& y) const {
735  Stokhos::multiply(A, x, y, tag);
736  }
737 };
738 
740  Kokkos_CrsMatrix_PCE, MeanMultiplyRank1, Scalar )
741 {
742  typedef typename Scalar::ordinal_type Ordinal;
743 
744  const Ordinal nGrid = 5;
745  const Ordinal stoch_dim = 5;
746  const Ordinal poly_ord = 3;
747  KokkosSparse::DeviceConfig dev_config;
748 
749  typedef typename Scalar::ordinal_type ordinal_type;
750  typedef typename Scalar::value_type scalar_type;
751  typedef typename Scalar::storage_type storage_type;
752  typedef typename Scalar::cijk_type cijk_type;
754  typedef Kokkos::LayoutLeft Layout;
755  typedef Kokkos::View< Scalar*, Layout, execution_space > block_vector_type;
756  typedef KokkosSparse::CrsMatrix< Scalar, ordinal_type, execution_space > block_matrix_type;
757  typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
758  typedef typename block_matrix_type::values_type matrix_values_type;
759 
760  // Build Cijk tensor
761  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
762  cijk_type mean_cijk =
763  Stokhos::create_mean_based_product_tensor<execution_space, ordinal_type, scalar_type>();
764  const ordinal_type stoch_length = cijk.dimension();
765  const ordinal_type align = 8;
766  const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
767 
768  // Generate FEM graph:
769  const ordinal_type fem_length = nGrid * nGrid * nGrid;
770  std::vector< std::vector<ordinal_type> > fem_graph;
771  const ordinal_type fem_graph_length = generate_fem_graph( nGrid, fem_graph );
772 
773  block_vector_type x =
774  Kokkos::make_view<block_vector_type>("x", cijk, fem_length, stoch_length_aligned);
775  block_vector_type y =
776  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, stoch_length_aligned);
777 
778  block_vector_type y_expected =
779  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, stoch_length_aligned);
780 
781  typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
782  typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
783  typename block_vector_type::HostMirror hy_expected =
784  Kokkos::create_mirror_view( y_expected );
785 
786  // View the block vector as an array of the embedded intrinsic type.
787  typename block_vector_type::HostMirror::array_type hax = hx ;
788  typename block_vector_type::HostMirror::array_type hay = hy ;
789  typename block_vector_type::HostMirror::array_type hay_expected =
790  hy_expected ;
791 
792  for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
793  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
794  hax(iRowStoch,iRowFEM) =
795  generate_vector_coefficient<scalar_type>(
796  fem_length, stoch_length, iRowFEM, iRowStoch );
797  hay(iRowStoch,iRowFEM) = 0.0;
798  hay_expected(iRowStoch,iRowFEM) = 0.0;
799  }
800  }
801  Kokkos::deep_copy( x, hx );
802  Kokkos::deep_copy( y, hy );
803  Kokkos::deep_copy( y_expected, hy_expected );
804 
805  //------------------------------
806  // Generate block matrix -- it is always LayoutRight (currently)
807  matrix_graph_type matrix_graph =
808  Kokkos::create_staticcrsgraph<matrix_graph_type>(
809  std::string("test crs graph"), fem_graph);
810  matrix_values_type matrix_values =
811  Kokkos::make_view<matrix_values_type>(
812  Kokkos::ViewAllocateWithoutInitializing("matrix"), mean_cijk, fem_graph_length, ordinal_type(1)); //instead of stoch_length
813  block_matrix_type matrix(
814  "block_matrix", fem_length, matrix_values, matrix_graph);
815  matrix.dev_config = dev_config;
816 
817  typename matrix_values_type::HostMirror hM =
818  Kokkos::create_mirror_view( matrix.values );
819  typename matrix_values_type::HostMirror::array_type haM = hM ;
820 
821  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
822  const ordinal_type row_size = fem_graph[iRowFEM].size();
823  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
824  ++iRowEntryFEM, ++iEntryFEM) {
825  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
826 
827  haM(iEntryFEM, 0) =
828  generate_matrix_coefficient<scalar_type>(
829  fem_length, 1, iRowFEM, iColFEM, 0);
830  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
831  hay_expected(iRowStoch,iRowFEM) +=
832  haM(iEntryFEM, 0) * hax(iRowStoch,iColFEM);
833  }
834  }
835  }
836  Kokkos::deep_copy( matrix.values, hM );
837  Kokkos::deep_copy( y_expected, hy_expected );
838 
839  /*
840  //Generate same matrix with stochastic dim = Kokkos::dimension_scalar(x) (i.e. not = 1)
841  matrix_values_type full_matrix_values =
842  Kokkos::make_view<matrix_values_type>(
843  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, stoch_length_aligned);
844  block_matrix_type full_matrix(
845  "block_matrix", fem_length, full_matrix_values, matrix_graph);
846  matrix.dev_config = dev_config;
847 
848  typename matrix_values_type::HostMirror full_hM =
849  Kokkos::create_mirror_view( full_matrix.values );
850  typename matrix_values_type::HostMirror::array_type full_haM = full_hM ;
851 
852  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
853  const ordinal_type row_size = fem_graph[iRowFEM].size();
854  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
855  ++iRowEntryFEM, ++iEntryFEM) {
856  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
857 
858  for (ordinal_type k=0; k<stoch_length; ++k) {
859  if (k == 0)
860  full_haM(iEntryFEM, k) =
861  generate_matrix_coefficient<scalar_type>(
862  fem_length, 1, iRowFEM, iColFEM, k);
863  else
864  full_haM(iEntryFEM, k) = 0.0;
865  }
866  }
867  }
868 
869  Kokkos::deep_copy( full_matrix.values, full_hM );
870  */
871 
872  //------------------------------
873  // multiply
874 
875  KokkosSparse::spmv("N", Scalar(1.0) , matrix, x, Scalar(0.0), y);
876 
877  //------------------------------
878  // multiply with same matrix but with sacado_size = x.sacado_size
879 
880  //Kokkos::MV_Multiply( y_expected, full_matrix, x );
881 
882  //------------------------------
883  // check
886  success = compareRank1(y, y_expected, rel_tol, abs_tol, out);
887 }
888 
889 
891  Kokkos_CrsMatrix_PCE, MeanMultiplyRank2, Scalar )
892 {
893  typedef typename Scalar::ordinal_type ordinal_type;
894  typedef typename Scalar::value_type scalar_type;
895  typedef typename Scalar::storage_type storage_type;
896  typedef typename Scalar::cijk_type cijk_type;
898  typedef Kokkos::LayoutLeft Layout;
899  typedef Kokkos::View< Scalar**, Layout, execution_space > block_vector_type;
900  typedef KokkosSparse::CrsMatrix< Scalar, ordinal_type, execution_space > block_matrix_type;
901  typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
902  typedef typename block_matrix_type::values_type matrix_values_type;
903 
904 
905  const ordinal_type nGrid = 5;
906  const ordinal_type stoch_dim = 2;
907  const ordinal_type poly_ord = 3;
908  KokkosSparse::DeviceConfig dev_config;
909 
910  // Build Cijk tensor
911  cijk_type cijk = build_cijk<cijk_type>(stoch_dim, poly_ord);
912  cijk_type mean_cijk =
913  Stokhos::create_mean_based_product_tensor<execution_space, ordinal_type, scalar_type>();
914  const ordinal_type stoch_length = cijk.dimension();
915  const ordinal_type align = 8;
916  const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
917  const ordinal_type num_cols = 2;
918  // Generate FEM graph:
919  const ordinal_type fem_length = nGrid * nGrid * nGrid;
920  std::vector< std::vector<ordinal_type> > fem_graph;
921  const ordinal_type fem_graph_length = generate_fem_graph( nGrid, fem_graph );
922 
923  block_vector_type x =
924  Kokkos::make_view<block_vector_type>("x", cijk, fem_length, num_cols, stoch_length_aligned);
925  block_vector_type y =
926  Kokkos::make_view<block_vector_type>("y", cijk, fem_length, num_cols, stoch_length_aligned);
927 
928  block_vector_type y_expected =
929  Kokkos::make_view<block_vector_type>("y_expected", cijk, fem_length, num_cols, stoch_length_aligned);
930 
931  typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
932  typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
933  typename block_vector_type::HostMirror hy_expected =
934  Kokkos::create_mirror_view( y_expected );
935 
936  for (ordinal_type i=0; i<num_cols; ++i){
937  for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
938  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
939  hx(iRowFEM,i).fastAccessCoeff(iRowStoch) =
940  generate_vector_coefficient<scalar_type>(
941  fem_length, stoch_length, iRowFEM, iRowStoch );
942  hy(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
943  hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
944  }
945  }
946  }
947  Kokkos::deep_copy( x, hx );
948  Kokkos::deep_copy( y, hy );
949  Kokkos::deep_copy( y_expected, hy_expected );
950 
951  //------------------------------
952  // Generate matrix with stochastic dimension 1
953  matrix_graph_type matrix_graph =
954  Kokkos::create_staticcrsgraph<matrix_graph_type>(
955  std::string("test crs graph"), fem_graph);
956  matrix_values_type matrix_values =
957  Kokkos::make_view<matrix_values_type>(
958  Kokkos::ViewAllocateWithoutInitializing("matrix"), mean_cijk, fem_graph_length, ordinal_type(1));
959  block_matrix_type matrix(
960  "block_matrix", fem_length, matrix_values, matrix_graph);
961  matrix.dev_config = dev_config;
962 
963  typename matrix_values_type::HostMirror hM =
964  Kokkos::create_mirror_view( matrix.values );
965 
966  typename matrix_values_type::HostMirror::array_type haM = hM ;
967 
968  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
969  const ordinal_type row_size = fem_graph[iRowFEM].size();
970  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
971  ++iRowEntryFEM, ++iEntryFEM) {
972  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
973 
974  haM(iEntryFEM, 0) =
975  generate_matrix_coefficient<scalar_type>(
976  fem_length, 1, iRowFEM, iColFEM, 0);
977  for (ordinal_type i=0; i<num_cols; ++i){
978  for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
979  hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) +=
980  haM(iEntryFEM, 0) * hx(iColFEM,i).fastAccessCoeff(iRowStoch);
981  }
982  }
983  }
984  }
985 
986  Kokkos::deep_copy( matrix.values, hM );
987  Kokkos::deep_copy( y_expected, hy_expected );
988 
989  /*
990  //Generate same matrix with stochastic dim = Kokkos::dimension_scalar(x) (i.e. not = 1)
991  matrix_values_type full_matrix_values =
992  Kokkos::make_view<matrix_values_type>(
993  Kokkos::ViewAllocateWithoutInitializing("matrix"), cijk, fem_graph_length, stoch_length_aligned);
994  block_matrix_type full_matrix(
995  "block_matrix", fem_length, full_matrix_values, matrix_graph);
996  matrix.dev_config = dev_config;
997 
998  typename matrix_values_type::HostMirror full_hM =
999  Kokkos::create_mirror_view( full_matrix.values );
1000 
1001  typename matrix_values_type::HostMirror::array_type full_haM = full_hM ;
1002 
1003  for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
1004  const ordinal_type row_size = fem_graph[iRowFEM].size();
1005  for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
1006  ++iRowEntryFEM, ++iEntryFEM) {
1007  const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
1008 
1009  for (ordinal_type k=0; k<stoch_length; ++k) {
1010  if (k == 0)
1011  full_haM(iEntryFEM, k) =
1012  generate_matrix_coefficient<scalar_type>(
1013  fem_length, 1, iRowFEM, iColFEM, k);
1014  else
1015  full_haM(iEntryFEM, k) = 0.0;
1016  }
1017  }
1018  }
1019 
1020  Kokkos::deep_copy( full_matrix.values, full_hM );
1021  */
1022 
1023  //------------------------------
1024  // multiply
1025 
1026  KokkosSparse::spmv("N", Scalar(1.0) , matrix, x, Scalar(0.0), y);
1027 
1028  //------------------------------
1029  // multiply with full matrix
1030 
1031  //Kokkos::MV_Multiply( y_expected, full_matrix, x );
1032 
1033  //------------------------------
1034  // check
1035 
1038  success = compareRank2(y, y_expected, rel_tol, abs_tol, out);
1039 }
1040 
1041 
1043 
1044 #define CRSMATRIX_UQ_PCE_TESTS_MATRIXSCALAR( SCALAR ) \
1045  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1046  Kokkos_CrsMatrix_PCE, ReplaceValues, SCALAR ) \
1047  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1048  Kokkos_CrsMatrix_PCE, SumIntoValues, SCALAR ) \
1049  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1050  Kokkos_CrsMatrix_PCE, SumIntoValuesAtomic, SCALAR )
1051 #define CRSMATRIX_UQ_PCE_MEAN_MULTIPLY_TESTS( SCALAR ) \
1052  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1053  Kokkos_CrsMatrix_PCE, MeanMultiplyRank1, SCALAR ) \
1054  TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1055  Kokkos_CrsMatrix_PCE, MeanMultiplyRank2, SCALAR )
1056 #define CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR_OP( SCALAR, OP ) \
1057  TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( \
1058  Kokkos_CrsMatrix_PCE, Multiply, SCALAR, OP )
1059 
1060 #define CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR( SCALAR ) \
1061  CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR_OP( SCALAR, KokkosMultiply )
1062 
1063 #define CRSMATRIX_UQ_PCE_TESTS_STORAGE( STORAGE ) \
1064  typedef Sacado::UQ::PCE<STORAGE> UQ_PCE_ ## STORAGE; \
1065  CRSMATRIX_UQ_PCE_TESTS_MATRIXSCALAR( UQ_PCE_ ## STORAGE ) \
1066  CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR( UQ_PCE_ ## STORAGE ) \
1067  CRSMATRIX_UQ_PCE_MEAN_MULTIPLY_TESTS( UQ_PCE_ ## STORAGE )
1068 
1069 #define CRSMATRIX_UQ_PCE_TESTS_ORDINAL_SCALAR_DEVICE( ORDINAL, SCALAR, DEVICE ) \
1070  typedef Stokhos::DynamicStorage<ORDINAL,SCALAR,DEVICE> DS; \
1071  CRSMATRIX_UQ_PCE_TESTS_STORAGE( DS )
1072 
1073 #define CRSMATRIX_UQ_PCE_TESTS_DEVICE( DEVICE ) \
1074  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.
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)
void operator()(const Matrix &A, const InputVector &x, OutputVector &y) const
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
Kokkos::DefaultExecutionSpace execution_space
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j) - expr2.val(j)
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)
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)
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)
expr val()
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