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 "Kokkos_Sparse.hpp" 78 #include "Kokkos_Blas1_MV.hpp" 86 using Teuchos::ParameterList;
88 template<
typename IntType >
95 return k + N * (
j + N * i );
98 template <
typename ordinal >
101 std::vector< std::vector<ordinal> > & graph )
103 graph.resize( N * N * N , std::vector<ordinal>() );
107 for (
int i = 0 ; i < (
int) N ; ++i ) {
108 for (
int j = 0 ;
j < (
int) N ; ++
j ) {
109 for (
int k = 0 ; k < (
int) N ; ++k ) {
113 graph[row].reserve(27);
115 for (
int ii = -1 ; ii < 2 ; ++ii ) {
116 for (
int jj = -1 ; jj < 2 ; ++jj ) {
117 for (
int kk = -1 ; kk < 2 ; ++kk ) {
118 if ( 0 <= i + ii && i + ii < (
int) N &&
119 0 <=
j + jj &&
j + jj < (
int) N &&
120 0 <= k + kk && k + kk < (
int) N ) {
123 graph[row].push_back(col);
126 total += graph[row].size();
132 template <
typename scalar,
typename ordinal>
135 const ordinal nStoch ,
136 const ordinal iRowFEM ,
137 const ordinal iColFEM ,
138 const ordinal iStoch )
140 const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
141 ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
143 const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
145 return A_fem + A_stoch ;
149 template <
typename scalar,
typename ordinal>
152 const ordinal nStoch ,
153 const ordinal iColFEM ,
154 const ordinal iStoch )
156 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
157 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
158 return X_fem + X_stoch ;
162 template <
typename Device>
183 void setup(
int p_ = 5,
int d_ = 2) {
193 globalComm = rcp(
new Epetra_MpiComm(MPI_COMM_WORLD));
206 Array< RCP<const abstract_basis_type> > bases(
d);
207 for (
int i=0; i<
d; i++)
211 Cijk =
basis->computeTripleProductTensor();
214 ParameterList parallelParams;
215 RCP<Stokhos::ParallelData> sg_parallel_data =
217 RCP<const EpetraExt::MultiComm> sg_comm =
218 sg_parallel_data->getMultiComm();
219 RCP<const Epetra_Comm> app_comm =
220 sg_parallel_data->getSpatialComm();
221 RCP<const Stokhos::EpetraSparse3Tensor> epetraCijk =
222 sg_parallel_data->getEpetraCijk();
223 RCP<const Epetra_BlockMap> stoch_row_map =
224 epetraCijk->getStochasticRowMap();
229 RCP<const Epetra_Map> x_map =
230 rcp(
new Epetra_Map(
fem_length, 0, *app_comm));
231 RCP<const Epetra_Map> sg_x_map =
232 rcp(EpetraExt::BlockUtility::GenerateBlockMap(
233 *x_map, *stoch_row_map, *sg_comm));
235 RCP<Epetra_CrsGraph> graph = rcp(
new Epetra_CrsGraph(Copy, *x_map, 27));
236 int *my_GIDs = x_map->MyGlobalElements();
237 int num_my_GIDs = x_map->NumMyElements();
238 for (
int i=0; i<num_my_GIDs; ++i) {
239 int row = my_GIDs[i];
242 graph->InsertGlobalIndices(row, num_indices, indices);
244 graph->FillComplete();
246 RCP<ParameterList> sg_op_params = rcp(
new ParameterList);
247 RCP<Stokhos::MatrixFreeOperator> sg_A =
249 x_map, x_map, sg_x_map, sg_x_map,
251 RCP<Epetra_BlockMap> sg_A_overlap_map =
252 rcp(
new Epetra_LocalMap(
253 stoch_length, 0, *(sg_parallel_data->getStochasticComm())));
254 RCP< Stokhos::EpetraOperatorOrthogPoly > A_sg_blocks =
256 basis, sg_A_overlap_map, x_map, x_map, sg_x_map, sg_comm));
258 RCP<Epetra_CrsMatrix> A = rcp(
new Epetra_CrsMatrix(Copy, *graph));
260 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
fem_length ; ++iRowFEM ) {
261 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
262 const int iColFEM =
fem_graph[iRowFEM][iRowEntryFEM] ;
264 double v = generate_matrix_coefficient<double>(
266 A->ReplaceGlobalValues(iRowFEM, 1, &v, &iColFEM);
270 A_sg_blocks->setCoeffPtr(i, A);
272 sg_A->setupOperator(A_sg_blocks);
275 basis, stoch_row_map, x_map, sg_x_map, sg_comm));
277 basis, stoch_row_map, x_map, sg_x_map, sg_comm));
280 for (
int iColFEM=0; iColFEM <
fem_length; ++iColFEM ) {
281 for (
int iColStoch=0 ; iColStoch <
stoch_length; ++iColStoch ) {
282 (*sg_x)[iColStoch][iColFEM] = generate_vector_coefficient<double>(
289 sg_A->Apply( *(
sg_x->getBlockVector()), *(
sg_y->getBlockVector()) );
293 x_map, stoch_row_map, sg_comm));
299 typedef typename Kokkos::ViewTraits<value_type,Device,void,void>::host_mirror_space host_device;
302 host_device > stoch_tensor_type ;
304 stoch_tensor_type tensor =
305 Stokhos::create_stochastic_product_tensor< tensor_type >( *
basis, *
Cijk );
312 for (
int j=0;
j<
d; ++
j)
313 term[
j] = tensor.bases_degree(i,
j);
314 int idx =
basis->index(term);
320 template <
typename vec_type>
322 Teuchos::FancyOStream& out)
const {
330 out <<
"y_expected[" << block <<
"][" << i <<
"] - " 331 <<
"y[" << block <<
"][" << i <<
"] = " << (*sg_y)[block][i]
332 <<
" - " <<
y[block][i] <<
" == " 333 << diff <<
" < " << tol <<
" : failed" 336 success = success && s;
343 template <
typename multi_vec_type>
345 Teuchos::FancyOStream& out)
const {
353 out <<
"y_expected[" << block <<
"][" << i <<
"] - " 354 <<
"y(" << i <<
"," << block <<
") = " << (*sg_y)[block][i]
355 <<
" - " <<
y(i,block) <<
" == " 356 << diff <<
" < " << tol <<
" : failed" 359 success = success && s;
366 template <
typename vec_type>
368 Teuchos::FancyOStream& out)
const {
377 out <<
"y_expected[" << block <<
"][" << i <<
"] - " 378 <<
"y(" << b <<
"," << i <<
") = " << (*sg_y)[block][i]
379 <<
" - " <<
y(b,i) <<
" == " 380 << diff <<
" < " << tol <<
" : failed" 383 success = success && s;
390 template <
typename vec_type>
392 Teuchos::FancyOStream& out)
const {
401 out <<
"y_expected[" << block <<
"][" << i <<
"] - " 402 <<
"y(" << b <<
"," << i <<
") = " << (*sg_y)[block][i] <<
" - " 404 << diff <<
" < " << tol <<
" : failed" 407 success = success && s;
414 template <
typename vec_type>
416 Teuchos::FancyOStream& out)
const {
425 out <<
"y_expected[" << block <<
"][" << i <<
"] - " 426 <<
"y(" << b <<
"," << i <<
") == " 427 << diff <<
" < " << tol <<
" : failed" 430 success = success && s;
437 template <
typename vec_type>
439 Teuchos::FancyOStream& out)
const {
448 out <<
"y_expected[" << block <<
"][" << i <<
"] - " 449 <<
"y(" << b <<
"," << i <<
") == " 450 << diff <<
" < " << tol <<
" : failed" 453 success = success && s;
462 template <
typename value_type,
typename Device,
typename SparseMatOps>
464 Teuchos::FancyOStream& out) {
466 typedef typename matrix_type::values_type matrix_values_type;
467 typedef typename matrix_type::graph_type matrix_graph_type;
468 typedef Kokkos::View<value_type*,Device>
vec_type;
472 std::vector<matrix_type> matrix(
setup.stoch_length ) ;
473 std::vector<vec_type>
x(
setup.stoch_length ) ;
474 std::vector<vec_type>
y(
setup.stoch_length ) ;
475 std::vector<vec_type> tmp(
setup.stoch_length ) ;
477 for (
int block=0; block<
setup.stoch_length; ++block) {
478 matrix[block].graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
479 std::string(
"testing") ,
setup.fem_graph );
481 matrix[block].values =
482 matrix_values_type(
"matrix" ,
setup.fem_graph_length );
488 typename matrix_values_type::HostMirror hM =
491 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
492 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
493 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
495 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
496 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , block );
502 typename vec_type::HostMirror hx =
504 typename vec_type::HostMirror hy =
507 for (
int i = 0 ; i <
setup.fem_length ; ++i ) {
508 hx(i) = generate_vector_coefficient<value_type>(
509 setup.fem_length ,
setup.stoch_length , i , block );
520 setup.Cijk->k_begin();
524 k_it!=k_end; ++k_it) {
525 int nj =
setup.Cijk->num_j(k_it);
529 setup.Cijk->j_begin(k_it);
531 setup.Cijk->j_end(k_it);
532 std::vector<vec_type> xx(nj), yy(nj);
545 j_it != j_end; ++j_it) {
547 setup.Cijk->i_begin(j_it);
549 setup.Cijk->i_end(j_it);
562 std::vector<typename vec_type::HostMirror> hy(
setup.stoch_length);
563 for (
int i=0; i<
setup.stoch_length; ++i) {
567 bool success =
setup.test_original(hy, out);
572 template <
typename value_type,
typename Device,
typename SparseMatOps>
574 Teuchos::FancyOStream& out) {
576 typedef typename matrix_type::values_type matrix_values_type;
577 typedef typename matrix_type::graph_type matrix_graph_type;
578 typedef Kokkos::View<value_type*, Kokkos::LayoutLeft, Device, Kokkos::MemoryUnmanaged>
vec_type;
579 typedef Kokkos::View<value_type**, Kokkos::LayoutLeft, Device> multi_vec_type;
583 std::vector<matrix_type> matrix(
setup.stoch_length ) ;
584 multi_vec_type
x(
"x",
setup.fem_length,
setup.stoch_length ) ;
585 multi_vec_type
y(
"y",
setup.fem_length,
setup.stoch_length ) ;
586 multi_vec_type tmp_x(
"tmp_x",
setup.fem_length,
setup.stoch_length ) ;
587 multi_vec_type tmp_y(
"tmp_y",
setup.fem_length,
setup.stoch_length ) ;
592 for (
int block=0; block<
setup.stoch_length; ++block) {
593 matrix[block].graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
594 std::string(
"testing") ,
setup.fem_graph );
596 matrix[block].values =
597 matrix_values_type(
"matrix" ,
setup.fem_graph_length );
599 typename matrix_values_type::HostMirror hM =
602 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
603 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
604 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
606 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
607 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , block );
613 for (
int i = 0 ; i <
setup.fem_length ; ++i ) {
614 hx(i, block) = generate_vector_coefficient<value_type>(
615 setup.fem_length ,
setup.stoch_length , i , block );
629 k_iterator k_begin =
setup.Cijk->k_begin();
630 k_iterator k_end =
setup.Cijk->k_end();
631 for (k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
632 unsigned nj =
setup.Cijk->num_j(k_it);
635 kj_iterator j_begin =
setup.Cijk->j_begin(k_it);
636 kj_iterator j_end =
setup.Cijk->j_end(k_it);
637 std::vector<int> j_indices(nj);
639 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
641 vec_type xx = Kokkos::subview(
x, Kokkos::ALL(),
j );
642 vec_type tt = Kokkos::subview( tmp_x, Kokkos::ALL(), jdx++ );
645 multi_vec_type tmp_x_view =
646 Kokkos::subview( tmp_x, Kokkos::ALL(),
647 std::make_pair(0u,nj));
648 multi_vec_type tmp_y_view =
649 Kokkos::subview( tmp_y, Kokkos::ALL(),
650 std::make_pair(0u,nj));
653 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
655 Kokkos::subview( tmp_y, Kokkos::ALL(), jdx++ );
656 kji_iterator i_begin =
setup.Cijk->i_begin(j_it);
657 kji_iterator i_end =
setup.Cijk->i_end(j_it);
658 for (kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
661 vec_type y_view = Kokkos::subview(
y, Kokkos::ALL(), i );
669 bool success =
setup.test_original(hy, out);
674 #ifdef HAVE_STOKHOS_KOKKOSLINALG 676 template <
typename value_type,
typename Device>
678 Teuchos::FancyOStream& out) {
680 typedef KokkosSparse::CrsMatrix<value_type,ordinal_type,Device> matrix_type;
681 typedef typename matrix_type::values_type matrix_values_type;
682 typedef typename matrix_type::StaticCrsGraphType matrix_graph_type;
683 typedef Kokkos::View<value_type*, Kokkos::LayoutLeft, Device, Kokkos::MemoryUnmanaged>
vec_type;
684 typedef Kokkos::View<value_type**, Kokkos::LayoutLeft, Device> multi_vec_type;
688 std::vector<matrix_type> matrix(
setup.stoch_length ) ;
689 multi_vec_type
x(
"x",
setup.fem_length,
setup.stoch_length ) ;
690 multi_vec_type
y(
"y",
setup.fem_length,
setup.stoch_length ) ;
691 multi_vec_type tmp_x(
"tmp_x",
setup.fem_length,
setup.stoch_length ) ;
692 multi_vec_type tmp_y(
"tmp_y",
setup.fem_length,
setup.stoch_length ) ;
697 for (
int block=0; block<
setup.stoch_length; ++block) {
698 matrix_graph_type matrix_graph =
699 Kokkos::create_staticcrsgraph<matrix_graph_type>(
700 std::string(
"test crs graph"),
setup.fem_graph);
701 matrix_values_type matrix_values =
702 matrix_values_type(
"matrix" ,
setup.fem_graph_length );
703 typename matrix_values_type::HostMirror hM =
706 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
707 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
708 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
710 hM(iEntryFEM) = generate_matrix_coefficient<value_type>(
711 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , block );
716 matrix[block] = matrix_type(
"matrix",
setup.fem_length, matrix_values,
719 for (
int i = 0 ; i <
setup.fem_length ; ++i ) {
720 hx(i, block) = generate_vector_coefficient<value_type>(
721 setup.fem_length ,
setup.stoch_length , i , block );
734 k_iterator k_begin =
setup.Cijk->k_begin();
735 k_iterator k_end =
setup.Cijk->k_end();
736 for (k_iterator k_it=k_begin; k_it!=k_end; ++k_it) {
737 int nj =
setup.Cijk->num_j(k_it);
740 kj_iterator j_begin =
setup.Cijk->j_begin(k_it);
741 kj_iterator j_end =
setup.Cijk->j_end(k_it);
743 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
745 vec_type xx = Kokkos::subview(
x, Kokkos::ALL(),
j );
746 vec_type tt = Kokkos::subview( tmp_x, Kokkos::ALL(), jdx++ );
749 multi_vec_type tmp_x_view =
750 Kokkos::subview( tmp_x, Kokkos::ALL(),
751 std::make_pair(0u,jdx));
752 multi_vec_type tmp_y_view =
753 Kokkos::subview( tmp_y, Kokkos::ALL(),
754 std::make_pair(0u,jdx));
757 for (kj_iterator j_it = j_begin; j_it != j_end; ++j_it) {
759 Kokkos::subview( tmp_y, Kokkos::ALL(), jdx++ );
760 kji_iterator i_begin =
setup.Cijk->i_begin(j_it);
761 kji_iterator i_end =
setup.Cijk->i_end(j_it);
762 for (kji_iterator i_it = i_begin; i_it != i_end; ++i_it) {
765 vec_type y_view = Kokkos::subview(
y, Kokkos::ALL(), i );
774 bool success =
setup.test_original(hy, out);
781 template<
typename ScalarType ,
class Device >
784 Teuchos::FancyOStream& out)
787 typedef Kokkos::View< value_type**, Kokkos::LayoutLeft , Device > block_vector_type ;
793 typedef typename matrix_type::graph_type graph_type ;
800 for (k_iterator k_it=
setup.Cijk->k_begin();
801 k_it!=
setup.Cijk->k_end(); ++k_it) {
803 for (kj_iterator j_it =
setup.Cijk->j_begin(k_it);
804 j_it !=
setup.Cijk->j_end(k_it); ++j_it) {
806 for (kji_iterator i_it =
setup.Cijk->i_begin(j_it);
807 i_it !=
setup.Cijk->i_end(j_it); ++i_it) {
809 double c = value(i_it);
810 dense_cijk(i,
j,k) = c;
818 block_vector_type
x =
819 block_vector_type(
"x" ,
setup.stoch_length ,
setup.fem_length );
820 block_vector_type
y =
821 block_vector_type(
"y" ,
setup.stoch_length ,
setup.fem_length );
825 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
826 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
827 hx(iColStoch,iColFEM) =
828 generate_vector_coefficient<ScalarType>(
829 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
842 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
843 std::string(
"test crs graph") ,
setup.fem_graph );
844 matrix.values = block_vector_type(
845 "matrix" , matrix.block.matrix_size() ,
setup.fem_graph_length );
848 typename block_vector_type::HostMirror hM =
851 for (
int iRowStoch = 0 ; iRowStoch <
setup.stoch_length ; ++iRowStoch ) {
852 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
854 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
855 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM];
857 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
859 const size_t offset =
860 matrix.block.matrix_offset( iRowStoch , iColStoch );
862 ScalarType value = 0 ;
864 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
865 value += dense_cijk( iRowStoch , iColStoch , k ) *
866 generate_matrix_coefficient<ScalarType>(
867 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , k );
870 hM( offset , iEntryFEM ) = value ;
887 bool success =
setup.test_commuted(hy, out);
892 template<
typename ScalarType ,
class Device >
895 Teuchos::FancyOStream& out)
898 typedef Kokkos::View< value_type* , Device > vector_type ;
903 typedef typename matrix_type::values_type matrix_values_type;
904 typedef typename matrix_type::graph_type matrix_graph_type;
907 std::vector< std::vector< int > > stoch_graph(
setup.stoch_length );
910 for (
int i = 0 ; i <
setup.stoch_length ; ++i ) {
911 int len = cijk_graph->NumGlobalIndices(i);
912 stoch_graph[i].resize(len);
914 cijk_graph->ExtractGlobalRowCopy(i, len, len2, &stoch_graph[i][0]);
922 for (k_iterator k_it=
setup.Cijk->k_begin();
923 k_it!=
setup.Cijk->k_end(); ++k_it) {
925 for (kj_iterator j_it =
setup.Cijk->j_begin(k_it);
926 j_it !=
setup.Cijk->j_end(k_it); ++j_it) {
928 for (kji_iterator i_it =
setup.Cijk->i_begin(j_it);
929 i_it !=
setup.Cijk->i_end(j_it); ++i_it) {
931 double c = value(i_it);
932 dense_cijk(i,
j,k) = c;
940 const int flat_length =
setup.fem_length *
setup.stoch_length ;
942 std::vector< std::vector<int> > flat_graph( flat_length );
944 for (
int iOuterRow = 0 ; iOuterRow <
setup.fem_length ; ++iOuterRow ) {
946 const size_t iOuterRowNZ =
setup.fem_graph[iOuterRow].size();
948 for (
int iInnerRow = 0 ; iInnerRow <
setup.stoch_length ; ++iInnerRow ) {
950 const size_t iInnerRowNZ = stoch_graph[ iInnerRow ].size(); ;
951 const int iFlatRowNZ = iOuterRowNZ * iInnerRowNZ ;
952 const int iFlatRow = iInnerRow + iOuterRow *
setup.stoch_length ;
954 flat_graph[iFlatRow].resize( iFlatRowNZ );
958 for (
size_t iOuterEntry = 0 ; iOuterEntry < iOuterRowNZ ; ++iOuterEntry ) {
960 const int iOuterCol =
setup.fem_graph[iOuterRow][iOuterEntry];
962 for (
size_t iInnerEntry = 0 ; iInnerEntry < iInnerRowNZ ; ++iInnerEntry ) {
964 const int iInnerCol = stoch_graph[iInnerRow][iInnerEntry] ;
965 const int iFlatColumn = iInnerCol + iOuterCol *
setup.stoch_length ;
967 flat_graph[iFlatRow][iFlatEntry] = iFlatColumn ;
977 vector_type
x = vector_type(
"x" , flat_length );
978 vector_type
y = vector_type(
"y" , flat_length );
982 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
983 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
984 hx(iColStoch + iColFEM*
setup.stoch_length) =
985 generate_vector_coefficient<ScalarType>(
986 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
996 matrix.graph = Kokkos::create_staticcrsgraph<matrix_graph_type>(
997 std::string(
"testing") , flat_graph );
999 const size_t flat_graph_length = matrix.graph.entries.dimension_0();
1001 matrix.values = matrix_values_type(
"matrix" , flat_graph_length );
1003 typename matrix_values_type::HostMirror hM =
1006 for (
int iRow = 0 , iEntry = 0 ; iRow < flat_length ; ++iRow ) {
1007 const int iRowFEM = iRow /
setup.stoch_length ;
1008 const int iRowStoch = iRow %
setup.stoch_length ;
1010 for (
size_t iRowEntry = 0 ; iRowEntry < flat_graph[ iRow ].size() ; ++iRowEntry , ++iEntry ) {
1011 const int iCol = flat_graph[ iRow ][ iRowEntry ];
1012 const int iColFEM = iCol /
setup.stoch_length ;
1013 const int iColStoch = iCol %
setup.stoch_length ;
1016 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1017 const double A_fem_k =
1018 generate_matrix_coefficient<ScalarType>(
1019 setup.fem_length ,
setup.stoch_length , iRowFEM, iColFEM, k );
1020 value += dense_cijk(iRowStoch,iColStoch,k) * A_fem_k ;
1022 hM( iEntry ) = value ;
1037 bool success =
setup.test_commuted_flat(hy, out);
1041 template<
typename ScalarType ,
class Device >
1044 Teuchos::FancyOStream& out)
1047 typedef Kokkos::View< value_type* , Device > vector_type ;
1052 typedef typename matrix_type::values_type matrix_values_type;
1053 typedef typename matrix_type::graph_type matrix_graph_type;
1056 std::vector< std::vector< int > > stoch_graph(
setup.stoch_length );
1059 for (
int i = 0 ; i <
setup.stoch_length ; ++i ) {
1060 int len = cijk_graph->NumGlobalIndices(i);
1061 stoch_graph[i].resize(len);
1063 cijk_graph->ExtractGlobalRowCopy(i, len, len2, &stoch_graph[i][0]);
1071 for (k_iterator k_it=
setup.Cijk->k_begin();
1072 k_it!=
setup.Cijk->k_end(); ++k_it) {
1073 int k = index(k_it);
1074 for (kj_iterator j_it =
setup.Cijk->j_begin(k_it);
1075 j_it !=
setup.Cijk->j_end(k_it); ++j_it) {
1076 int j = index(j_it);
1077 for (kji_iterator i_it =
setup.Cijk->i_begin(j_it);
1078 i_it !=
setup.Cijk->i_end(j_it); ++i_it) {
1079 int i = index(i_it);
1080 double c = value(i_it);
1081 dense_cijk(i,
j,k) = c;
1089 const size_t flat_length =
setup.fem_length *
setup.stoch_length ;
1091 std::vector< std::vector<int> > flat_graph( flat_length );
1093 for (
int iOuterRow = 0 ; iOuterRow <
setup.stoch_length ; ++iOuterRow ) {
1095 const size_t iOuterRowNZ = stoch_graph[iOuterRow].size();
1097 for (
int iInnerRow = 0 ; iInnerRow <
setup.fem_length ; ++iInnerRow ) {
1099 const size_t iInnerRowNZ =
setup.fem_graph[iInnerRow].size();
1100 const int iFlatRowNZ = iOuterRowNZ * iInnerRowNZ ;
1101 const int iFlatRow = iInnerRow + iOuterRow *
setup.fem_length ;
1103 flat_graph[iFlatRow].resize( iFlatRowNZ );
1105 int iFlatEntry = 0 ;
1107 for (
size_t iOuterEntry = 0 ; iOuterEntry < iOuterRowNZ ; ++iOuterEntry ) {
1109 const int iOuterCol = stoch_graph[ iOuterRow ][ iOuterEntry ];
1111 for (
size_t iInnerEntry = 0 ; iInnerEntry < iInnerRowNZ ; ++iInnerEntry ) {
1113 const int iInnerCol =
setup.fem_graph[ iInnerRow][iInnerEntry];
1114 const int iFlatColumn = iInnerCol + iOuterCol *
setup.fem_length ;
1116 flat_graph[iFlatRow][iFlatEntry] = iFlatColumn ;
1125 vector_type
x = vector_type(
"x" , flat_length );
1126 vector_type
y = vector_type(
"y" , flat_length );
1130 for (
size_t iCol = 0 ; iCol < flat_length ; ++iCol ) {
1131 const int iColStoch = iCol /
setup.fem_length ;
1132 const int iColFEM = iCol %
setup.fem_length ;
1134 hx(iCol) = generate_vector_coefficient<ScalarType>(
1135 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
1142 matrix_type matrix ;
1144 matrix.graph = Kokkos::create_staticcrsgraph<matrix_graph_type>( std::string(
"testing") , flat_graph );
1146 const size_t flat_graph_length = matrix.graph.entries.dimension_0();
1148 matrix.values = matrix_values_type(
"matrix" , flat_graph_length );
1150 typename matrix_values_type::HostMirror hM =
1153 for (
size_t iRow = 0 , iEntry = 0 ; iRow < flat_length ; ++iRow ) {
1154 const int iRowStoch = iRow /
setup.fem_length ;
1155 const int iRowFEM = iRow %
setup.fem_length ;
1157 for (
size_t iRowEntry = 0 ; iRowEntry < flat_graph[ iRow ].size() ; ++iRowEntry , ++iEntry ) {
1158 const int iCol = flat_graph[ iRow ][ iRowEntry ];
1159 const int iColStoch = iCol /
setup.fem_length ;
1160 const int iColFEM = iCol %
setup.fem_length ;
1163 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1164 const double A_fem_k =
1165 generate_matrix_coefficient<ScalarType>(
1167 iRowFEM , iColFEM , k );
1168 value += dense_cijk(iRowStoch,iColStoch,k) * A_fem_k ;
1170 hM( iEntry ) = value ;
1184 bool success =
setup.test_original_flat(hy, out);
1188 template<
typename ScalarType ,
typename TensorType,
class Device >
1191 Teuchos::FancyOStream& out,
1192 const Teuchos::ParameterList& params = Teuchos::ParameterList()) {
1194 typedef Kokkos::View<
value_type** , Kokkos::LayoutLeft ,
1195 Device > block_vector_type ;
1200 typedef typename matrix_type::graph_type graph_type ;
1205 block_vector_type
x =
1206 block_vector_type(
"x" ,
setup.stoch_length_aligned ,
setup.fem_length );
1207 block_vector_type
y =
1208 block_vector_type(
"y" ,
setup.stoch_length_aligned ,
setup.fem_length );
1210 typename block_vector_type::HostMirror hx =
1213 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
1214 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
1215 hx(
setup.perm[iColStoch],iColFEM) =
1216 generate_vector_coefficient<ScalarType>(
1217 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
1225 matrix_type matrix ;
1228 Stokhos::create_stochastic_product_tensor< TensorType >( *
setup.basis,
1232 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1233 std::string(
"test crs graph") ,
setup.fem_graph );
1235 matrix.values = block_vector_type(
1236 "matrix" ,
setup.stoch_length_aligned ,
setup.fem_graph_length );
1238 typename block_vector_type::HostMirror hM =
1241 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
1242 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1243 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1245 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1246 hM(
setup.perm[k],iEntryFEM) =
1247 generate_matrix_coefficient<ScalarType>(
1248 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , k );
1263 bool success =
setup.test_commuted_perm(hy, out);
1268 template<
typename ScalarType ,
class Device ,
int BlockSize >
1270 Teuchos::FancyOStream& out,
1271 const bool symmetric) {
1273 typedef Kokkos::View<
value_type** , Kokkos::LayoutLeft ,
1274 Device > block_vector_type ;
1279 typedef typename matrix_type::graph_type graph_type ;
1282 matrix_type matrix ;
1284 Teuchos::ParameterList params;
1285 params.set(
"Symmetric",symmetric);
1287 Stokhos::create_stochastic_product_tensor< TensorType >( *
setup.basis,
1290 int aligned_stoch_length = matrix.block.tensor().aligned_dimension();
1295 block_vector_type
x =
1296 block_vector_type(
"x" , aligned_stoch_length ,
setup.fem_length );
1297 block_vector_type
y =
1298 block_vector_type(
"y" , aligned_stoch_length ,
setup.fem_length );
1300 typename block_vector_type::HostMirror hx =
1303 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
1304 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
1305 hx(iColStoch,iColFEM) =
1306 generate_vector_coefficient<ScalarType>(
1307 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
1316 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1317 std::string(
"test crs graph") ,
setup.fem_graph );
1319 matrix.values = block_vector_type(
1320 "matrix" , aligned_stoch_length ,
setup.fem_graph_length );
1322 typename block_vector_type::HostMirror hM =
1325 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
1326 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1327 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1329 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1331 generate_matrix_coefficient<ScalarType>(
1332 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , k );
1347 bool success =
setup.test_commuted(hy, out);
1352 template<
typename ScalarType ,
class Device >
1354 Teuchos::FancyOStream& out) {
1357 typedef Kokkos::View<
value_type** , Kokkos::LayoutLeft ,
1358 Device > block_vector_type ;
1363 typedef typename matrix_type::graph_type graph_type ;
1368 block_vector_type
x =
1369 block_vector_type(
"x" ,
setup.stoch_length ,
setup.fem_length );
1370 block_vector_type
y =
1371 block_vector_type(
"y" ,
setup.stoch_length ,
setup.fem_length );
1373 typename block_vector_type::HostMirror hx =
1376 for (
int iColFEM = 0 ; iColFEM <
setup.fem_length ; ++iColFEM ) {
1377 for (
int iColStoch = 0 ; iColStoch <
setup.stoch_length ; ++iColStoch ) {
1378 hx(iColStoch,iColFEM) =
1379 generate_vector_coefficient<ScalarType>(
1380 setup.fem_length ,
setup.stoch_length , iColFEM , iColStoch );
1388 matrix_type matrix ;
1400 const bool symmetric =
true;
1401 Teuchos::RCP< Stokhos::LTBSparse3Tensor<ordinal_type, value_type> > Cijk =
1405 Stokhos::create_stochastic_product_tensor< TensorType >( *
setup.basis,
1408 matrix.graph = Kokkos::create_staticcrsgraph<graph_type>(
1409 std::string(
"test crs graph") ,
setup.fem_graph );
1411 matrix.values = block_vector_type(
1412 "matrix" ,
setup.stoch_length ,
setup.fem_graph_length );
1414 typename block_vector_type::HostMirror hM =
1417 for (
int iRowFEM = 0 , iEntryFEM = 0 ; iRowFEM <
setup.fem_length ; ++iRowFEM ) {
1418 for (
size_t iRowEntryFEM = 0 ; iRowEntryFEM <
setup.fem_graph[iRowFEM].size() ; ++iRowEntryFEM , ++iEntryFEM ) {
1419 const int iColFEM =
setup.fem_graph[iRowFEM][iRowEntryFEM] ;
1421 for (
int k = 0 ; k <
setup.stoch_length ; ++k ) {
1423 generate_matrix_coefficient<ScalarType>(
1424 setup.fem_length ,
setup.stoch_length , iRowFEM , iColFEM , k );
1438 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