42#include "Teuchos_UnitTestHelpers.hpp"
46#include "KokkosSparse_CrsMatrix.hpp"
47#include "KokkosSparse_spmv.hpp"
55#include "Kokkos_Core.hpp"
58template<
typename IntType >
65 return k + N * (
j + N * i );
68template <
typename ordinal >
71 std::vector< std::vector<ordinal> >& graph )
73 graph.resize( N * N * N, std::vector<ordinal>() );
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 ) {
83 graph[row].reserve(27);
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 ) {
93 graph[row].push_back(col);
96 total += graph[row].size();
102template <
typename scalar,
typename ordinal>
105 const ordinal nStoch,
106 const ordinal iRowFEM,
107 const ordinal iColFEM,
108 const ordinal iStoch )
110 const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
111 ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
113 const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
115 return A_fem + A_stoch;
119template <
typename scalar,
typename ordinal>
122 const ordinal nStoch,
123 const ordinal iColFEM,
124 const ordinal iStoch )
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;
132template <
typename kokkos_cijk_type,
typename ordinal_type>
138 using Teuchos::Array;
140 typedef typename kokkos_cijk_type::value_type value_type;
148 Array< RCP<const one_d_basis> > bases(
stoch_dim);
150 bases[i] = rcp(
new legendre_basis(
poly_ord,
true));
151 RCP<const product_basis> basis = rcp(
new product_basis(bases));
154 RCP<Cijk> cijk = basis->computeTripleProductTensor();
157 kokkos_cijk_type kokkos_cijk =
158 Stokhos::create_product_tensor<execution_space>(*basis, *cijk);
164template <
typename Scalar>
struct ScalarTol {};
165template <>
struct ScalarTol<float> {
static float tol() {
return 1e-4; } };
169template <
typename array_type,
typename scalar_type>
171 const array_type& y_exp,
172 const scalar_type rel_tol,
173 const scalar_type abs_tol,
174 Teuchos::FancyOStream& out)
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);
182 size_type num_rows = y.extent(0);
183 size_type num_cols = y.extent(1);
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;
190 out <<
"y_expected(" << i <<
"," <<
j <<
") - "
191 <<
"y(" << i <<
"," <<
j <<
") = " << hy_exp(i,
j)
192 <<
" - " << hy(i,
j) <<
" == "
193 << diff <<
" < " << tol <<
" : ";
199 success = success && s;
206template <
typename vector_type,
typename scalar_type>
208 const vector_type& y_exp,
209 const scalar_type rel_tol,
210 const scalar_type abs_tol,
211 Teuchos::FancyOStream& out)
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);
219 size_type num_rows = y.extent(0);
221 for (size_type i=0; i<num_rows; ++i) {
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 <<
" : ";
235 success = success && s;
241template <
typename vector_type,
typename scalar_type>
243 const vector_type& y_exp,
244 const scalar_type rel_tol,
245 const scalar_type abs_tol,
246 Teuchos::FancyOStream& out)
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);
254 size_type num_rows = y.extent(0);
255 size_type num_cols = y.extent(1);
258 for (size_type col = 0; col < num_cols; ++col){
259 for (size_type i=0; i<num_rows; ++i) {
262 scalar_type tol = rel_tol*std::abs(hy_exp(i,col).
fastAccessCoeff(
j)) + abs_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 <<
" : ";
273 success = success && s;
284template <
typename MatrixType,
typename CijkType>
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;
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;
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);
303 MatrixType matrix(
"matrix", nrow, matrix_values, matrix_graph);
312template <
typename MatrixType>
323 KOKKOS_INLINE_FUNCTION
328 m_matrix.replaceValues(row, &col, 1, &
val,
false,
true);
332 static void apply(
const MatrixType matrix) {
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);
349 bool s = compareVecs(host_matrix_values(row),
350 "matrix_values(row)",
354 success = success && s;
361template <
typename MatrixType>
372 KOKKOS_INLINE_FUNCTION
377 m_matrix.sumIntoValues(row, &col, 1, &
val,
false,
true);
381 static void apply(
const MatrixType matrix) {
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);
398 bool s = compareVecs(host_matrix_values(row),
399 "matrix_values(row)",
403 success = success && s;
411template <
typename MatrixType>
422 KOKKOS_INLINE_FUNCTION
427 m_matrix.sumIntoValues(row, &col, 1, &
val,
false,
true);
431 static void apply(
const MatrixType matrix) {
437 static bool check(
const MatrixType matrix,
438 Teuchos::FancyOStream& out) {
439 typedef typename MatrixType::values_type matrix_values_type;
440 typename matrix_values_type::HostMirror host_matrix_values =
441 Kokkos::create_mirror_view(matrix.values);
442 Kokkos::deep_copy(host_matrix_values, matrix.values);
449 val_expected = nrow*(nrow-1)/2;
452 bool s = compareVecs(host_matrix_values(row),
453 "matrix_values(row)",
457 success = success && s;
464 Kokkos_CrsMatrix_PCE, ReplaceValues, MatrixScalar )
466 typedef typename MatrixScalar::ordinal_type
Ordinal;
467 typedef typename MatrixScalar::execution_space Device;
468 typedef typename MatrixScalar::cijk_type Cijk;
469 typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
478 const Ordinal pce_size = cijk.dimension();
479 Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
483 kernel::apply(matrix);
486 success = kernel::check(matrix, out);
490 Kokkos_CrsMatrix_PCE, SumIntoValues, MatrixScalar )
492 typedef typename MatrixScalar::ordinal_type
Ordinal;
493 typedef typename MatrixScalar::execution_space Device;
494 typedef typename MatrixScalar::cijk_type Cijk;
495 typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
504 const Ordinal pce_size = cijk.dimension();
505 Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
509 kernel::apply(matrix);
512 success = kernel::check(matrix, out);
516 Kokkos_CrsMatrix_PCE, SumIntoValuesAtomic, MatrixScalar )
518 typedef typename MatrixScalar::ordinal_type
Ordinal;
519 typedef typename MatrixScalar::execution_space Device;
520 typedef typename MatrixScalar::cijk_type Cijk;
521 typedef KokkosSparse::CrsMatrix<MatrixScalar,Ordinal,Device> Matrix;
530 const Ordinal pce_size = cijk.dimension();
531 Matrix matrix = buildDiagonalMatrix<Matrix>(nrow, pce_size, cijk);
535 kernel::apply(matrix);
538 success = kernel::check(matrix, out);
541template <
typename PCEType,
typename Multiply>
543 const typename PCEType::ordinal_type
stoch_dim,
544 const typename PCEType::ordinal_type
poly_ord,
545 KokkosSparse::DeviceConfig dev_config,
546 Multiply multiply_op,
547 Teuchos::FancyOStream& out)
549 typedef typename PCEType::ordinal_type ordinal_type;
550 typedef typename PCEType::value_type scalar_type;
552 typedef typename PCEType::cijk_type cijk_type;
554 typedef Kokkos::LayoutLeft Layout;
555 typedef Kokkos::View< PCEType*, Layout, execution_space > block_vector_type;
556 typedef KokkosSparse::CrsMatrix< PCEType, ordinal_type, execution_space > block_matrix_type;
557 typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
558 typedef typename block_matrix_type::values_type matrix_values_type;
562 const ordinal_type stoch_length = cijk.dimension();
565 const ordinal_type stoch_length_aligned = stoch_length;
568 TEUCHOS_TEST_FOR_EXCEPTION(
569 storage_type::is_static && storage_type::static_size != stoch_length,
571 "Static storage size must equal pce size");
574 const ordinal_type fem_length = nGrid * nGrid * nGrid;
575 std::vector< std::vector<ordinal_type> > fem_graph;
582 block_vector_type x =
583 Kokkos::make_view<block_vector_type>(
"x", cijk, fem_length, stoch_length_aligned);
584 block_vector_type y =
585 Kokkos::make_view<block_vector_type>(
"y", cijk, fem_length, stoch_length_aligned);
587 typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
588 typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
591 typename block_vector_type::HostMirror::array_type hax = hx ;
592 typename block_vector_type::HostMirror::array_type hay = hy ;
594 for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
595 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
596 hax(iRowStoch, iRowFEM) =
597 generate_vector_coefficient<scalar_type>(
598 fem_length, stoch_length, iRowFEM, iRowStoch );
599 hay(iRowStoch, iRowFEM) = 0.0;
603 Kokkos::deep_copy( x, hx );
604 Kokkos::deep_copy( y, hy );
609 matrix_graph_type matrix_graph =
610 Kokkos::create_staticcrsgraph<matrix_graph_type>(
611 std::string(
"test crs graph"), fem_graph);
612 matrix_values_type matrix_values =
613 Kokkos::make_view<matrix_values_type>(
614 Kokkos::ViewAllocateWithoutInitializing(
"matrix"), cijk, fem_graph_length, stoch_length_aligned);
615 block_matrix_type matrix(
616 "block_matrix", fem_length, matrix_values, matrix_graph);
617 matrix.dev_config = dev_config;
619 typename matrix_values_type::HostMirror hM =
620 Kokkos::create_mirror_view( matrix.values );
622 typename matrix_values_type::HostMirror::array_type haM = hM ;
624 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
625 const ordinal_type row_size = fem_graph[iRowFEM].size();
626 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
627 ++iRowEntryFEM, ++iEntryFEM) {
628 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
630 for (ordinal_type k=0; k<stoch_length; ++k) {
632 generate_matrix_coefficient<scalar_type>(
633 fem_length, stoch_length, iRowFEM, iColFEM, k);
638 Kokkos::deep_copy( matrix.values, hM );
643 multiply_op( matrix, x, y );
648 typedef typename block_vector_type::array_type array_type;
649 array_type ay_expected =
650 array_type(
"ay_expected", stoch_length_aligned, fem_length);
651 typename array_type::HostMirror hay_expected =
652 Kokkos::create_mirror_view(ay_expected);
653 typename cijk_type::HostMirror host_cijk =
654 Kokkos::create_mirror_view(cijk);
655 Kokkos::deep_copy(host_cijk, cijk);
656 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
657 const ordinal_type row_size = fem_graph[iRowFEM].size();
658 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
659 ++iRowEntryFEM, ++iEntryFEM) {
660 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
661 for (ordinal_type i=0; i<stoch_length; ++i) {
662 const ordinal_type num_entry = host_cijk.num_entry(i);
663 const ordinal_type entry_beg = host_cijk.entry_begin(i);
664 const ordinal_type entry_end = entry_beg + num_entry;
666 for (ordinal_type entry = entry_beg; entry < entry_end; ++entry) {
667 const ordinal_type
j = host_cijk.coord(entry,0);
668 const ordinal_type k = host_cijk.coord(entry,1);
669 const scalar_type a_j =
670 generate_matrix_coefficient<scalar_type>(
671 fem_length, stoch_length, iRowFEM, iColFEM,
j);
672 const scalar_type a_k =
673 generate_matrix_coefficient<scalar_type>(
674 fem_length, stoch_length, iRowFEM, iColFEM, k);
675 const scalar_type x_j =
676 generate_vector_coefficient<scalar_type>(
677 fem_length, stoch_length, iColFEM,
j);
678 const scalar_type x_k =
679 generate_vector_coefficient<scalar_type>(
680 fem_length, stoch_length, iColFEM, k);
681 tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
683 hay_expected(i, iRowFEM) += tmp;
687 Kokkos::deep_copy( ay_expected, hay_expected );
692 typename block_vector_type::array_type ay = y;
701 Kokkos_CrsMatrix_PCE, Multiply,
Scalar, MultiplyOp )
703 typedef typename Scalar::ordinal_type
Ordinal;
708 KokkosSparse::DeviceConfig dev_config;
710 success = test_embedded_pce<Scalar>(
715 template <
typename Matrix,
typename InputVector,
typename OutputVector>
717 const InputVector& x,
718 OutputVector& y)
const {
719 KokkosSparse::spmv(
"N",
typename Matrix::value_type(1.0) , A, x,
typename Matrix::value_type(0.0), y);
723template <
typename Tag>
728 template <
typename Matrix,
typename InputVector,
typename OutputVector>
730 const InputVector& x,
731 OutputVector& y)
const {
737 Kokkos_CrsMatrix_PCE, MeanMultiplyRank1,
Scalar )
739 typedef typename Scalar::ordinal_type
Ordinal;
744 KokkosSparse::DeviceConfig dev_config;
746 typedef typename Scalar::ordinal_type ordinal_type;
747 typedef typename Scalar::value_type scalar_type;
749 typedef typename Scalar::cijk_type cijk_type;
751 typedef Kokkos::LayoutLeft Layout;
752 typedef Kokkos::View< Scalar*, Layout, execution_space > block_vector_type;
753 typedef KokkosSparse::CrsMatrix< Scalar, ordinal_type, execution_space > block_matrix_type;
754 typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
755 typedef typename block_matrix_type::values_type matrix_values_type;
759 cijk_type mean_cijk =
760 Stokhos::create_mean_based_product_tensor<execution_space, ordinal_type, scalar_type>();
761 const ordinal_type stoch_length = cijk.dimension();
762 const ordinal_type align = 8;
763 const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
766 const ordinal_type fem_length = nGrid * nGrid * nGrid;
767 std::vector< std::vector<ordinal_type> > fem_graph;
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);
775 block_vector_type y_expected =
776 Kokkos::make_view<block_vector_type>(
"y", cijk, fem_length, stoch_length_aligned);
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 );
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 =
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;
798 Kokkos::deep_copy( x, hx );
799 Kokkos::deep_copy( y, hy );
800 Kokkos::deep_copy( y_expected, hy_expected );
804 matrix_graph_type matrix_graph =
805 Kokkos::create_staticcrsgraph<matrix_graph_type>(
806 std::string(
"test crs graph"), fem_graph);
807 matrix_values_type matrix_values =
808 Kokkos::make_view<matrix_values_type>(
809 Kokkos::ViewAllocateWithoutInitializing(
"matrix"), mean_cijk, fem_graph_length, ordinal_type(1));
810 block_matrix_type matrix(
811 "block_matrix", fem_length, matrix_values, matrix_graph);
812 matrix.dev_config = dev_config;
814 typename matrix_values_type::HostMirror hM =
815 Kokkos::create_mirror_view( matrix.values );
816 typename matrix_values_type::HostMirror::array_type haM = hM ;
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];
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);
833 Kokkos::deep_copy( matrix.values, hM );
834 Kokkos::deep_copy( y_expected, hy_expected );
883 success =
compareRank1(y, y_expected, rel_tol, abs_tol, out);
888 Kokkos_CrsMatrix_PCE, MeanMultiplyRank2,
Scalar )
890 typedef typename Scalar::ordinal_type ordinal_type;
891 typedef typename Scalar::value_type scalar_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;
902 const ordinal_type nGrid = 5;
905 KokkosSparse::DeviceConfig dev_config;
909 cijk_type mean_cijk =
910 Stokhos::create_mean_based_product_tensor<execution_space, ordinal_type, scalar_type>();
911 const ordinal_type stoch_length = cijk.dimension();
912 const ordinal_type align = 8;
913 const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
914 const ordinal_type num_cols = 2;
916 const ordinal_type fem_length = nGrid * nGrid * nGrid;
917 std::vector< std::vector<ordinal_type> > fem_graph;
920 block_vector_type x =
921 Kokkos::make_view<block_vector_type>(
"x", cijk, fem_length, num_cols, stoch_length_aligned);
922 block_vector_type y =
923 Kokkos::make_view<block_vector_type>(
"y", cijk, fem_length, num_cols, stoch_length_aligned);
925 block_vector_type y_expected =
926 Kokkos::make_view<block_vector_type>(
"y_expected", cijk, fem_length, num_cols, stoch_length_aligned);
928 typename block_vector_type::HostMirror hx = Kokkos::create_mirror_view( x );
929 typename block_vector_type::HostMirror hy = Kokkos::create_mirror_view( y );
930 typename block_vector_type::HostMirror hy_expected =
931 Kokkos::create_mirror_view( y_expected );
933 for (ordinal_type i=0; i<num_cols; ++i){
934 for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
935 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
936 hx(iRowFEM,i).fastAccessCoeff(iRowStoch) =
937 generate_vector_coefficient<scalar_type>(
938 fem_length, stoch_length, iRowFEM, iRowStoch );
939 hy(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
940 hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
944 Kokkos::deep_copy( x, hx );
945 Kokkos::deep_copy( y, hy );
946 Kokkos::deep_copy( y_expected, hy_expected );
950 matrix_graph_type matrix_graph =
951 Kokkos::create_staticcrsgraph<matrix_graph_type>(
952 std::string(
"test crs graph"), fem_graph);
953 matrix_values_type matrix_values =
954 Kokkos::make_view<matrix_values_type>(
955 Kokkos::ViewAllocateWithoutInitializing(
"matrix"), mean_cijk, fem_graph_length, ordinal_type(1));
956 block_matrix_type matrix(
957 "block_matrix", fem_length, matrix_values, matrix_graph);
958 matrix.dev_config = dev_config;
960 typename matrix_values_type::HostMirror hM =
961 Kokkos::create_mirror_view( matrix.values );
963 typename matrix_values_type::HostMirror::array_type haM = hM ;
965 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
966 const ordinal_type row_size = fem_graph[iRowFEM].size();
967 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
968 ++iRowEntryFEM, ++iEntryFEM) {
969 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
972 generate_matrix_coefficient<scalar_type>(
973 fem_length, 1, iRowFEM, iColFEM, 0);
974 for (ordinal_type i=0; i<num_cols; ++i){
975 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
976 hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) +=
977 haM(iEntryFEM, 0) * hx(iColFEM,i).fastAccessCoeff(iRowStoch);
983 Kokkos::deep_copy( matrix.values, hM );
984 Kokkos::deep_copy( y_expected, hy_expected );
1035 success =
compareRank2(y, y_expected, rel_tol, abs_tol, out);
1041#define CRSMATRIX_UQ_PCE_TESTS_MATRIXSCALAR( SCALAR ) \
1042 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1043 Kokkos_CrsMatrix_PCE, ReplaceValues, SCALAR ) \
1044 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1045 Kokkos_CrsMatrix_PCE, SumIntoValues, SCALAR ) \
1046 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1047 Kokkos_CrsMatrix_PCE, SumIntoValuesAtomic, SCALAR )
1048#define CRSMATRIX_UQ_PCE_MEAN_MULTIPLY_TESTS( SCALAR ) \
1049 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1050 Kokkos_CrsMatrix_PCE, MeanMultiplyRank1, SCALAR ) \
1051 TEUCHOS_UNIT_TEST_TEMPLATE_1_INSTANT( \
1052 Kokkos_CrsMatrix_PCE, MeanMultiplyRank2, SCALAR )
1053#define CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR_OP( SCALAR, OP ) \
1054 TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT( \
1055 Kokkos_CrsMatrix_PCE, Multiply, SCALAR, OP )
1057#define CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR( SCALAR ) \
1058 CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR_OP( SCALAR, KokkosMultiply )
1060#define CRSMATRIX_UQ_PCE_TESTS_STORAGE( STORAGE ) \
1061 typedef Sacado::UQ::PCE<STORAGE> UQ_PCE_ ## STORAGE; \
1062 CRSMATRIX_UQ_PCE_TESTS_MATRIXSCALAR( UQ_PCE_ ## STORAGE ) \
1063 CRS_MATRIX_UQ_PCE_MULTIPLY_TESTS_SCALAR( UQ_PCE_ ## STORAGE ) \
1064 CRSMATRIX_UQ_PCE_MEAN_MULTIPLY_TESTS( UQ_PCE_ ## STORAGE )
1066#define CRSMATRIX_UQ_PCE_TESTS_ORDINAL_SCALAR_DEVICE( ORDINAL, SCALAR, DEVICE ) \
1067 typedef Stokhos::DynamicStorage<ORDINAL,SCALAR,DEVICE> DS; \
1068 CRSMATRIX_UQ_PCE_TESTS_STORAGE( DS )
1070#define CRSMATRIX_UQ_PCE_TESTS_DEVICE( DEVICE ) \
1071 CRSMATRIX_UQ_PCE_TESTS_ORDINAL_SCALAR_DEVICE( int, double, DEVICE )
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j) - expr2.val(j)
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)
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
kokkos_cijk_type build_cijk(ordinal_type stoch_dim, ordinal_type poly_ord)
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)
TEUCHOS_UNIT_TEST_TEMPLATE_1_DECL(Kokkos_CrsMatrix_PCE, ReplaceValues, MatrixScalar)
Kokkos_MV_Multiply_Op KokkosMultiply
MatrixType buildDiagonalMatrix(typename MatrixType::ordinal_type nrow, typename MatrixType::ordinal_type pce_size, const CijkType &cijk)
TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(Kokkos_CrsMatrix_PCE, Multiply, Scalar, MultiplyOp)
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)
bool compareRank2(const vector_type &y, const vector_type &y_exp, const scalar_type rel_tol, const scalar_type abs_tol, Teuchos::FancyOStream &out)
bool compareRank1(const vector_type &y, const vector_type &y_exp, const scalar_type rel_tol, const scalar_type abs_tol, Teuchos::FancyOStream &out)
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
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Abstract base class for 1-D orthogonal polynomials.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
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)
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typenameCijkType< view_type >::type >::type cijk(const view_type &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
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())