42 #ifndef STOKHOS_KOKKOS_CORE_KERNELS_UNIT_TEST_HPP
43 #define STOKHOS_KOKKOS_CORE_KERNELS_UNIT_TEST_HPP
51 #include "EpetraExt_BlockUtility.h"
60 #include "Kokkos_Macros.hpp"
76 #ifdef HAVE_STOKHOS_KOKKOSLINALG
77 #include "KokkosSparse_CrsMatrix.hpp"
78 #include "KokkosSparse_spmv.hpp"
79 #include "KokkosBlas1_update.hpp"
89 template<
typename IntType >
96 return k + N * (
j + N * i );
99 template <
typename ordinal >
102 std::vector< std::vector<ordinal> > & graph )
104 graph.resize( N * N * N , std::vector<ordinal>() );
108 for (
int i = 0 ; i < (
int) N ; ++i ) {
109 for (
int j = 0 ;
j < (
int) N ; ++
j ) {
110 for (
int k = 0 ; k < (
int) N ; ++k ) {
114 graph[row].reserve(27);
116 for (
int ii = -1 ; ii < 2 ; ++ii ) {
117 for (
int jj = -1 ; jj < 2 ; ++jj ) {
118 for (
int kk = -1 ; kk < 2 ; ++kk ) {
119 if ( 0 <= i + ii && i + ii < (
int) N &&
120 0 <=
j + jj &&
j + jj < (
int) N &&
121 0 <= k + kk && k + kk < (
int) N ) {
124 graph[row].push_back(col);
127 total += graph[row].size();
133 template <
typename scalar,
typename ordinal>
136 const ordinal nStoch ,
137 const ordinal iRowFEM ,
138 const ordinal iColFEM ,
139 const ordinal iStoch )
141 const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
142 ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
144 const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
146 return A_fem + A_stoch ;
150 template <
typename scalar,
typename ordinal>
153 const ordinal nStoch ,
154 const ordinal iColFEM ,
155 const ordinal iStoch )
157 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
158 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
159 return X_fem + X_stoch ;
163 template <
typename Device>
184 void setup(
int p_ = 5,
int d_ = 2) {
208 for (
int i=0; i<
d; i++)
212 Cijk =
basis->computeTripleProductTensor();
233 rcp(EpetraExt::BlockUtility::GenerateBlockMap(
234 *x_map, *stoch_row_map, *sg_comm));
239 for (
int i=0; i<num_my_GIDs; ++i) {
240 int row = my_GIDs[i];
250 x_map, x_map, sg_x_map, sg_x_map,
257 basis, sg_A_overlap_map, x_map, x_map, sg_x_map, sg_comm));
261 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
fem_length ; ++iRowFEM ) {
262 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
263 const int iColFEM =
fem_graph[iRowFEM][iRowEntryFEM] ;
265 double v = generate_matrix_coefficient<double>(
267 A->ReplaceGlobalValues(iRowFEM, 1, &v, &iColFEM);
271 A_sg_blocks->setCoeffPtr(i,
A);
276 basis, stoch_row_map, x_map, sg_x_map, sg_comm));
278 basis, stoch_row_map, x_map, sg_x_map, sg_comm));
281 for (
int iColFEM=0; iColFEM <
fem_length; ++iColFEM ) {
282 for (
int iColStoch=0 ; iColStoch <
stoch_length; ++iColStoch ) {
283 (*sg_x)[iColStoch][iColFEM] = generate_vector_coefficient<double>(
294 x_map, stoch_row_map, sg_comm));
300 typedef typename Kokkos::ViewTraits<value_type,Device,void,void>::host_mirror_space host_device;
303 host_device > stoch_tensor_type ;
305 stoch_tensor_type tensor =
306 Stokhos::create_stochastic_product_tensor< tensor_type >( *
basis, *
Cijk );
313 for (
int j=0;
j<
d; ++
j)
314 term[
j] = tensor.bases_degree(i,
j);
315 int idx =
basis->index(term);
321 template <
typename vec_type>
331 out <<
"y_expected[" << block <<
"][" << i <<
"] - "
332 <<
"y[" << block <<
"][" << i <<
"] = " << (*sg_y)[block][i]
333 <<
" - " <<
y[block][i] <<
" == "
334 << diff <<
" < " << tol <<
" : failed"
337 success = success && s;
344 template <
typename multi_vec_type>
354 out <<
"y_expected[" << block <<
"][" << i <<
"] - "
355 <<
"y(" << i <<
"," << block <<
") = " << (*sg_y)[block][i]
356 <<
" - " <<
y(i,block) <<
" == "
357 << diff <<
" < " << tol <<
" : failed"
360 success = success && s;
367 template <
typename vec_type>
378 out <<
"y_expected[" << block <<
"][" << i <<
"] - "
379 <<
"y(" << b <<
"," << i <<
") = " << (*sg_y)[block][i]
380 <<
" - " <<
y(b,i) <<
" == "
381 << diff <<
" < " << tol <<
" : failed"
384 success = success && s;
391 template <
typename vec_type>
402 out <<
"y_expected[" << block <<
"][" << i <<
"] - "
403 <<
"y(" << b <<
"," << i <<
") = " << (*sg_y)[block][i] <<
" - "
405 << diff <<
" < " << tol <<
" : failed"
408 success = success && s;
415 template <
typename vec_type>
426 out <<
"y_expected[" << block <<
"][" << i <<
"] - "
427 <<
"y(" << b <<
"," << i <<
") == "
428 << diff <<
" < " << tol <<
" : failed"
431 success = success && s;
438 template <
typename vec_type>
449 out <<
"y_expected[" << block <<
"][" << i <<
"] - "
450 <<
"y(" << b <<
"," << i <<
") == "
451 << diff <<
" < " << tol <<
" : failed"
454 success = success && s;
463 template <
typename value_type,
typename Device,
typename SparseMatOps>
467 typedef typename matrix_type::values_type matrix_values_type;
468 typedef typename matrix_type::graph_type matrix_graph_type;
469 typedef Kokkos::View<value_type*,Device>
vec_type;
473 std::vector<matrix_type> matrix(
setup.stoch_length ) ;
474 std::vector<vec_type>
x(
setup.stoch_length ) ;
475 std::vector<vec_type>
y(
setup.stoch_length ) ;
476 std::vector<vec_type> tmp(
setup.stoch_length ) ;
478 for (
int block=0; block<
setup.stoch_length; ++block) {
479 matrix[block].graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
480 std::string(
"testing") ,
setup.fem_graph );
482 matrix[block].values =
483 matrix_values_type(
"matrix" ,
setup.fem_graph_length );
489 typename matrix_values_type::HostMirror hM =
492 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
493 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
494 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
496 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
497 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , block );
503 typename vec_type::HostMirror hx =
505 typename vec_type::HostMirror hy =
508 for (
int i = 0 ; i <
setup.fem_length ; ++i ) {
509 hx(i) = generate_vector_coefficient<value_type>(
510 setup.fem_length ,
setup.stoch_length , i , block );
521 setup.Cijk->k_begin();
525 k_it!=k_end; ++k_it) {
526 int nj =
setup.Cijk->num_j(k_it);
530 setup.Cijk->j_begin(k_it);
532 setup.Cijk->j_end(k_it);
533 std::vector<vec_type> xx(nj), yy(nj);
546 j_it != j_end; ++j_it) {
548 setup.Cijk->i_begin(j_it);
550 setup.Cijk->i_end(j_it);
563 std::vector<typename vec_type::HostMirror> hy(
setup.stoch_length);
564 for (
int i=0; i<
setup.stoch_length; ++i) {
568 bool success =
setup.test_original(hy, out);
573 template <
typename value_type,
typename Device,
typename SparseMatOps>
577 typedef typename matrix_type::values_type matrix_values_type;
578 typedef typename matrix_type::graph_type matrix_graph_type;
579 typedef Kokkos::View<value_type*, Kokkos::LayoutLeft, Device, Kokkos::MemoryUnmanaged>
vec_type;
580 typedef Kokkos::View<value_type**, Kokkos::LayoutLeft, Device> multi_vec_type;
584 std::vector<matrix_type> matrix(
setup.stoch_length ) ;
585 multi_vec_type
x(
"x",
setup.fem_length,
setup.stoch_length ) ;
586 multi_vec_type
y(
"y",
setup.fem_length,
setup.stoch_length ) ;
587 multi_vec_type tmp_x(
"tmp_x",
setup.fem_length,
setup.stoch_length ) ;
588 multi_vec_type tmp_y(
"tmp_y",
setup.fem_length,
setup.stoch_length ) ;
593 for (
int block=0; block<
setup.stoch_length; ++block) {
594 matrix[block].graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
595 std::string(
"testing") ,
setup.fem_graph );
597 matrix[block].values =
598 matrix_values_type(
"matrix" ,
setup.fem_graph_length );
600 typename matrix_values_type::HostMirror hM =
603 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
604 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
605 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
607 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
608 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , block );
614 for (
int i = 0 ; i <
setup.fem_length ; ++i ) {
615 hx(i, block) = generate_vector_coefficient<value_type>(
616 setup.fem_length ,
setup.stoch_length , i , block );
630 k_iterator k_begin =
setup.Cijk->k_begin();
631 k_iterator k_end =
setup.Cijk->k_end();
632 for (k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
633 unsigned nj =
setup.Cijk->num_j(k_it);
636 kj_iterator j_begin =
setup.Cijk->j_begin(k_it);
637 kj_iterator j_end =
setup.Cijk->j_end(k_it);
638 std::vector<int> j_indices(nj);
640 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
642 vec_type xx = Kokkos::subview(
x, Kokkos::ALL(),
j );
643 vec_type tt = Kokkos::subview( tmp_x, Kokkos::ALL(), jdx++ );
646 multi_vec_type tmp_x_view =
647 Kokkos::subview( tmp_x, Kokkos::ALL(),
648 std::make_pair(0u,nj));
649 multi_vec_type tmp_y_view =
650 Kokkos::subview( tmp_y, Kokkos::ALL(),
651 std::make_pair(0u,nj));
654 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
656 Kokkos::subview( tmp_y, Kokkos::ALL(), jdx++ );
657 kji_iterator i_begin =
setup.Cijk->i_begin(j_it);
658 kji_iterator i_end =
setup.Cijk->i_end(j_it);
659 for (kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
662 vec_type y_view = Kokkos::subview(
y, Kokkos::ALL(), i );
670 bool success =
setup.test_original(hy, out);
675 #ifdef HAVE_STOKHOS_KOKKOSLINALG
677 template <
typename value_type,
typename Device>
681 typedef KokkosSparse::CrsMatrix<value_type,ordinal_type,Device> matrix_type;
682 typedef typename matrix_type::values_type matrix_values_type;
683 typedef typename matrix_type::StaticCrsGraphType matrix_graph_type;
684 typedef Kokkos::View<value_type*, Kokkos::LayoutLeft, Device, Kokkos::MemoryUnmanaged>
vec_type;
685 typedef Kokkos::View<value_type**, Kokkos::LayoutLeft, Device> multi_vec_type;
689 std::vector<matrix_type> matrix(
setup.stoch_length ) ;
690 multi_vec_type
x(
"x",
setup.fem_length,
setup.stoch_length ) ;
691 multi_vec_type
y(
"y",
setup.fem_length,
setup.stoch_length ) ;
692 multi_vec_type tmp_x(
"tmp_x",
setup.fem_length,
setup.stoch_length ) ;
693 multi_vec_type tmp_y(
"tmp_y",
setup.fem_length,
setup.stoch_length ) ;
698 for (
int block=0; block<
setup.stoch_length; ++block) {
699 matrix_graph_type matrix_graph =
700 Kokkos::create_staticcrsgraph<matrix_graph_type>(
701 std::string(
"test crs graph"),
setup.fem_graph);
702 matrix_values_type matrix_values =
703 matrix_values_type(
"matrix" ,
setup.fem_graph_length );
704 typename matrix_values_type::HostMirror hM =
707 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
708 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
709 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
711 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
712 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , block );
717 matrix[block] = matrix_type(
"matrix",
setup.fem_length, matrix_values,
720 for (
int i = 0 ; i <
setup.fem_length ; ++i ) {
721 hx(i, block) = generate_vector_coefficient<value_type>(
722 setup.fem_length ,
setup.stoch_length , i , block );
735 k_iterator k_begin =
setup.Cijk->k_begin();
736 k_iterator k_end =
setup.Cijk->k_end();
737 for (k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
738 int nj =
setup.Cijk->num_j(k_it);
741 kj_iterator j_begin =
setup.Cijk->j_begin(k_it);
742 kj_iterator j_end =
setup.Cijk->j_end(k_it);
744 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
746 vec_type xx = Kokkos::subview(
x, Kokkos::ALL(),
j );
747 vec_type tt = Kokkos::subview( tmp_x, Kokkos::ALL(), jdx++ );
750 multi_vec_type tmp_x_view =
751 Kokkos::subview( tmp_x, Kokkos::ALL(),
752 std::make_pair(0u,jdx));
753 multi_vec_type tmp_y_view =
754 Kokkos::subview( tmp_y, Kokkos::ALL(),
755 std::make_pair(0u,jdx));
758 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
760 Kokkos::subview( tmp_y, Kokkos::ALL(), jdx++ );
761 kji_iterator i_begin =
setup.Cijk->i_begin(j_it);
762 kji_iterator i_end =
setup.Cijk->i_end(j_it);
763 for (kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
766 vec_type y_view = Kokkos::subview(
y, Kokkos::ALL(), i );
775 bool success =
setup.test_original(hy, out);
782 template<
typename ScalarType ,
class Device >
788 typedef Kokkos::View< value_type**, Kokkos::LayoutLeft , Device > block_vector_type ;
794 typedef typename matrix_type::graph_type graph_type ;
801 for (k_iterator k_it=
setup.Cijk->k_begin();
802 k_it!=
setup.Cijk->k_end(); ++k_it) {
804 for (kj_iterator j_it =
setup.Cijk->j_begin(k_it);
805 j_it !=
setup.Cijk->j_end(k_it); ++j_it) {
807 for (kji_iterator i_it =
setup.Cijk->i_begin(j_it);
808 i_it !=
setup.Cijk->i_end(j_it); ++i_it) {
810 double c = value(i_it);
811 dense_cijk(i,
j,k) = c;
819 block_vector_type
x =
820 block_vector_type(
"x" ,
setup.stoch_length ,
setup.fem_length );
821 block_vector_type
y =
822 block_vector_type(
"y" ,
setup.stoch_length ,
setup.fem_length );
826 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
827 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
828 hx(iColStoch,iColFEM) =
829 generate_vector_coefficient<ScalarType>(
830 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
843 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
844 std::string(
"test crs graph") ,
setup.fem_graph );
845 matrix.values = block_vector_type(
846 "matrix" , matrix.block.matrix_size() ,
setup.fem_graph_length );
849 typename block_vector_type::HostMirror hM =
852 for (
int iRowStoch = 0 ; iRowStoch <
setup.stoch_length ; ++iRowStoch ) {
853 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
855 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
856 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM];
858 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
860 const size_t offset =
861 matrix.block.matrix_offset( iRowStoch , iColStoch );
863 ScalarType value = 0 ;
865 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
866 value += dense_cijk( iRowStoch , iColStoch , k ) *
867 generate_matrix_coefficient<ScalarType>(
868 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , k );
871 hM( offset , iEntryFEM ) = value ;
888 bool success =
setup.test_commuted(hy, out);
893 template<
typename ScalarType ,
class Device >
899 typedef Kokkos::View< value_type* , Device > vector_type ;
904 typedef typename matrix_type::values_type matrix_values_type;
905 typedef typename matrix_type::graph_type matrix_graph_type;
908 std::vector< std::vector< int > > stoch_graph(
setup.stoch_length );
911 for (
int i = 0 ; i <
setup.stoch_length ; ++i ) {
913 stoch_graph[i].resize(len);
923 for (k_iterator k_it=
setup.Cijk->k_begin();
924 k_it!=
setup.Cijk->k_end(); ++k_it) {
926 for (kj_iterator j_it =
setup.Cijk->j_begin(k_it);
927 j_it !=
setup.Cijk->j_end(k_it); ++j_it) {
929 for (kji_iterator i_it =
setup.Cijk->i_begin(j_it);
930 i_it !=
setup.Cijk->i_end(j_it); ++i_it) {
932 double c = value(i_it);
933 dense_cijk(i,
j,k) = c;
941 const int flat_length =
setup.fem_length *
setup.stoch_length ;
943 std::vector< std::vector<int> > flat_graph( flat_length );
945 for (
int iOuterRow = 0 ; iOuterRow <
setup.fem_length ; ++iOuterRow ) {
947 const size_t iOuterRowNZ =
setup.fem_graph[iOuterRow].size();
949 for (
int iInnerRow = 0 ; iInnerRow <
setup.stoch_length ; ++iInnerRow ) {
951 const size_t iInnerRowNZ = stoch_graph[ iInnerRow ].size(); ;
952 const int iFlatRowNZ = iOuterRowNZ * iInnerRowNZ ;
953 const int iFlatRow = iInnerRow + iOuterRow *
setup.stoch_length ;
955 flat_graph[iFlatRow].resize( iFlatRowNZ );
959 for (
size_t iOuterEntry = 0 ; iOuterEntry < iOuterRowNZ ; ++iOuterEntry ) {
961 const int iOuterCol =
setup.fem_graph[iOuterRow][iOuterEntry];
963 for (
size_t iInnerEntry = 0 ; iInnerEntry < iInnerRowNZ ; ++iInnerEntry ) {
965 const int iInnerCol = stoch_graph[iInnerRow][iInnerEntry] ;
966 const int iFlatColumn = iInnerCol + iOuterCol *
setup.stoch_length ;
968 flat_graph[iFlatRow][iFlatEntry] = iFlatColumn ;
978 vector_type
x = vector_type(
"x" , flat_length );
979 vector_type
y = vector_type(
"y" , flat_length );
983 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
984 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
985 hx(iColStoch + iColFEM*
setup.stoch_length) =
986 generate_vector_coefficient<ScalarType>(
987 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
997 matrix.graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
998 std::string(
"testing") , flat_graph );
1000 const size_t flat_graph_length = matrix.graph.entries.extent(0);
1002 matrix.values = matrix_values_type(
"matrix" , flat_graph_length );
1004 typename matrix_values_type::HostMirror hM =
1007 for (
int iRow = 0 , iEntry = 0 ; iRow < flat_length ; ++iRow ) {
1008 const int iRowFEM = iRow /
setup.stoch_length ;
1009 const int iRowStoch = iRow %
setup.stoch_length ;
1011 for (
size_t iRowEntry = 0 ; iRowEntry < flat_graph[ iRow ].size() ; ++iRowEntry , ++iEntry ) {
1012 const int iCol = flat_graph[ iRow ][ iRowEntry ];
1013 const int iColFEM = iCol /
setup.stoch_length ;
1014 const int iColStoch = iCol %
setup.stoch_length ;
1017 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1018 const double A_fem_k =
1019 generate_matrix_coefficient<ScalarType>(
1020 setup.fem_length ,
setup.stoch_length , iRowFEM, iColFEM, k );
1021 value += dense_cijk(iRowStoch,iColStoch,k) * A_fem_k ;
1023 hM( iEntry ) = value ;
1038 bool success =
setup.test_commuted_flat(hy, out);
1042 template<
typename ScalarType ,
class Device >
1048 typedef Kokkos::View< value_type* , Device > vector_type ;
1053 typedef typename matrix_type::values_type matrix_values_type;
1054 typedef typename matrix_type::graph_type matrix_graph_type;
1057 std::vector< std::vector< int > > stoch_graph(
setup.stoch_length );
1060 for (
int i = 0 ; i <
setup.stoch_length ; ++i ) {
1062 stoch_graph[i].resize(len);
1072 for (k_iterator k_it=
setup.Cijk->k_begin();
1073 k_it!=
setup.Cijk->k_end(); ++k_it) {
1074 int k = index(k_it);
1075 for (kj_iterator j_it =
setup.Cijk->j_begin(k_it);
1076 j_it !=
setup.Cijk->j_end(k_it); ++j_it) {
1077 int j = index(j_it);
1078 for (kji_iterator i_it =
setup.Cijk->i_begin(j_it);
1079 i_it !=
setup.Cijk->i_end(j_it); ++i_it) {
1080 int i = index(i_it);
1081 double c = value(i_it);
1082 dense_cijk(i,
j,k) = c;
1090 const size_t flat_length =
setup.fem_length *
setup.stoch_length ;
1092 std::vector< std::vector<int> > flat_graph( flat_length );
1094 for (
int iOuterRow = 0 ; iOuterRow <
setup.stoch_length ; ++iOuterRow ) {
1096 const size_t iOuterRowNZ = stoch_graph[iOuterRow].size();
1098 for (
int iInnerRow = 0 ; iInnerRow <
setup.fem_length ; ++iInnerRow ) {
1100 const size_t iInnerRowNZ =
setup.fem_graph[iInnerRow].size();
1101 const int iFlatRowNZ = iOuterRowNZ * iInnerRowNZ ;
1102 const int iFlatRow = iInnerRow + iOuterRow *
setup.fem_length ;
1104 flat_graph[iFlatRow].resize( iFlatRowNZ );
1106 int iFlatEntry = 0 ;
1108 for (
size_t iOuterEntry = 0 ; iOuterEntry < iOuterRowNZ ; ++iOuterEntry ) {
1110 const int iOuterCol = stoch_graph[ iOuterRow ][ iOuterEntry ];
1112 for (
size_t iInnerEntry = 0 ; iInnerEntry < iInnerRowNZ ; ++iInnerEntry ) {
1114 const int iInnerCol =
setup.fem_graph[ iInnerRow][iInnerEntry];
1115 const int iFlatColumn = iInnerCol + iOuterCol *
setup.fem_length ;
1117 flat_graph[iFlatRow][iFlatEntry] = iFlatColumn ;
1126 vector_type
x = vector_type(
"x" , flat_length );
1127 vector_type
y = vector_type(
"y" , flat_length );
1131 for (
size_t iCol = 0 ; iCol < flat_length ; ++iCol ) {
1132 const int iColStoch = iCol /
setup.fem_length ;
1133 const int iColFEM = iCol %
setup.fem_length ;
1135 hx(iCol) = generate_vector_coefficient<ScalarType>(
1136 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
1143 matrix_type matrix ;
1145 matrix.graph = Kokkos::create_staticcrsgraph<matrix_graph_type>( std::string(
"testing") , flat_graph );
1147 const size_t flat_graph_length = matrix.graph.entries.extent(0);
1149 matrix.values = matrix_values_type(
"matrix" , flat_graph_length );
1151 typename matrix_values_type::HostMirror hM =
1154 for (
size_t iRow = 0 , iEntry = 0 ; iRow < flat_length ; ++iRow ) {
1155 const int iRowStoch = iRow /
setup.fem_length ;
1156 const int iRowFEM = iRow %
setup.fem_length ;
1158 for (
size_t iRowEntry = 0 ; iRowEntry < flat_graph[ iRow ].size() ; ++iRowEntry , ++iEntry ) {
1159 const int iCol = flat_graph[ iRow ][ iRowEntry ];
1160 const int iColStoch = iCol /
setup.fem_length ;
1161 const int iColFEM = iCol %
setup.fem_length ;
1164 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1165 const double A_fem_k =
1166 generate_matrix_coefficient<ScalarType>(
1168 iRowFEM , iColFEM , k );
1169 value += dense_cijk(iRowStoch,iColStoch,k) * A_fem_k ;
1171 hM( iEntry ) = value ;
1185 bool success =
setup.test_original_flat(hy, out);
1189 template<
typename ScalarType ,
typename TensorType,
class Device >
1195 typedef Kokkos::View<
value_type** , Kokkos::LayoutLeft ,
1196 Device > block_vector_type ;
1201 typedef typename matrix_type::graph_type graph_type ;
1206 block_vector_type
x =
1207 block_vector_type(
"x" ,
setup.stoch_length_aligned ,
setup.fem_length );
1208 block_vector_type
y =
1209 block_vector_type(
"y" ,
setup.stoch_length_aligned ,
setup.fem_length );
1211 typename block_vector_type::HostMirror hx =
1214 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
1215 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
1216 hx(
setup.perm[iColStoch],iColFEM) =
1217 generate_vector_coefficient<ScalarType>(
1218 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
1226 matrix_type matrix ;
1229 Stokhos::create_stochastic_product_tensor< TensorType >( *
setup.basis,
1233 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1234 std::string(
"test crs graph") ,
setup.fem_graph );
1236 matrix.values = block_vector_type(
1237 "matrix" ,
setup.stoch_length_aligned ,
setup.fem_graph_length );
1239 typename block_vector_type::HostMirror hM =
1242 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
1243 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1244 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1246 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1247 hM(
setup.perm[k],iEntryFEM) =
1248 generate_matrix_coefficient<ScalarType>(
1249 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , k );
1264 bool success =
setup.test_commuted_perm(hy, out);
1269 template<
typename ScalarType ,
class Device ,
int BlockSize >
1272 const bool symmetric) {
1274 typedef Kokkos::View<
value_type** , Kokkos::LayoutLeft ,
1275 Device > block_vector_type ;
1280 typedef typename matrix_type::graph_type graph_type ;
1283 matrix_type matrix ;
1286 params.
set(
"Symmetric",symmetric);
1288 Stokhos::create_stochastic_product_tensor< TensorType >( *
setup.basis,
1291 int aligned_stoch_length = matrix.block.tensor().aligned_dimension();
1296 block_vector_type
x =
1297 block_vector_type(
"x" , aligned_stoch_length ,
setup.fem_length );
1298 block_vector_type
y =
1299 block_vector_type(
"y" , aligned_stoch_length ,
setup.fem_length );
1301 typename block_vector_type::HostMirror hx =
1304 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
1305 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
1306 hx(iColStoch,iColFEM) =
1307 generate_vector_coefficient<ScalarType>(
1308 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
1317 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1318 std::string(
"test crs graph") ,
setup.fem_graph );
1320 matrix.values = block_vector_type(
1321 "matrix" , aligned_stoch_length ,
setup.fem_graph_length );
1323 typename block_vector_type::HostMirror hM =
1326 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
1327 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1328 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1330 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1332 generate_matrix_coefficient<ScalarType>(
1333 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , k );
1348 bool success =
setup.test_commuted(hy, out);
1353 template<
typename ScalarType ,
class Device >
1358 typedef Kokkos::View<
value_type** , Kokkos::LayoutLeft ,
1359 Device > block_vector_type ;
1364 typedef typename matrix_type::graph_type graph_type ;
1369 block_vector_type
x =
1370 block_vector_type(
"x" ,
setup.stoch_length ,
setup.fem_length );
1371 block_vector_type
y =
1372 block_vector_type(
"y" ,
setup.stoch_length ,
setup.fem_length );
1374 typename block_vector_type::HostMirror hx =
1377 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
1378 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
1379 hx(iColStoch,iColFEM) =
1380 generate_vector_coefficient<ScalarType>(
1381 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
1389 matrix_type matrix ;
1401 const bool symmetric =
true;
1406 Stokhos::create_stochastic_product_tensor< TensorType >( *
setup.basis,
1409 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1410 std::string(
"test crs graph") ,
setup.fem_graph );
1412 matrix.values = block_vector_type(
1413 "matrix" ,
setup.stoch_length ,
setup.fem_graph_length );
1415 typename block_vector_type::HostMirror hM =
1418 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
1419 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1420 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1422 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1424 generate_matrix_coefficient<ScalarType>(
1425 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , k );
1439 bool success =
setup.test_commuted(hy, out);