42#include "Teuchos_UnitTestHelpers.hpp"
46#include "Stokhos_Ensemble_Sizes.hpp"
51#include "Kokkos_Core.hpp"
54template<
typename IntType >
61 return k + N * (
j + N * i );
64template <
typename ordinal >
67 std::vector< std::vector<ordinal> >& graph )
69 graph.resize( N * N * N, std::vector<ordinal>() );
73 for (
int i = 0; i < (int) N; ++i ) {
74 for (
int j = 0;
j < (int) N; ++
j ) {
75 for (
int k = 0; k < (int) N; ++k ) {
79 graph[row].reserve(27);
81 for (
int ii = -1; ii < 2; ++ii ) {
82 for (
int jj = -1; jj < 2; ++jj ) {
83 for (
int kk = -1; kk < 2; ++kk ) {
84 if ( 0 <= i + ii && i + ii < (
int) N &&
85 0 <=
j + jj &&
j + jj < (
int) N &&
86 0 <= k + kk && k + kk < (
int) N ) {
89 graph[row].push_back(col);
92 total += graph[row].size();
98template <
typename scalar,
typename ordinal>
101 const ordinal nStoch,
102 const ordinal iRowFEM,
103 const ordinal iColFEM,
104 const ordinal iStoch )
106 const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
107 ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
109 const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
111 return A_fem + A_stoch;
115template <
typename scalar,
typename ordinal>
118 const ordinal nStoch,
119 const ordinal iColFEM,
120 const ordinal iStoch )
122 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
123 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
124 return X_fem + X_stoch;
130template <>
struct ScalarTol<float> {
static float tol() {
return 1e-4; } };
134template <
typename array_type,
typename scalar_type>
136 const array_type& y_exp,
137 const scalar_type rel_tol,
138 const scalar_type abs_tol,
139 Teuchos::FancyOStream& out)
141 typedef typename array_type::size_type size_type;
142 typename array_type::HostMirror hy = Kokkos::create_mirror_view(y);
143 typename array_type::HostMirror hy_exp = Kokkos::create_mirror_view(y_exp);
144 Kokkos::deep_copy(hy, y);
145 Kokkos::deep_copy(hy_exp, y_exp);
147 size_type num_rows = y.extent(0);
148 size_type num_cols = y.extent(1);
150 for (size_type i=0; i<num_rows; ++i) {
151 for (size_type
j=0;
j<num_cols; ++
j) {
152 scalar_type diff = std::abs( hy(i,
j) - hy_exp(i,
j) );
153 scalar_type tol = rel_tol*std::abs(hy_exp(i,
j)) + abs_tol;
155 out <<
"y_expected(" << i <<
"," <<
j <<
") - "
156 <<
"y(" << i <<
"," <<
j <<
") = " << hy_exp(i,
j)
157 <<
" - " << hy(i,
j) <<
" == "
158 << diff <<
" < " << tol <<
" : ";
164 success = success && s;
172template <
typename MatrixType>
175 typename MatrixType::ordinal_type mp_vector_size) {
176 typedef typename MatrixType::ordinal_type ordinal_type;
177 typedef typename MatrixType::StaticCrsGraphType matrix_graph_type;
178 typedef typename MatrixType::values_type matrix_values_type;
180 std::vector< std::vector<ordinal_type> > graph(nrow);
181 for (ordinal_type i=0; i<nrow; ++i)
182 graph[i] = std::vector<ordinal_type>(1, i);
183 ordinal_type graph_length = nrow;
185 matrix_graph_type matrix_graph =
186 Kokkos::create_staticcrsgraph<matrix_graph_type>(
"graph", graph);
187 matrix_values_type matrix_values =
188 matrix_values_type(
"values", graph_length, mp_vector_size);
190 MatrixType matrix(
"matrix", nrow, matrix_values, matrix_graph);
199template <
typename MatrixType>
210 KOKKOS_INLINE_FUNCTION
215 m_matrix.replaceValues(row, &col, 1, &
val,
false,
true);
219 static void apply(
const MatrixType matrix) {
225 static bool check(
const MatrixType matrix,
226 Teuchos::FancyOStream& out) {
227 typedef typename MatrixType::values_type matrix_values_type;
228 typename matrix_values_type::HostMirror host_matrix_values =
229 Kokkos::create_mirror_view(matrix.values);
230 Kokkos::deep_copy(host_matrix_values, matrix.values);
235 bool s = compareVecs(host_matrix_values(row),
236 "matrix_values(row)",
240 success = success && s;
247template <
typename MatrixType>
258 KOKKOS_INLINE_FUNCTION
263 m_matrix.sumIntoValues(row, &col, 1, &
val,
false,
true);
267 static void apply(
const MatrixType matrix) {
273 static bool check(
const MatrixType matrix,
274 Teuchos::FancyOStream& out) {
275 typedef typename MatrixType::values_type matrix_values_type;
276 typename matrix_values_type::HostMirror host_matrix_values =
277 Kokkos::create_mirror_view(matrix.values);
278 Kokkos::deep_copy(host_matrix_values, matrix.values);
283 bool s = compareVecs(host_matrix_values(row),
284 "matrix_values(row)",
288 success = success && s;
296template <
typename MatrixType>
307 KOKKOS_INLINE_FUNCTION
312 m_matrix.sumIntoValues(row, &col, 1, &
val,
false,
true);
316 static void apply(
const MatrixType matrix) {
322 static bool check(
const MatrixType matrix,
323 Teuchos::FancyOStream& out) {
324 typedef typename MatrixType::values_type matrix_values_type;
325 typename matrix_values_type::HostMirror host_matrix_values =
326 Kokkos::create_mirror_view(matrix.values);
327 Kokkos::deep_copy(host_matrix_values, matrix.values);
337 bool s = compareVecs(host_matrix_values(row),
338 "matrix_values(row)",
342 success = success && s;
351 Kokkos_CrsMatrix_MP, ReplaceValues, MatrixScalar )
353 typedef typename MatrixScalar::ordinal_type
Ordinal;
355 typedef Kokkos::Device<execution_space, typename execution_space::memory_space> Device;
356 typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
360 Matrix matrix = buildDiagonalMatrix<Matrix>(nrow,
VectorSize);
364 kernel::apply(matrix);
367 success = kernel::check(matrix, out);
371 Kokkos_CrsMatrix_MP, SumIntoValues, MatrixScalar )
373 typedef typename MatrixScalar::ordinal_type
Ordinal;
375 typedef Kokkos::Device<execution_space, typename execution_space::memory_space> Device;
376 typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
380 Matrix matrix = buildDiagonalMatrix<Matrix>(nrow,
VectorSize);
384 kernel::apply(matrix);
387 success = kernel::check(matrix, out);
391 Kokkos_CrsMatrix_MP, SumIntoValuesAtomic, MatrixScalar )
393 typedef typename MatrixScalar::ordinal_type
Ordinal;
395 typedef Kokkos::Device<execution_space, typename execution_space::memory_space> Device;
396 typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
400 Matrix matrix = buildDiagonalMatrix<Matrix>(nrow,
VectorSize);
404 kernel::apply(matrix);
407 success = kernel::check(matrix, out);
414template <
class ViewType,
class OrdinalType,
size_t I >
416 static ViewType
create_view(
const std::string & name,
const OrdinalType & fem_length,
const OrdinalType & stoch_length ) {
417 return ViewType(name, fem_length, stoch_length);
421template <
class ViewType,
class OrdinalType>
423 static ViewType
create_view(
const std::string & name,
const OrdinalType & fem_length,
const OrdinalType & stoch_length ) {
425 return ViewType(name, fem_length);
429template <
class ViewType,
class OrdinalType>
431 static ViewType
create_view(
const std::string & name,
const OrdinalType & fem_length,
const OrdinalType & stoch_length ) {
434 return ViewType(name);
439template <
typename VectorType,
typename Multiply>
441 const typename VectorType::ordinal_type stoch_length,
442 KokkosSparse::DeviceConfig dev_config,
443 Multiply multiply_op,
444 Teuchos::FancyOStream& out)
446 typedef typename VectorType::ordinal_type ordinal_type;
447 typedef typename VectorType::value_type scalar_type;
450 typedef Kokkos::Device<execution_space, typename execution_space::memory_space> device_type;
451 typedef Kokkos::LayoutRight Layout;
452 typedef Kokkos::View< VectorType*, Layout, execution_space > block_vector_type;
453 typedef KokkosSparse::CrsMatrix< VectorType, ordinal_type, device_type > block_matrix_type;
454 typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
455 typedef typename block_matrix_type::values_type matrix_values_type;
458 TEUCHOS_TEST_FOR_EXCEPTION(
459 storage_type::is_static && storage_type::static_size != stoch_length,
461 "Static storage size must equal ensemble size");
464 ordinal_type fem_length = nGrid * nGrid * nGrid;
465 std::vector< std::vector<ordinal_type> > fem_graph;
478 block_vector_type x =
479 block_vector_type(
"x", fem_length, stoch_length);
480 block_vector_type y =
481 block_vector_type(
"y", fem_length, stoch_length);
483 typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
484 typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
487 typename block_vector_type::HostMirror::array_type hax = hx ;
488 typename block_vector_type::HostMirror::array_type hay = hy ;
490 for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
491 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
492 hax(iRowFEM,iRowStoch) =
493 generate_vector_coefficient<scalar_type>(
494 fem_length, stoch_length, iRowFEM, iRowStoch );
495 hay(iRowFEM,iRowStoch) = 0.0;
499 Kokkos::deep_copy( x, hx );
500 Kokkos::deep_copy( y, hy );
505 matrix_graph_type matrix_graph =
506 Kokkos::create_staticcrsgraph<matrix_graph_type>(
507 std::string(
"test crs graph"), fem_graph);
512 matrix_values_type matrix_values =
513 matrix_values_type(
"matrix", fem_graph_length, stoch_length);
514 block_matrix_type matrix(
515 "block_matrix", fem_length, matrix_values, matrix_graph);
516 matrix.dev_config = dev_config;
518 typename matrix_values_type::HostMirror hM =
519 Kokkos::create_mirror_view( matrix.values );
521 typename matrix_values_type::HostMirror::array_type haM = hM ;
523 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
524 const ordinal_type row_size = fem_graph[iRowFEM].size();
525 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
526 ++iRowEntryFEM, ++iEntryFEM) {
527 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
529 for (ordinal_type k=0; k<stoch_length; ++k) {
531 generate_matrix_coefficient<scalar_type>(
532 fem_length, stoch_length, iRowFEM, iColFEM, k);
537 Kokkos::deep_copy( matrix.values, hM );
542 multiply_op( matrix, x, y );
547 typedef typename block_vector_type::array_type array_type;
548#ifdef KOKKOS_ENABLE_DEPRECATED_CODE
549 array_type ay_expected =
550 array_type(
"ay_expected", fem_length, stoch_length);
552 array_type ay_expected =
555 typename array_type::HostMirror hay_expected =
556 Kokkos::create_mirror_view(ay_expected);
557 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
558 const ordinal_type row_size = fem_graph[iRowFEM].size();
559 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
560 ++iRowEntryFEM, ++iEntryFEM) {
561 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
562 for (ordinal_type k=0; k<stoch_length; ++k) {
563 hay_expected(iRowFEM, k) +=
564 generate_matrix_coefficient<scalar_type>(
565 fem_length, stoch_length, iRowFEM, iColFEM, k) *
566 generate_vector_coefficient<scalar_type>(
567 fem_length, stoch_length, iColFEM, k );
571 Kokkos::deep_copy( ay_expected, hay_expected );
576 typename block_vector_type::array_type ay = y;
585 template <
typename Matrix,
typename InputVector,
typename OutputVector>
587 const InputVector& x,
588 OutputVector& y)
const {
589 KokkosSparse::spmv(
"N",
typename Matrix::value_type(1.0) , A, x,
typename Matrix::value_type(0.0), y);
593template <
typename Tag>
598 template <
typename Matrix,
typename InputVector,
typename OutputVector>
600 const InputVector& x,
601 OutputVector& y)
const {
609#define CRSMATRIX_MP_VECTOR_TESTS_MATRIXSCALAR( SCALAR ) \
610 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
611 Kokkos_CrsMatrix_MP, ReplaceValues, SCALAR ) \
612 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
613 Kokkos_CrsMatrix_MP, SumIntoValues, SCALAR ) \
614 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
615 Kokkos_CrsMatrix_MP, SumIntoValuesAtomic, SCALAR )
617#define CRSMATRIX_MP_VECTOR_TESTS_STORAGE( STORAGE ) \
618 typedef Sacado::MP::Vector<STORAGE> MP_Vector_ ## STORAGE; \
619 CRSMATRIX_MP_VECTOR_TESTS_MATRIXSCALAR( MP_Vector_ ## STORAGE )
621#define CRSMATRIX_MP_VECTOR_TESTS_ORDINAL_SCALAR_DEVICE( ORDINAL, SCALAR, DEVICE ) \
622 typedef Stokhos::StaticFixedStorage<ORDINAL,SCALAR,VectorSize,DEVICE> SFS; \
623 typedef Stokhos::DynamicStorage<ORDINAL,SCALAR,DEVICE> DS; \
624 CRSMATRIX_MP_VECTOR_TESTS_STORAGE( SFS ) \
625 CRSMATRIX_MP_VECTOR_TESTS_STORAGE( DS )
627#define CRSMATRIX_MP_VECTOR_TESTS_DEVICE( DEVICE ) \
628 CRSMATRIX_MP_VECTOR_TESTS_ORDINAL_SCALAR_DEVICE( int, double, DEVICE )
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
bool test_embedded_vector(const typename VectorType::ordinal_type nGrid, const typename VectorType::ordinal_type stoch_length, KokkosSparse::DeviceConfig dev_config, Multiply multiply_op, Teuchos::FancyOStream &out)
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)
MatrixType buildDiagonalMatrix(typename MatrixType::ordinal_type nrow, typename MatrixType::ordinal_type mp_vector_size)
Kokkos_MV_Multiply_Op KokkosMultiply
const unsigned VectorSize
TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(Kokkos_CrsMatrix_MP, ReplaceValues, MatrixScalar)
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
Stokhos_MV_Multiply_Op< Stokhos::DefaultMultiply > DefaultMultiply
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
Stokhos::StandardStorage< int, double > storage_type
Kokkos::DefaultExecutionSpace execution_space
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)
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< View< T, P... > >::value, unsigned >::type dimension_scalar(const View< T, P... > &view)
void multiply(const CrsMatrix< MatrixValue, Device, Layout > &A, const InputMultiVectorType &x, OutputMultiVectorType &y, const std::vector< OrdinalType > &col_indices, SingleColumnMultivectorMultiply)
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
static void apply(const MatrixType matrix)
MatrixType::size_type size_type
MatrixType::ordinal_type ordinal_type
const MatrixType m_matrix
MatrixType::value_type value_type
MatrixType::execution_space execution_space
AddDiagonalValuesAtomicKernel(const MatrixType matrix)
MatrixType::value_type value_type
static void apply(const MatrixType matrix)
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
const MatrixType m_matrix
MatrixType::size_type size_type
AddDiagonalValuesKernel(const MatrixType matrix)
MatrixType::execution_space execution_space
MatrixType::ordinal_type ordinal_type
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
void operator()(const Matrix &A, const InputVector &x, OutputVector &y) const
static ViewType create_view(const std::string &name, const OrdinalType &fem_length, const OrdinalType &stoch_length)
static ViewType create_view(const std::string &name, const OrdinalType &fem_length, const OrdinalType &stoch_length)
static ViewType create_view(const std::string &name, const OrdinalType &fem_length, const OrdinalType &stoch_length)
MatrixType::execution_space execution_space
static void apply(const MatrixType matrix)
MatrixType::size_type size_type
MatrixType::ordinal_type ordinal_type
ReplaceDiagonalValuesKernel(const MatrixType matrix)
KOKKOS_INLINE_FUNCTION void operator()(const size_type i) const
const MatrixType m_matrix
static bool check(const MatrixType matrix, Teuchos::FancyOStream &out)
MatrixType::value_type value_type
void operator()(const Matrix &A, const InputVector &x, OutputVector &y) const
Stokhos_MV_Multiply_Op(const Tag &tg=Tag())