287 typename MatrixType::ordinal_type pce_size,
288 const CijkType& cijk) {
289 typedef typename MatrixType::ordinal_type ordinal_type;
290 typedef typename MatrixType::StaticCrsGraphType matrix_graph_type;
291 typedef typename MatrixType::values_type matrix_values_type;
293 std::vector< std::vector<ordinal_type> > graph(nrow);
294 for (ordinal_type i=0; i<nrow; ++i)
295 graph[i] = std::vector<ordinal_type>(1, i);
296 ordinal_type graph_length = nrow;
298 matrix_graph_type matrix_graph =
299 Kokkos::create_staticcrsgraph<matrix_graph_type>(
"graph", graph);
300 matrix_values_type matrix_values =
301 Kokkos::make_view<matrix_values_type>(
"values", cijk, graph_length, pce_size);
303 MatrixType matrix(
"matrix", nrow, matrix_values, matrix_graph);
543 const typename PCEType::ordinal_type
stoch_dim,
544 const typename PCEType::ordinal_type
poly_ord,
545 KokkosSparse::DeviceConfig dev_config,
546 Multiply multiply_op,
547 Teuchos::FancyOStream& out)
549 typedef typename PCEType::ordinal_type ordinal_type;
550 typedef typename PCEType::value_type scalar_type;
552 typedef typename PCEType::cijk_type cijk_type;
554 typedef Kokkos::LayoutLeft Layout;
555 typedef Kokkos::View< PCEType*, Layout, execution_space > block_vector_type;
556 typedef KokkosSparse::CrsMatrix< PCEType, ordinal_type, execution_space > block_matrix_type;
557 typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
558 typedef typename block_matrix_type::values_type matrix_values_type;
562 const ordinal_type stoch_length = cijk.dimension();
565 const ordinal_type stoch_length_aligned = stoch_length;
568 TEUCHOS_TEST_FOR_EXCEPTION(
569 storage_type::is_static && storage_type::static_size != stoch_length,
571 "Static storage size must equal pce size");
574 const ordinal_type fem_length = nGrid * nGrid * nGrid;
575 std::vector< std::vector<ordinal_type> > fem_graph;
582 block_vector_type x =
583 Kokkos::make_view<block_vector_type>(
"x", cijk, fem_length, stoch_length_aligned);
584 block_vector_type y =
585 Kokkos::make_view<block_vector_type>(
"y", cijk, fem_length, stoch_length_aligned);
591 typename block_vector_type::HostMirror::array_type hax = hx ;
592 typename block_vector_type::HostMirror::array_type hay = hy ;
594 for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
595 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
596 hax(iRowStoch, iRowFEM) =
597 generate_vector_coefficient<scalar_type>(
598 fem_length, stoch_length, iRowFEM, iRowStoch );
599 hay(iRowStoch, iRowFEM) = 0.0;
609 matrix_graph_type matrix_graph =
610 Kokkos::create_staticcrsgraph<matrix_graph_type>(
611 std::string(
"test crs graph"), fem_graph);
612 matrix_values_type matrix_values =
613 Kokkos::make_view<matrix_values_type>(
614 Kokkos::ViewAllocateWithoutInitializing(
"matrix"), cijk, fem_graph_length, stoch_length_aligned);
615 block_matrix_type matrix(
616 "block_matrix", fem_length, matrix_values, matrix_graph);
617 matrix.dev_config = dev_config;
619 typename matrix_values_type::HostMirror hM =
622 typename matrix_values_type::HostMirror::array_type haM = hM ;
624 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
625 const ordinal_type row_size = fem_graph[iRowFEM].size();
626 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
627 ++iRowEntryFEM, ++iEntryFEM) {
628 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
630 for (ordinal_type k=0; k<stoch_length; ++k) {
632 generate_matrix_coefficient<scalar_type>(
633 fem_length, stoch_length, iRowFEM, iColFEM, k);
643 multiply_op( matrix, x, y );
648 typedef typename block_vector_type::array_type array_type;
649 array_type ay_expected =
650 array_type(
"ay_expected", stoch_length_aligned, fem_length);
651 typename array_type::HostMirror hay_expected =
653 typename cijk_type::HostMirror host_cijk =
656 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
657 const ordinal_type row_size = fem_graph[iRowFEM].size();
658 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
659 ++iRowEntryFEM, ++iEntryFEM) {
660 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
661 for (ordinal_type i=0; i<stoch_length; ++i) {
662 const ordinal_type num_entry = host_cijk.num_entry(i);
663 const ordinal_type entry_beg = host_cijk.entry_begin(i);
664 const ordinal_type entry_end = entry_beg + num_entry;
666 for (ordinal_type entry = entry_beg; entry < entry_end; ++entry) {
667 const ordinal_type
j = host_cijk.coord(entry,0);
668 const ordinal_type k = host_cijk.coord(entry,1);
669 const scalar_type a_j =
670 generate_matrix_coefficient<scalar_type>(
671 fem_length, stoch_length, iRowFEM, iColFEM,
j);
672 const scalar_type a_k =
673 generate_matrix_coefficient<scalar_type>(
674 fem_length, stoch_length, iRowFEM, iColFEM, k);
675 const scalar_type x_j =
676 generate_vector_coefficient<scalar_type>(
677 fem_length, stoch_length, iColFEM,
j);
678 const scalar_type x_k =
679 generate_vector_coefficient<scalar_type>(
680 fem_length, stoch_length, iColFEM, k);
681 tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
683 hay_expected(i, iRowFEM) += tmp;
692 typename block_vector_type::array_type ay = y;
737 Kokkos_CrsMatrix_PCE, MeanMultiplyRank1,
Scalar )
739 typedef typename Scalar::ordinal_type
Ordinal;
744 KokkosSparse::DeviceConfig dev_config;
746 typedef typename Scalar::ordinal_type ordinal_type;
747 typedef typename Scalar::value_type scalar_type;
749 typedef typename Scalar::cijk_type cijk_type;
751 typedef Kokkos::LayoutLeft Layout;
752 typedef Kokkos::View< Scalar*, Layout, execution_space > block_vector_type;
753 typedef KokkosSparse::CrsMatrix< Scalar, ordinal_type, execution_space > block_matrix_type;
754 typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
755 typedef typename block_matrix_type::values_type matrix_values_type;
759 cijk_type mean_cijk =
760 Stokhos::create_mean_based_product_tensor<execution_space, ordinal_type, scalar_type>();
761 const ordinal_type stoch_length = cijk.dimension();
762 const ordinal_type align = 8;
763 const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
766 const ordinal_type fem_length = nGrid * nGrid * nGrid;
767 std::vector< std::vector<ordinal_type> > fem_graph;
770 block_vector_type x =
771 Kokkos::make_view<block_vector_type>(
"x", cijk, fem_length, stoch_length_aligned);
772 block_vector_type y =
773 Kokkos::make_view<block_vector_type>(
"y", cijk, fem_length, stoch_length_aligned);
775 block_vector_type y_expected =
776 Kokkos::make_view<block_vector_type>(
"y", cijk, fem_length, stoch_length_aligned);
780 typename block_vector_type::HostMirror hy_expected =
784 typename block_vector_type::HostMirror::array_type hax = hx ;
785 typename block_vector_type::HostMirror::array_type hay = hy ;
786 typename block_vector_type::HostMirror::array_type hay_expected =
789 for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
790 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
791 hax(iRowStoch,iRowFEM) =
792 generate_vector_coefficient<scalar_type>(
793 fem_length, stoch_length, iRowFEM, iRowStoch );
794 hay(iRowStoch,iRowFEM) = 0.0;
795 hay_expected(iRowStoch,iRowFEM) = 0.0;
804 matrix_graph_type matrix_graph =
805 Kokkos::create_staticcrsgraph<matrix_graph_type>(
806 std::string(
"test crs graph"), fem_graph);
807 matrix_values_type matrix_values =
808 Kokkos::make_view<matrix_values_type>(
809 Kokkos::ViewAllocateWithoutInitializing(
"matrix"), mean_cijk, fem_graph_length, ordinal_type(1));
810 block_matrix_type matrix(
811 "block_matrix", fem_length, matrix_values, matrix_graph);
812 matrix.dev_config = dev_config;
814 typename matrix_values_type::HostMirror hM =
816 typename matrix_values_type::HostMirror::array_type haM = hM ;
818 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
819 const ordinal_type row_size = fem_graph[iRowFEM].size();
820 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
821 ++iRowEntryFEM, ++iEntryFEM) {
822 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
825 generate_matrix_coefficient<scalar_type>(
826 fem_length, 1, iRowFEM, iColFEM, 0);
827 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
828 hay_expected(iRowStoch,iRowFEM) +=
829 haM(iEntryFEM, 0) * hax(iRowStoch,iColFEM);
883 success =
compareRank1(y, y_expected, rel_tol, abs_tol, out);
888 Kokkos_CrsMatrix_PCE, MeanMultiplyRank2,
Scalar )
890 typedef typename Scalar::ordinal_type ordinal_type;
891 typedef typename Scalar::value_type scalar_type;
893 typedef typename Scalar::cijk_type cijk_type;
895 typedef Kokkos::LayoutLeft Layout;
896 typedef Kokkos::View< Scalar**, Layout, execution_space > block_vector_type;
897 typedef KokkosSparse::CrsMatrix< Scalar, ordinal_type, execution_space > block_matrix_type;
898 typedef typename block_matrix_type::StaticCrsGraphType matrix_graph_type;
899 typedef typename block_matrix_type::values_type matrix_values_type;
902 const ordinal_type nGrid = 5;
905 KokkosSparse::DeviceConfig dev_config;
909 cijk_type mean_cijk =
910 Stokhos::create_mean_based_product_tensor<execution_space, ordinal_type, scalar_type>();
911 const ordinal_type stoch_length = cijk.dimension();
912 const ordinal_type align = 8;
913 const ordinal_type stoch_length_aligned = (stoch_length+align-1) & ~(align-1);
914 const ordinal_type num_cols = 2;
916 const ordinal_type fem_length = nGrid * nGrid * nGrid;
917 std::vector< std::vector<ordinal_type> > fem_graph;
920 block_vector_type x =
921 Kokkos::make_view<block_vector_type>(
"x", cijk, fem_length, num_cols, stoch_length_aligned);
922 block_vector_type y =
923 Kokkos::make_view<block_vector_type>(
"y", cijk, fem_length, num_cols, stoch_length_aligned);
925 block_vector_type y_expected =
926 Kokkos::make_view<block_vector_type>(
"y_expected", cijk, fem_length, num_cols, stoch_length_aligned);
930 typename block_vector_type::HostMirror hy_expected =
933 for (ordinal_type i=0; i<num_cols; ++i){
934 for (ordinal_type iRowFEM=0; iRowFEM<fem_length; ++iRowFEM) {
935 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
936 hx(iRowFEM,i).fastAccessCoeff(iRowStoch) =
937 generate_vector_coefficient<scalar_type>(
938 fem_length, stoch_length, iRowFEM, iRowStoch );
939 hy(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
940 hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) = 0.0;
950 matrix_graph_type matrix_graph =
951 Kokkos::create_staticcrsgraph<matrix_graph_type>(
952 std::string(
"test crs graph"), fem_graph);
953 matrix_values_type matrix_values =
954 Kokkos::make_view<matrix_values_type>(
955 Kokkos::ViewAllocateWithoutInitializing(
"matrix"), mean_cijk, fem_graph_length, ordinal_type(1));
956 block_matrix_type matrix(
957 "block_matrix", fem_length, matrix_values, matrix_graph);
958 matrix.dev_config = dev_config;
960 typename matrix_values_type::HostMirror hM =
963 typename matrix_values_type::HostMirror::array_type haM = hM ;
965 for (ordinal_type iRowFEM=0, iEntryFEM=0; iRowFEM<fem_length; ++iRowFEM) {
966 const ordinal_type row_size = fem_graph[iRowFEM].size();
967 for (ordinal_type iRowEntryFEM=0; iRowEntryFEM<row_size;
968 ++iRowEntryFEM, ++iEntryFEM) {
969 const ordinal_type iColFEM = fem_graph[iRowFEM][iRowEntryFEM];
972 generate_matrix_coefficient<scalar_type>(
973 fem_length, 1, iRowFEM, iColFEM, 0);
974 for (ordinal_type i=0; i<num_cols; ++i){
975 for (ordinal_type iRowStoch=0; iRowStoch<stoch_length; ++iRowStoch) {
976 hy_expected(iRowFEM,i).fastAccessCoeff(iRowStoch) +=
977 haM(iEntryFEM, 0) * hx(iColFEM,i).fastAccessCoeff(iRowStoch);
1035 success =
compareRank2(y, y_expected, rel_tol, abs_tol, out);