52 #include "Tpetra_Core.hpp"
53 #include "Tpetra_Map.hpp"
54 #include "Tpetra_MultiVector.hpp"
55 #include "Tpetra_Vector.hpp"
56 #include "Tpetra_CrsGraph.hpp"
57 #include "Tpetra_CrsMatrix.hpp"
61 #ifdef HAVE_STOKHOS_BELOS
63 #include "BelosLinearProblem.hpp"
64 #include "BelosPseudoBlockGmresSolMgr.hpp"
65 #include "BelosPseudoBlockCGSolMgr.hpp"
69 #ifdef HAVE_STOKHOS_IFPACK2
71 #include "Ifpack2_Factory.hpp"
75 #ifdef HAVE_STOKHOS_MUELU
77 #include "MueLu_CreateTpetraPreconditioner.hpp"
81 #ifdef HAVE_STOKHOS_AMESOS2
83 #include "Amesos2_Factory.hpp"
95 #define USE_SCALAR_MEAN_BASED_PREC 1
97 template <
typename scalar,
typename ordinal>
100 const ordinal nStoch,
101 const ordinal iColFEM,
102 const ordinal iStoch )
104 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
105 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
106 return X_fem + X_stoch;
110 template <
typename scalar,
typename ordinal>
114 const ordinal nStoch,
115 const ordinal iColFEM,
117 const ordinal iStoch)
119 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
120 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
121 const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
122 return X_fem + X_stoch + X_col;
126 template <
typename scalar,
typename ordinal>
129 const ordinal nStoch,
130 const ordinal iRowFEM,
131 const ordinal iColFEM,
132 const ordinal iStoch )
134 const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
135 ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
137 const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
139 return A_fem + A_stoch;
143 template <
typename kokkos_cijk_type,
typename ordinal_type>
159 Array< RCP<const one_d_basis> > bases(stoch_dim);
161 bases[i] =
rcp(
new legendre_basis(poly_ord,
true));
162 RCP<const product_basis> basis =
rcp(
new product_basis(bases));
165 RCP<Cijk>
cijk = basis->computeTripleProductTensor();
168 kokkos_cijk_type kokkos_cijk =
169 Stokhos::create_product_tensor<execution_space>(*basis, *
cijk);
196 typedef typename Scalar::cijk_type Cijk;
199 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
200 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
203 if ( !Kokkos::is_initialized() )
204 Kokkos::initialize();
212 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
216 RCP<const Tpetra_Map> map =
217 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
219 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
220 const size_t num_my_row = myGIDs.size();
223 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
224 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
225 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
226 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
227 Scalar val1(cijk), val2(cijk);
228 for (
size_t i=0; i<num_my_row; ++i) {
231 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
232 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
243 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
244 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
250 ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
252 BaseScalar tol = 1.0e-14;
253 for (
size_t i=0; i<num_my_row; ++i) {
256 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
257 nrow, pce_size, row,
j);
258 val.fastAccessCoeff(
j) = alpha.coeff(0)*v + 0.12345*beta.coeff(0)*v;
283 typedef typename Scalar::cijk_type Cijk;
286 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
287 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
288 typedef typename Tpetra_Vector::dot_type dot_type;
291 if ( !Kokkos::is_initialized() )
292 Kokkos::initialize();
300 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
304 RCP<const Tpetra_Map> map =
305 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
307 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
308 const size_t num_my_row = myGIDs.size();
311 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
312 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
313 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
314 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
315 Scalar val1(cijk), val2(cijk);
316 for (
size_t i=0; i<num_my_row; ++i) {
319 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
320 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
329 dot_type
dot = x1->dot(*x2);
334 dot_type local_val(0);
335 for (
size_t i=0; i<num_my_row; ++i) {
338 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
339 nrow, pce_size, row,
j);
340 local_val += 0.12345 * v * v;
347 Teuchos::outArg(val));
349 out <<
"dot = " << dot <<
" expected = " << val << std::endl;
351 BaseScalar tol = 1.0e-14;
372 typedef typename Scalar::cijk_type Cijk;
375 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
376 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
379 if ( !Kokkos::is_initialized() )
380 Kokkos::initialize();
388 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
392 RCP<const Tpetra_Map> map =
393 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
395 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
396 const size_t num_my_row = myGIDs.size();
400 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
401 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
402 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
403 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
404 Scalar val1(cijk), val2(cijk);
405 for (
size_t i=0; i<num_my_row; ++i) {
407 for (
size_t j=0;
j<ncol; ++
j) {
410 generate_multi_vector_coefficient<BaseScalar,size_t>(
411 nrow, ncol, pce_size, row,
j, k);
412 val1.fastAccessCoeff(k) = v;
413 val2.fastAccessCoeff(k) = 0.12345 * v;
415 x1_view[
j][i] = val1;
416 x2_view[
j][i] = val2;
425 RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
426 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
432 ArrayRCP< ArrayRCP<Scalar> > y_view = y->get2dViewNonConst();
434 BaseScalar tol = 1.0e-14;
435 for (
size_t i=0; i<num_my_row; ++i) {
437 for (
size_t j=0;
j<ncol; ++
j) {
439 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
440 nrow, ncol, pce_size, row,
j, k);
441 val.fastAccessCoeff(k) = alpha.coeff(0)*v + 0.12345*beta.coeff(0)*v;
446 val.fastAccessCoeff(k), tol );
468 typedef typename Scalar::cijk_type Cijk;
471 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
472 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
473 typedef typename Tpetra_MultiVector::dot_type dot_type;
476 if ( !Kokkos::is_initialized() )
477 Kokkos::initialize();
485 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
489 RCP<const Tpetra_Map> map =
490 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
492 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
493 const size_t num_my_row = myGIDs.size();
497 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
498 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
499 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
500 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
501 Scalar val1(cijk), val2(cijk);
502 for (
size_t i=0; i<num_my_row; ++i) {
504 for (
size_t j=0;
j<ncol; ++
j) {
507 generate_multi_vector_coefficient<BaseScalar,size_t>(
508 nrow, ncol, pce_size, row,
j, k);
509 val1.fastAccessCoeff(k) = v;
510 val2.fastAccessCoeff(k) = 0.12345 * v;
512 x1_view[
j][i] = val1;
513 x2_view[
j][i] = val2;
520 Array<dot_type> dots(ncol);
521 x1->dot(*x2, dots());
526 Array<dot_type> local_vals(ncol, dot_type(0));
527 for (
size_t i=0; i<num_my_row; ++i) {
529 for (
size_t j=0;
j<ncol; ++
j) {
531 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
532 nrow, ncol, pce_size, row,
j, k);
533 local_vals[
j] += 0.12345 * v * v;
539 Array<dot_type> vals(ncol, dot_type(0));
541 local_vals.getRawPtr(), vals.getRawPtr());
543 BaseScalar tol = 1.0e-14;
544 for (
size_t j=0;
j<ncol; ++
j) {
545 out <<
"dots(" <<
j <<
") = " << dots[
j]
546 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
568 typedef typename Scalar::cijk_type Cijk;
571 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
572 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
573 typedef typename Tpetra_MultiVector::dot_type dot_type;
576 if ( !Kokkos::is_initialized() )
577 Kokkos::initialize();
585 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
589 RCP<const Tpetra_Map> map =
590 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
592 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
593 const size_t num_my_row = myGIDs.size();
597 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
598 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
599 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
600 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
601 Scalar val1(cijk), val2(cijk);
602 for (
size_t i=0; i<num_my_row; ++i) {
604 for (
size_t j=0;
j<ncol; ++
j) {
607 generate_multi_vector_coefficient<BaseScalar,size_t>(
608 nrow, ncol, pce_size, row,
j, k);
609 val1.fastAccessCoeff(k) = v;
610 val2.fastAccessCoeff(k) = 0.12345 * v;
612 x1_view[
j][i] = val1;
613 x2_view[
j][i] = val2;
622 cols[0] = 4; cols[1] = 2;
623 RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
624 RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
627 Array<dot_type> dots(ncol_sub);
628 x1_sub->dot(*x2_sub, dots());
633 Array<dot_type> local_vals(ncol_sub, dot_type(0));
634 for (
size_t i=0; i<num_my_row; ++i) {
636 for (
size_t j=0;
j<ncol_sub; ++
j) {
638 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
639 nrow, ncol, pce_size, row, cols[
j], k);
640 local_vals[
j] += 0.12345 * v * v;
646 Array<dot_type> vals(ncol_sub, dot_type(0));
648 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
651 BaseScalar tol = 1.0e-14;
652 for (
size_t j=0;
j<ncol_sub; ++
j) {
653 out <<
"dots(" <<
j <<
") = " << dots[
j]
654 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
676 typedef typename Scalar::cijk_type Cijk;
679 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
680 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
681 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
682 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
685 if ( !Kokkos::is_initialized() )
686 Kokkos::initialize();
695 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
696 RCP<const Tpetra_Map> map =
697 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
699 RCP<Tpetra_CrsGraph> graph =
700 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
701 Array<GlobalOrdinal> columnIndices(2);
702 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
703 const size_t num_my_row = myGIDs.size();
704 for (
size_t i=0; i<num_my_row; ++i) {
706 columnIndices[0] = row;
709 columnIndices[1] = row+1;
712 graph->insertGlobalIndices(row, columnIndices(0,ncol));
714 graph->fillComplete();
715 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
718 Array<Scalar> vals(2);
720 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
722 const size_t num_col = row == nrow - 1 ? 1 : 2;
723 for (
size_t local_col=0; local_col<num_col; ++local_col) {
725 columnIndices[local_col] = col;
727 val.fastAccessCoeff(k) =
728 generate_matrix_coefficient<BaseScalar,size_t>(
729 nrow, pce_size, row, col, k);
730 vals[local_col] =
val;
732 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
734 matrix->fillComplete();
737 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
738 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
739 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
742 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
743 nrow, pce_size, row,
j);
744 x_view[local_row] =
val;
755 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
756 matrix->apply(*x, *y);
762 ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
763 BaseScalar tol = 1.0e-14;
764 typename Cijk::HostMirror host_cijk =
767 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
769 const size_t num_col = row == nrow - 1 ? 1 : 2;
771 for (
size_t local_col=0; local_col<num_col; ++local_col) {
775 const LocalOrdinal entry_beg = host_cijk.entry_begin(i);
778 for (
LocalOrdinal entry = entry_beg; entry < entry_end; ++entry) {
781 const BaseScalar a_j =
782 generate_matrix_coefficient<BaseScalar,size_t>(
783 nrow, pce_size, row, col,
j);
784 const BaseScalar a_k =
785 generate_matrix_coefficient<BaseScalar,size_t>(
786 nrow, pce_size, row, col, k);
787 const BaseScalar x_j =
788 generate_vector_coefficient<BaseScalar,size_t>(
789 nrow, pce_size, col,
j);
790 const BaseScalar x_k =
791 generate_vector_coefficient<BaseScalar,size_t>(
792 nrow, pce_size, col, k);
793 tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
795 val.fastAccessCoeff(i) += tmp;
801 val.fastAccessCoeff(i), tol );
822 typedef typename Scalar::cijk_type Cijk;
825 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
826 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
827 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
828 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
831 if ( !Kokkos::is_initialized() )
832 Kokkos::initialize();
841 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
842 RCP<const Tpetra_Map> map =
843 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
845 RCP<Tpetra_CrsGraph> graph =
846 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
847 Array<GlobalOrdinal> columnIndices(2);
848 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
849 const size_t num_my_row = myGIDs.size();
850 for (
size_t i=0; i<num_my_row; ++i) {
852 columnIndices[0] = row;
855 columnIndices[1] = row+1;
858 graph->insertGlobalIndices(row, columnIndices(0,ncol));
860 graph->fillComplete();
861 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
864 Array<Scalar> vals(2);
866 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
868 const size_t num_col = row == nrow - 1 ? 1 : 2;
869 for (
size_t local_col=0; local_col<num_col; ++local_col) {
871 columnIndices[local_col] = col;
873 val.fastAccessCoeff(k) =
874 generate_matrix_coefficient<BaseScalar,size_t>(
875 nrow, pce_size, row, col, k);
876 vals[local_col] =
val;
878 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
880 matrix->fillComplete();
884 RCP<Tpetra_MultiVector> x = Tpetra::createMultiVector<Scalar>(map, ncol);
885 ArrayRCP< ArrayRCP<Scalar> > x_view = x->get2dViewNonConst();
886 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
888 for (
size_t col=0; col<ncol; ++col) {
891 generate_multi_vector_coefficient<BaseScalar,size_t>(
892 nrow, ncol, pce_size, row, col, k);
893 val.fastAccessCoeff(k) = v;
895 x_view[col][local_row] =
val;
907 RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
908 matrix->apply(*x, *y);
914 ArrayRCP< ArrayRCP<Scalar> > y_view = y->get2dViewNonConst();
915 BaseScalar tol = 1.0e-14;
916 typename Cijk::HostMirror host_cijk =
919 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
921 for (
size_t xcol=0; xcol<ncol; ++xcol) {
922 const size_t num_col = row == nrow - 1 ? 1 : 2;
924 for (
size_t local_col=0; local_col<num_col; ++local_col) {
928 const LocalOrdinal entry_beg = host_cijk.entry_begin(i);
931 for (
LocalOrdinal entry = entry_beg; entry < entry_end; ++entry) {
934 const BaseScalar a_j =
935 generate_matrix_coefficient<BaseScalar,size_t>(
936 nrow, pce_size, row, col,
j);
937 const BaseScalar a_k =
938 generate_matrix_coefficient<BaseScalar,size_t>(
939 nrow, pce_size, row, col, k);
940 const BaseScalar x_j =
941 generate_multi_vector_coefficient<BaseScalar,size_t>(
942 nrow, ncol, pce_size, col, xcol,
j);
943 const BaseScalar x_k =
944 generate_multi_vector_coefficient<BaseScalar,size_t>(
945 nrow, ncol, pce_size, col, xcol, k);
946 tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
948 val.fastAccessCoeff(i) += tmp;
954 val.fastAccessCoeff(i), tol );
976 typedef typename Scalar::cijk_type Cijk;
979 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
980 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
981 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
982 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
984 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
985 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
988 if ( !Kokkos::is_initialized() )
989 Kokkos::initialize();
998 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
999 RCP<const Tpetra_Map> map =
1000 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1002 RCP<Tpetra_CrsGraph> graph =
1003 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1004 Array<GlobalOrdinal> columnIndices(2);
1005 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1006 const size_t num_my_row = myGIDs.size();
1007 for (
size_t i=0; i<num_my_row; ++i) {
1009 columnIndices[0] = row;
1011 if (row != nrow-1) {
1012 columnIndices[1] = row+1;
1015 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1017 graph->fillComplete();
1018 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1021 Array<Scalar> vals(2);
1023 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1025 const size_t num_col = row == nrow - 1 ? 1 : 2;
1026 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1028 columnIndices[local_col] = col;
1030 val.fastAccessCoeff(k) =
1031 generate_matrix_coefficient<BaseScalar,size_t>(
1032 nrow, pce_size, row, col, k);
1033 vals[local_col] =
val;
1035 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1037 matrix->fillComplete();
1040 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1041 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1042 for (
size_t i=0; i<num_my_row; ++i) {
1045 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
1046 nrow, pce_size, row,
j);
1051 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
1052 matrix->apply(*x, *y);
1069 RCP<const Tpetra_Map> flat_x_map, flat_y_map;
1070 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1072 Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_y_map,
1073 cijk_graph, pce_size);
1074 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1078 RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
1079 RCP<Flat_Tpetra_Vector> flat_x =
1081 RCP<Flat_Tpetra_Vector> flat_y =
1083 flat_matrix->apply(*flat_x, *flat_y);
1103 BaseScalar tol = 1.0e-14;
1104 ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
1105 ArrayRCP<Scalar> y2_view = y2->get1dViewNonConst();
1106 for (
size_t i=0; i<num_my_row; ++i) {
1133 typedef typename Scalar::cijk_type Cijk;
1136 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1137 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1138 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1139 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1142 if ( !Kokkos::is_initialized() )
1143 Kokkos::initialize();
1152 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1153 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1154 RCP<const Tpetra_Map> map =
1155 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1157 RCP<Tpetra_CrsGraph> graph =
1158 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1159 Array<GlobalOrdinal> columnIndices(3);
1160 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1161 const size_t num_my_row = myGIDs.size();
1162 for (
size_t i=0; i<num_my_row; ++i) {
1164 if (row == 0 || row == nrow-1) {
1165 columnIndices[0] = row;
1166 graph->insertGlobalIndices(row, columnIndices(0,1));
1169 columnIndices[0] = row-1;
1170 columnIndices[1] = row;
1171 columnIndices[2] = row+1;
1172 graph->insertGlobalIndices(row, columnIndices(0,3));
1175 graph->fillComplete();
1176 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1179 Array<Scalar> vals(3);
1182 a_val.fastAccessCoeff(
j) =
1183 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1185 for (
size_t i=0; i<num_my_row; ++i) {
1187 if (row == 0 || row == nrow-1) {
1188 columnIndices[0] = row;
1190 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1193 columnIndices[0] = row-1;
1194 columnIndices[1] = row;
1195 columnIndices[2] = row+1;
1196 vals[0] =
Scalar(1.0) * a_val;
1197 vals[1] =
Scalar(-2.0) * a_val;
1198 vals[2] =
Scalar(1.0) * a_val;
1199 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1202 matrix->fillComplete();
1208 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1209 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1212 b_val.fastAccessCoeff(
j) =
1213 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1215 for (
size_t i=0; i<num_my_row; ++i) {
1217 if (row == 0 || row == nrow-1)
1220 b_view[i] = b_val * (h*h);
1227 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1228 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1229 typedef typename BST::mag_type base_mag_type;
1230 typedef typename Tpetra_Vector::mag_type mag_type;
1231 base_mag_type btol = 1e-9;
1232 mag_type tol = btol;
1235 out.getOStream().get());
1242 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1243 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1244 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1245 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1247 Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1248 cijk_graph, pce_size);
1249 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1251 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1252 RCP<Flat_Tpetra_Vector> flat_x =
1254 RCP<Flat_Tpetra_Vector> flat_b =
1257 tol, max_its, out.getOStream().get());
1258 TEST_EQUALITY_CONST( solved_flat,
true );
1261 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1262 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1263 for (
size_t i=0; i<num_my_row; ++i) {
1268 Scalar v = x_view[i];
1269 Scalar v2 = x2_view[i];
1271 if (
j < v.size() &&
BST::abs(v.coeff(
j)) < btol)
1272 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1273 if (
j < v2.size() &&
BST::abs(v2.coeff(
j)) < btol)
1274 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1286 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION)
1303 using Teuchos::getParametersFromXmlFile;
1307 typedef typename Scalar::cijk_type Cijk;
1310 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1311 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1312 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1313 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1316 if ( !Kokkos::is_initialized() )
1317 Kokkos::initialize();
1326 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1327 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1328 RCP<const Tpetra_Map> map =
1329 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1331 RCP<Tpetra_CrsGraph> graph =
1332 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1333 Array<GlobalOrdinal> columnIndices(3);
1334 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1335 const size_t num_my_row = myGIDs.size();
1336 for (
size_t i=0; i<num_my_row; ++i) {
1338 if (row == 0 || row == nrow-1) {
1339 columnIndices[0] = row;
1340 graph->insertGlobalIndices(row, columnIndices(0,1));
1343 columnIndices[0] = row-1;
1344 columnIndices[1] = row;
1345 columnIndices[2] = row+1;
1346 graph->insertGlobalIndices(row, columnIndices(0,3));
1349 graph->fillComplete();
1350 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1353 Array<Scalar> vals(3);
1356 a_val.fastAccessCoeff(
j) =
1357 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1359 for (
size_t i=0; i<num_my_row; ++i) {
1361 if (row == 0 || row == nrow-1) {
1362 columnIndices[0] = row;
1364 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1367 columnIndices[0] = row-1;
1368 columnIndices[1] = row;
1369 columnIndices[2] = row+1;
1370 vals[0] =
Scalar(1.0) * a_val;
1371 vals[1] =
Scalar(-2.0) * a_val;
1372 vals[2] =
Scalar(1.0) * a_val;
1373 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1376 matrix->fillComplete();
1379 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1380 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1383 b_val.fastAccessCoeff(
j) =
1384 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1386 for (
size_t i=0; i<num_my_row; ++i) {
1388 if (row == 0 || row == nrow-1)
1391 b_view[i] = b_val * (h*h);
1395 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1396 RCP<ParameterList> muelu_params =
1397 getParametersFromXmlFile(
"muelu_cheby.xml");
1398 #if USE_SCALAR_MEAN_BASED_PREC
1399 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_OP;
1400 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1401 RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1403 RCP<Scalar_OP> mean_matrix_op = mean_matrix;
1404 RCP<Scalar_OP> M_s =
1405 MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1409 Stokhos::create_mean_based_product_tensor<typename Storage::execution_space,typename Storage::ordinal_type,BaseScalar>();
1412 RCP<OP> mean_matrix_op = mean_matrix;
1414 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1419 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1420 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1421 typedef typename BST::mag_type base_mag_type;
1422 typedef typename Tpetra_Vector::mag_type mag_type;
1423 base_mag_type btol = 1e-9;
1424 mag_type tol = btol;
1427 out.getOStream().get());
1434 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1435 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1436 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1437 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1439 Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1440 cijk_graph, pce_size);
1441 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1443 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1444 RCP<Flat_Tpetra_Vector> flat_x =
1446 RCP<Flat_Tpetra_Vector> flat_b =
1454 tol, max_its, out.getOStream().get());
1458 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1459 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1460 for (
size_t i=0; i<num_my_row; ++i) {
1465 Scalar v = x_view[i];
1466 Scalar v2 = x2_view[i];
1468 if (
j < v.size() &&
BST::abs(v.coeff(
j)) < btol)
1469 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1470 if (
j < v2.size() &&
BST::abs(v2.coeff(
j)) < btol)
1471 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1490 #if defined(HAVE_STOKHOS_BELOS)
1507 typedef typename Scalar::cijk_type Cijk;
1510 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1511 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1512 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1513 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1516 if ( !Kokkos::is_initialized() )
1517 Kokkos::initialize();
1526 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1527 RCP<const Tpetra_Map> map =
1528 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1530 RCP<Tpetra_CrsGraph> graph =
1531 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1532 Array<GlobalOrdinal> columnIndices(2);
1533 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1534 const size_t num_my_row = myGIDs.size();
1535 for (
size_t i=0; i<num_my_row; ++i) {
1537 columnIndices[0] = row;
1539 if (row != nrow-1) {
1540 columnIndices[1] = row+1;
1543 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1545 graph->fillComplete();
1546 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1549 Array<Scalar> vals(2);
1551 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1553 const size_t num_col = row == nrow - 1 ? 1 : 2;
1554 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1556 columnIndices[local_col] = col;
1558 val.fastAccessCoeff(k) =
1559 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1560 vals[local_col] =
val;
1562 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1564 matrix->fillComplete();
1567 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1568 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1569 for (
size_t i=0; i<num_my_row; ++i) {
1575 typedef BaseScalar BelosScalar;
1576 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1577 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1578 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1579 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1580 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
1581 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
1582 typename ST::magnitudeType tol = 1e-9;
1583 belosParams->set(
"Flexible Gmres",
false);
1584 belosParams->set(
"Num Blocks", 100);
1585 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1586 belosParams->set(
"Maximum Iterations", 100);
1587 belosParams->set(
"Verbosity", 33);
1588 belosParams->set(
"Output Style", 1);
1589 belosParams->set(
"Output Frequency", 1);
1590 belosParams->set(
"Output Stream", out.getOStream());
1591 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1592 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1593 problem->setProblem();
1594 Belos::ReturnType ret = solver->solve();
1601 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1602 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1603 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1604 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1606 Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1607 cijk_graph, pce_size);
1608 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1610 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1611 RCP<Flat_Tpetra_Vector> flat_x =
1613 RCP<Flat_Tpetra_Vector> flat_b =
1615 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
1616 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1617 typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1618 RCP< FBLinProb > flat_problem =
1619 rcp(
new FBLinProb(flat_matrix, flat_x, flat_b));
1620 RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1621 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1623 flat_problem->setProblem();
1624 Belos::ReturnType flat_ret = flat_solver->solve();
1627 typename ST::magnitudeType btol = 100*tol;
1628 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1629 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1630 for (
size_t i=0; i<num_my_row; ++i) {
1635 Scalar v = x_view[i];
1636 Scalar v2 = x2_view[i];
1638 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
1639 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1640 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
1641 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1662 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2)
1680 typedef typename Scalar::cijk_type Cijk;
1683 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1684 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1685 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1686 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1689 if ( !Kokkos::is_initialized() )
1690 Kokkos::initialize();
1699 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1700 RCP<const Tpetra_Map> map =
1701 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1703 RCP<Tpetra_CrsGraph> graph =
1704 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1705 Array<GlobalOrdinal> columnIndices(2);
1706 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1707 const size_t num_my_row = myGIDs.size();
1708 for (
size_t i=0; i<num_my_row; ++i) {
1710 columnIndices[0] = row;
1712 if (row != nrow-1) {
1713 columnIndices[1] = row+1;
1716 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1718 graph->fillComplete();
1719 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1722 Array<Scalar> vals(2);
1724 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1726 const size_t num_col = row == nrow - 1 ? 1 : 2;
1727 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1729 columnIndices[local_col] = col;
1731 val.fastAccessCoeff(k) =
1732 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1733 vals[local_col] =
val;
1735 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1737 matrix->fillComplete();
1743 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1744 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1745 for (
size_t i=0; i<num_my_row; ++i) {
1750 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
1751 Ifpack2::Factory factory;
1752 RCP<Prec> M = factory.create<Tpetra_CrsMatrix>(
"RILUK", mean_matrix);
1758 typedef BaseScalar BelosScalar;
1759 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1760 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1761 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1762 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1763 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
1764 problem->setRightPrec(M);
1765 problem->setProblem();
1766 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
1767 typename ST::magnitudeType tol = 1e-9;
1768 belosParams->set(
"Flexible Gmres",
false);
1769 belosParams->set(
"Num Blocks", 100);
1770 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1771 belosParams->set(
"Maximum Iterations", 100);
1772 belosParams->set(
"Verbosity", 33);
1773 belosParams->set(
"Output Style", 1);
1774 belosParams->set(
"Output Frequency", 1);
1775 belosParams->set(
"Output Stream", out.getOStream());
1777 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1778 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1779 Belos::ReturnType ret = solver->solve();
1783 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1784 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1785 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1786 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1788 Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
1789 cijk_graph, pce_size);
1790 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1792 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1793 RCP<Flat_Tpetra_Vector> flat_x =
1795 RCP<Flat_Tpetra_Vector> flat_b =
1797 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
1798 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1799 typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1800 RCP< FBLinProb > flat_problem =
1801 rcp(
new FBLinProb(flat_matrix, flat_x, flat_b));
1802 RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1803 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1805 flat_problem->setProblem();
1806 Belos::ReturnType flat_ret = flat_solver->solve();
1809 typename ST::magnitudeType btol = 100*tol;
1810 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1811 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1812 for (
size_t i=0; i<num_my_row; ++i) {
1817 Scalar v = x_view[i];
1818 Scalar v2 = x2_view[i];
1820 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
1821 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1822 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
1823 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1842 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION)
1859 using Teuchos::getParametersFromXmlFile;
1863 typedef typename Scalar::cijk_type Cijk;
1866 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1867 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1868 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1869 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1872 if ( !Kokkos::is_initialized() )
1873 Kokkos::initialize();
1882 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1883 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1884 RCP<const Tpetra_Map> map =
1885 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1887 RCP<Tpetra_CrsGraph> graph =
1888 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1889 Array<GlobalOrdinal> columnIndices(3);
1890 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1891 const size_t num_my_row = myGIDs.size();
1892 for (
size_t i=0; i<num_my_row; ++i) {
1894 if (row == 0 || row == nrow-1) {
1895 columnIndices[0] = row;
1896 graph->insertGlobalIndices(row, columnIndices(0,1));
1899 columnIndices[0] = row-1;
1900 columnIndices[1] = row;
1901 columnIndices[2] = row+1;
1902 graph->insertGlobalIndices(row, columnIndices(0,3));
1905 graph->fillComplete();
1906 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1909 Array<Scalar> vals(3);
1912 a_val.fastAccessCoeff(
j) =
1913 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1915 for (
size_t i=0; i<num_my_row; ++i) {
1917 if (row == 0 || row == nrow-1) {
1918 columnIndices[0] = row;
1920 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1923 columnIndices[0] = row-1;
1924 columnIndices[1] = row;
1925 columnIndices[2] = row+1;
1926 vals[0] =
Scalar(1.0) * a_val;
1927 vals[1] =
Scalar(-2.0) * a_val;
1928 vals[2] =
Scalar(1.0) * a_val;
1929 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1932 matrix->fillComplete();
1935 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1936 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1939 b_val.fastAccessCoeff(
j) =
1940 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1942 for (
size_t i=0; i<num_my_row; ++i) {
1944 if (row == 0 || row == nrow-1)
1947 b_view[i] = b_val * (h*h);
1951 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1952 RCP<ParameterList> muelu_params =
1953 getParametersFromXmlFile(
"muelu_cheby.xml");
1954 #if USE_SCALAR_MEAN_BASED_PREC
1955 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_OP;
1956 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1957 RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1959 RCP<Scalar_OP> mean_matrix_op = mean_matrix;
1960 RCP<Scalar_OP> M_s =
1961 MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1965 Stokhos::create_mean_based_product_tensor<typename Storage::execution_space,typename Storage::ordinal_type,BaseScalar>();
1968 RCP<OP> mean_matrix_op = mean_matrix;
1970 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1976 typedef BaseScalar BelosScalar;
1977 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1978 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1979 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1980 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
1981 problem->setRightPrec(M);
1982 problem->setProblem();
1983 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
1984 typename ST::magnitudeType tol = 1e-9;
1985 belosParams->set(
"Num Blocks", 100);
1986 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1987 belosParams->set(
"Maximum Iterations", 100);
1988 belosParams->set(
"Verbosity", 33);
1989 belosParams->set(
"Output Style", 1);
1990 belosParams->set(
"Output Frequency", 1);
1991 belosParams->set(
"Output Stream", out.getOStream());
1993 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1994 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1997 Belos::ReturnType ret = solver->solve();
2004 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
2005 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2006 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2007 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2009 Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
2010 cijk_graph, pce_size);
2011 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2013 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2014 RCP<Flat_Tpetra_Vector> flat_x =
2016 RCP<Flat_Tpetra_Vector> flat_b =
2018 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
2019 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
2020 typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
2021 RCP< FBLinProb > flat_problem =
2022 rcp(
new FBLinProb(flat_matrix, flat_x, flat_b));
2023 RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
2024 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
2029 flat_problem->setProblem();
2030 Belos::ReturnType flat_ret = flat_solver->solve();
2033 typename ST::magnitudeType btol = 100*tol;
2034 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2035 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2036 for (
size_t i=0; i<num_my_row; ++i) {
2041 Scalar v = x_view[i];
2042 Scalar v2 = x2_view[i];
2044 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
2045 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2046 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
2047 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
2067 #if defined(HAVE_STOKHOS_AMESOS2)
2084 typedef typename Scalar::cijk_type Cijk;
2087 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2088 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2089 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
2090 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2091 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2094 if ( !Kokkos::is_initialized() )
2095 Kokkos::initialize();
2104 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
2105 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2106 RCP<const Tpetra_Map> map =
2107 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2109 RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map,
size_t(3));
2110 Array<GlobalOrdinal> columnIndices(3);
2111 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
2112 const size_t num_my_row = myGIDs.size();
2113 for (
size_t i=0; i<num_my_row; ++i) {
2115 if (row == 0 || row == nrow-1) {
2116 columnIndices[0] = row;
2117 graph->insertGlobalIndices(row, columnIndices(0,1));
2120 columnIndices[0] = row-1;
2121 columnIndices[1] = row;
2122 columnIndices[2] = row+1;
2123 graph->insertGlobalIndices(row, columnIndices(0,3));
2126 graph->fillComplete();
2127 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
2130 Array<Scalar> vals(3);
2133 a_val.fastAccessCoeff(
j) =
2134 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
2136 for (
size_t i=0; i<num_my_row; ++i) {
2138 if (row == 0 || row == nrow-1) {
2139 columnIndices[0] = row;
2141 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2144 columnIndices[0] = row-1;
2145 columnIndices[1] = row;
2146 columnIndices[2] = row+1;
2147 vals[0] =
Scalar(1.0) * a_val;
2148 vals[1] =
Scalar(-2.0) * a_val;
2149 vals[2] =
Scalar(1.0) * a_val;
2150 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2153 matrix->fillComplete();
2156 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2157 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2160 b_val.fastAccessCoeff(
j) =
2161 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
2163 for (
size_t i=0; i<num_my_row; ++i) {
2165 if (row == 0 || row == nrow-1)
2168 b_view[i] = b_val * (h*h);
2172 typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
2173 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2174 std::string solver_name;
2175 #if defined(HAVE_AMESOS2_BASKER)
2176 solver_name =
"basker";
2177 #elif defined(HAVE_AMESOS2_KLU2)
2178 solver_name =
"klu2";
2179 #elif defined(HAVE_AMESOS2_SUPERLUDIST)
2180 solver_name =
"superlu_dist";
2181 #elif defined(HAVE_AMESOS2_SUPERLUMT)
2182 solver_name =
"superlu_mt";
2183 #elif defined(HAVE_AMESOS2_SUPERLU)
2184 solver_name =
"superlu";
2185 #elif defined(HAVE_AMESOS2_PARDISO_MKL)
2186 solver_name =
"pardisomkl";
2187 #elif defined(HAVE_AMESOS2_LAPACK)
2188 solver_name =
"lapack";
2189 #elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
2190 solver_name =
"lapack";
2196 out <<
"Solving linear system with " << solver_name << std::endl;
2197 RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2198 solver_name, matrix, x, b);
2205 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
2206 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_MultiVector;
2207 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2208 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2209 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2211 Stokhos::create_flat_pce_graph(*graph, cijk, flat_x_map, flat_b_map,
2212 cijk_graph, pce_size);
2213 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2215 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2216 RCP<Flat_Tpetra_Vector> flat_x =
2218 RCP<Flat_Tpetra_Vector> flat_b =
2220 typedef Amesos2::Solver<Flat_Tpetra_CrsMatrix,Flat_Tpetra_MultiVector> Flat_Solver;
2221 RCP<Flat_Solver> flat_solver =
2222 Amesos2::create<Flat_Tpetra_CrsMatrix,Flat_Tpetra_MultiVector>(
2223 solver_name, flat_matrix, flat_x, flat_b);
2224 flat_solver->solve();
2226 typedef Kokkos::Details::ArithTraits<BaseScalar> ST;
2227 typename ST::mag_type btol = 1e-12;
2228 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
2229 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2230 for (
size_t i=0; i<num_my_row; ++i) {
2235 Scalar v = x_view[i];
2236 Scalar v2 = x2_view[i];
2238 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
2239 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2240 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
2241 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
2260 #define CRSMATRIX_UQ_PCE_TESTS_SLGN(S, LO, GO, N) \
2261 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorAdd, S, LO, GO, N ) \
2262 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorDot, S, LO, GO, N ) \
2263 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorAdd, S, LO, GO, N ) \
2264 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDot, S, LO, GO, N ) \
2265 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDotSub, S, LO, GO, N ) \
2266 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixVectorMultiply, S, LO, GO, N ) \
2267 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixMultiVectorMultiply, S, LO, GO, N ) \
2268 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Flatten, S, LO, GO, N ) \
2269 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimpleCG, S, LO, GO, N ) \
2270 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, S, LO, GO, N ) \
2271 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES, S, LO, GO, N ) \
2272 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, S, LO, GO, N ) \
2273 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosCG_Muelu, S, LO, GO, N ) \
2274 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Amesos2, S, LO, GO, N )
2276 #define CRSMATRIX_UQ_PCE_TESTS_N(N) \
2277 typedef Stokhos::DeviceForNode2<N>::type Device; \
2278 typedef Stokhos::DynamicStorage<int,double,Device::execution_space> DS; \
2279 CRSMATRIX_UQ_PCE_TESTS_SLGN(DS, int, int, N)
Teuchos::RCP< Tpetra::CrsMatrix< Scalar, LO, GO, N > > build_mean_matrix(const Tpetra::CrsMatrix< Scalar, LO, GO, N > &A)
bool CG_Solve(const Matrix &A, Vector &x, const Vector &b, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
scalar generate_multi_vector_coefficient(const ordinal nFEM, const ordinal nVec, const ordinal nStoch, const ordinal iColFEM, const ordinal iVec, const ordinal iStoch)
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node)
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
#define TEST_EQUALITY_CONST(v1, v2)
Kokkos::DefaultExecutionSpace execution_space
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
std::enable_if< Kokkos::is_view_uq_pce< Kokkos::View< XD, XP...> >::value &&Kokkos::is_view_uq_pce< Kokkos::View< YD, YP...> >::value, typename Kokkos::Details::InnerProductSpaceTraits< typename Kokkos::View< XD, XP...>::non_const_value_type >::dot_type >::type dot(const Kokkos::View< XD, XP...> &x, const Kokkos::View< YD, YP...> &y)
Teuchos::RCP< Tpetra::CrsMatrix< typename Scalar::value_type, LO, GO, N > > build_mean_scalar_matrix(const Tpetra::CrsMatrix< Scalar, LO, GO, N > &A)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec_const, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
bool PCG_Solve(const Matrix &A, Vector &x, const Vector &b, const Prec &M, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
KokkosClassic::DefaultNode::DefaultNodeType Node
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
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)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &flat_graph, const LocalOrdinal block_size)
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j)-expr2.val(j)
void setGlobalCijkTensor(const cijk_type &cijk)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
kokkos_cijk_type build_cijk(ordinal_type stoch_dim, ordinal_type poly_ord)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Legendre polynomial basis.
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.
#define TEST_EQUALITY(v1, v2)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)