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"
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"
87using Teuchos::ParameterList;
89template<
typename IntType >
96 return k + N * (
j + N * i );
99template <
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();
133template <
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 ;
150template <
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 ;
163template <
typename Device>
184 void setup(
int p_ = 5,
int d_ = 2) {
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 =
232 RCP<const Epetra_Map> sg_x_map =
233 rcp(EpetraExt::BlockUtility::GenerateBlockMap(
234 *x_map, *stoch_row_map, *sg_comm));
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 =
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));
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 {
327 double diff = std::abs( (*
sg_y)[block][i] - y[block][i] );
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 {
350 double diff = std::abs( (*
sg_y)[block][i] - y(i,block) );
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 {
374 double diff = std::abs( (*
sg_y)[block][i] - y(b,i) );
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 {
398 double diff = std::abs( (*
sg_y)[block][i] - y(b,i) );
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;
463template <
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 );
485 x[block] = vec_type(
"x" ,
setup.fem_length );
486 y[block] = vec_type(
"y" ,
setup.fem_length );
487 tmp[block] = vec_type(
"tmp" ,
setup.fem_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);
555 value_type c = value(i_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);
573template <
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) {
655 vec_type tmp_y_view =
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) {
661 value_type c = value(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
677template <
typename value_type,
typename Device>
678bool test_crs_matrix_free_kokkos(
const UnitTestSetup<Device>&
setup,
679 Teuchos::FancyOStream& out) {
680 typedef int ordinal_type;
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);
782template<
typename ScalarType ,
class Device >
785 Teuchos::FancyOStream& out)
787 typedef ScalarType value_type ;
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);
893template<
typename ScalarType ,
class Device >
896 Teuchos::FancyOStream& out)
898 typedef ScalarType value_type ;
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);
1042template<
typename ScalarType ,
class Device >
1045 Teuchos::FancyOStream& out)
1047 typedef ScalarType value_type ;
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);
1189template<
typename ScalarType ,
typename TensorType,
class Device >
1191 const UnitTestSetup<Device>&
setup,
1192 Teuchos::FancyOStream& out,
1193 const Teuchos::ParameterList& params = Teuchos::ParameterList()) {
1194 typedef ScalarType value_type ;
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);
1269template<
typename ScalarType ,
class Device ,
int BlockSize >
1271 Teuchos::FancyOStream& out,
1272 const bool symmetric) {
1273 typedef ScalarType value_type ;
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);
1353template<
typename ScalarType ,
class Device >
1355 Teuchos::FancyOStream& out) {
1356 typedef ScalarType value_type ;
1357 typedef int ordinal_type;
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);
UnitTestSetup< Kokkos::Cuda > setup
CRS matrix of dense blocks.
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
Data structure storing a dense 3-tensor C(i,j,k).
A container class storing an orthogonal polynomial whose coefficients are vectors,...
A container class storing an orthogonal polynomial whose coefficients are vectors,...
Legendre polynomial basis.
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
A comparison functor implementing a strict weak ordering based lexographic ordering.
Sparse product tensor with replicated entries to provide subsets with a given coordinate.
An Epetra operator representing the block stochastic Galerkin operator.
A multidimensional index.
Abstract base class for 1-D orthogonal polynomials.
A container class for products of Epetra_Vector's.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
kji_sparse_array::const_iterator k_iterator
Iterator for looping over k entries.
j_sparse_array::const_iterator kji_iterator
Iterator for looping over i entries given k and j.
Bases defined by combinatorial product of polynomial bases.
Symmetric diagonal storage for a dense matrix.
Multivariate orthogonal polynomial basis generated from a total order tensor product of univariate po...
Sacado::MP::Vector< storage_type > vec_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)
bool test_crs_flat_original(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
ordinal generate_fem_graph(ordinal N, std::vector< std::vector< ordinal > > &graph)
bool test_crs_product_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out, const Teuchos::ParameterList ¶ms=Teuchos::ParameterList())
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_crs_flat_commuted(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
bool test_linear_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out, const bool symmetric)
bool test_crs_matrix_free(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
IntType map_fem_graph_coord(const IntType &N, const IntType &i, const IntType &j, const IntType &k)
bool test_lexo_block_tensor(const UnitTestSetup< Device > &setup, Teuchos::FancyOStream &out)
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
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)
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)
Teuchos::RCP< LTBSparse3Tensor< ordinal_type, value_type > > computeTripleProductTensorLTBBlockLeaf(const TotalOrderBasis< ordinal_type, value_type, LexographicLess< MultiIndex< ordinal_type > > > &product_basis, bool symmetric=false)
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.
void update(const ValueType &alpha, VectorType &x, const ValueType &beta, const VectorType &y)
void multiply(const CrsMatrix< MatrixValue, Device, Layout > &A, const InputMultiVectorType &x, OutputMultiVectorType &y, const std::vector< OrdinalType > &col_indices, SingleColumnMultivectorMultiply)
Stokhos::TotalOrderBasis< int, value_type, less_type > product_basis_type
bool test_original(const multi_vec_type &y, Teuchos::FancyOStream &out) const
bool test_commuted_perm(const vec_type &y, Teuchos::FancyOStream &out) const
bool test_commuted_flat(const vec_type &y, Teuchos::FancyOStream &out) const
Stokhos::Sparse3Tensor< int, value_type > Cijk_type
RCP< Stokhos::EpetraVectorOrthogPoly > sg_y
void setup(int p_=5, int d_=2)
std::vector< std::vector< int > > fem_graph
Teuchos::Array< int > perm
bool test_original(const std::vector< vec_type > &y, Teuchos::FancyOStream &out) const
RCP< Stokhos::ProductEpetraVector > sg_y_commuted
RCP< const Epetra_Comm > globalComm
bool test_original_flat(const vec_type &y, Teuchos::FancyOStream &out) const
bool test_commuted(const vec_type &y, Teuchos::FancyOStream &out) const
Stokhos::LexographicLess< Stokhos::MultiIndex< int > > less_type
RCP< product_basis_type > basis
Stokhos::LegendreBasis< int, value_type > basis_type
RCP< Stokhos::EpetraVectorOrthogPoly > sg_x
Stokhos::OneDOrthogPolyBasis< int, value_type > abstract_basis_type
Teuchos::Array< int > inv_perm
Bi-directional iterator for traversing a sparse array.