42 #ifndef STOKHOS_KOKKOS_CORE_KERNELS_UNIT_TEST_HPP 43 #define STOKHOS_KOKKOS_CORE_KERNELS_UNIT_TEST_HPP 45 #include "Teuchos_TestingHelpers.hpp" 46 #include "Teuchos_UnitTestHelpers.hpp" 47 #include "Teuchos_TestForException.hpp" 51 #include "EpetraExt_BlockUtility.h" 55 #include "Epetra_MpiComm.h" 57 #include "Epetra_SerialComm.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" 87 using Teuchos::ParameterList;
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) {
194 globalComm = rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
207 Array< RCP<const abstract_basis_type> > bases(
d);
208 for (
int i=0; i<
d; i++)
212 Cijk =
basis->computeTripleProductTensor();
215 ParameterList parallelParams;
216 RCP<Stokhos::ParallelData> sg_parallel_data =
218 RCP<const EpetraExt::MultiComm> sg_comm =
219 sg_parallel_data->getMultiComm();
220 RCP<const Epetra_Comm> app_comm =
221 sg_parallel_data->getSpatialComm();
222 RCP<const Stokhos::EpetraSparse3Tensor> epetraCijk =
223 sg_parallel_data->getEpetraCijk();
224 RCP<const Epetra_BlockMap> stoch_row_map =
225 epetraCijk->getStochasticRowMap();
230 RCP<const Epetra_Map> x_map =
231 rcp(
new Epetra_Map(
fem_length, 0, *app_comm));
232 RCP<const Epetra_Map> sg_x_map =
233 rcp(EpetraExt::BlockUtility::GenerateBlockMap(
234 *x_map, *stoch_row_map, *sg_comm));
236 RCP<Epetra_CrsGraph> graph = rcp(
new Epetra_CrsGraph(Copy, *x_map, 27));
237 int *my_GIDs = x_map->MyGlobalElements();
238 int num_my_GIDs = x_map->NumMyElements();
239 for (
int i=0; i<num_my_GIDs; ++i) {
240 int row = my_GIDs[i];
243 graph->InsertGlobalIndices(row, num_indices, indices);
245 graph->FillComplete();
247 RCP<ParameterList> sg_op_params = rcp(
new ParameterList);
248 RCP<Stokhos::MatrixFreeOperator> sg_A =
250 x_map, x_map, sg_x_map, sg_x_map,
252 RCP<Epetra_BlockMap> sg_A_overlap_map =
253 rcp(
new Epetra_LocalMap(
254 stoch_length, 0, *(sg_parallel_data->getStochasticComm())));
255 RCP< Stokhos::EpetraOperatorOrthogPoly > A_sg_blocks =
257 basis, sg_A_overlap_map, x_map, x_map, sg_x_map, sg_comm));
259 RCP<Epetra_CrsMatrix> A = rcp(
new Epetra_CrsMatrix(Copy, *graph));
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);
273 sg_A->setupOperator(A_sg_blocks);
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>(
290 sg_A->Apply( *(
sg_x->getBlockVector()), *(
sg_y->getBlockVector()) );
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>
323 Teuchos::FancyOStream& out)
const {
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>
346 Teuchos::FancyOStream& out)
const {
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>
369 Teuchos::FancyOStream& out)
const {
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>
393 Teuchos::FancyOStream& out)
const {
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>
417 Teuchos::FancyOStream& out)
const {
426 out <<
"y_expected[" << block <<
"][" << i <<
"] - " 427 <<
"y(" << b <<
"," << i <<
") == " 428 << diff <<
" < " << tol <<
" : failed" 431 success = success && s;
438 template <
typename vec_type>
440 Teuchos::FancyOStream& out)
const {
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>
465 Teuchos::FancyOStream& out) {
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>
575 Teuchos::FancyOStream& out) {
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>
679 Teuchos::FancyOStream& out) {
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 >
785 Teuchos::FancyOStream& out)
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 >
896 Teuchos::FancyOStream& out)
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 ) {
912 int len = cijk_graph->NumGlobalIndices(i);
913 stoch_graph[i].resize(len);
915 cijk_graph->ExtractGlobalRowCopy(i, len, len2, &stoch_graph[i][0]);
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 >
1045 Teuchos::FancyOStream& out)
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 ) {
1061 int len = cijk_graph->NumGlobalIndices(i);
1062 stoch_graph[i].resize(len);
1064 cijk_graph->ExtractGlobalRowCopy(i, len, len2, &stoch_graph[i][0]);
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 >
1192 Teuchos::FancyOStream& out,
1193 const Teuchos::ParameterList& params = Teuchos::ParameterList()) {
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 >
1271 Teuchos::FancyOStream& out,
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 ;
1285 Teuchos::ParameterList params;
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 >
1355 Teuchos::FancyOStream& out) {
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;
1402 Teuchos::RCP< Stokhos::LTBSparse3Tensor<ordinal_type, value_type> > Cijk =
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);
RCP< const Epetra_Comm > globalComm
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
Stokhos::Sparse3Tensor< int, value_type > Cijk_type
An Epetra operator representing the block stochastic Galerkin operator.
Bases defined by combinatorial product of polynomial bases.
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
bool test_lexo_block_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
bool test_crs_product_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out, const Teuchos::ParameterList ¶ms=Teuchos::ParameterList())
RCP< Stokhos::EpetraVectorOrthogPoly > sg_y
bool test_crs_flat_commuted(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
Stokhos::OneDOrthogPolyBasis< int, value_type > abstract_basis_type
Symmetric diagonal storage for a dense matrix.
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
Stokhos::TotalOrderBasis< int, value_type, less_type > product_basis_type
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTBBlockLeaf(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
Bi-directional iterator for traversing a sparse array.
bool test_linear_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out, const bool symmetric)
void multiply(const CrsMatrix< MatrixValue, Device, Layout > &A, const InputMultiVectorType &x, OutputMultiVectorType &y, const std::vector< OrdinalType > &col_indices, SingleColumnMultivectorMultiply)
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
void setup(int p_=5, int d_=2)
A multidimensional index.
bool test_original_flat(const vec_type &y, Teuchos::FancyOStream &out) const
RCP< Stokhos::EpetraVectorOrthogPoly > sg_x
RCP< Stokhos::ProductEpetraVector > sg_y_commuted
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
bool test_commuted_flat(const vec_type &y, Teuchos::FancyOStream &out) const
A container class storing an orthogonal polynomial whose coefficients are vectors, operators, or in general any type that would have an expensive copy constructor.
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
UnitTestSetup< int, double > setup
A container class for products of Epetra_Vector's.
bool test_crs_flat_original(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
bool test_commuted(const vec_type &y, Teuchos::FancyOStream &out) const
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)
bool test_original(const multi_vec_type &y, Teuchos::FancyOStream &out) const
Stokhos::LexographicLess< Stokhos::MultiIndex< int > > less_type
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< XD, XP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< YD, YP... > >::value &&Kokkos::is_view_uq_pce< Kokkos::View< ZD, ZP... > >::value >::type update(const typename Kokkos::View< XD, XP... >::array_type::non_const_value_type &alpha, const Kokkos::View< XD, XP... > &x, const typename Kokkos::View< YD, YP... >::array_type::non_const_value_type &beta, const Kokkos::View< YD, YP... > &y, const typename Kokkos::View< ZD, ZP... >::array_type::non_const_value_type &gamma, const Kokkos::View< ZD, ZP... > &z)
Legendre polynomial basis.
bool test_original(const std::vector< vec_type > &y, Teuchos::FancyOStream &out) const
RCP< product_basis_type > basis
std::vector< std::vector< int > > fem_graph
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
Abstract base class for 1-D orthogonal polynomials.
Teuchos::Array< int > inv_perm
Teuchos::RCP< Epetra_CrsGraph > sparse3Tensor2CrsGraph(const Stokhos::OrthogPolyBasis< ordinal_type, value_type > &basis, const Stokhos::Sparse3Tensor< ordinal_type, value_type > &Cijk, const Epetra_Comm &comm)
Build an Epetra_CrsGraph from a sparse 3 tensor.
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
CRS matrix of dense blocks.
A comparison functor implementing a strict weak ordering based lexographic ordering.
Data structure storing a dense 3-tensor C(i,j,k).
Stokhos::LegendreBasis< int, value_type > basis_type
Sacado::MP::Vector< storage_type > vec_type
Teuchos::Array< int > perm
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
void update(const ValueType &alpha, VectorType &x, const ValueType &beta, const VectorType &y)
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)
bool test_crs_matrix_free(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
bool test_crs_matrix_free_view(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
bool test_crs_dense_block(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
bool test_commuted_perm(const vec_type &y, Teuchos::FancyOStream &out) const