12 #include "Stokhos_Ensemble_Sizes.hpp"
20 #include "Tpetra_Core.hpp"
21 #include "Tpetra_Map.hpp"
22 #include "Tpetra_MultiVector.hpp"
23 #include "Tpetra_Vector.hpp"
24 #include "Tpetra_CrsGraph.hpp"
25 #include "Tpetra_CrsMatrix.hpp"
26 #include "Tpetra_Details_WrappedDualView.hpp"
30 #ifdef HAVE_STOKHOS_BELOS
32 #include "BelosLinearProblem.hpp"
33 #include "BelosPseudoBlockGmresSolMgr.hpp"
34 #include "BelosPseudoBlockCGSolMgr.hpp"
38 #ifdef HAVE_STOKHOS_IFPACK2
40 #include "Ifpack2_Factory.hpp"
44 #ifdef HAVE_STOKHOS_MUELU
46 #include "MueLu_CreateTpetraPreconditioner.hpp"
50 #ifdef HAVE_STOKHOS_AMESOS2
52 #include "Amesos2_Factory.hpp"
55 template <
typename scalar,
typename ordinal>
59 const ordinal iColFEM,
60 const ordinal iStoch )
62 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
63 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
64 return X_fem + X_stoch;
68 template <
typename scalar,
typename ordinal>
73 const ordinal iColFEM,
77 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
78 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
79 const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
80 return X_fem + X_stoch + X_col;
106 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
107 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
110 if ( !Kokkos::is_initialized() )
111 Kokkos::initialize();
114 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
118 RCP<const Tpetra_Map> map =
119 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
121 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
122 const size_t num_my_row = myGIDs.size();
125 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
126 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
128 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
129 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
131 for (
size_t i=0; i<num_my_row; ++i) {
134 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
135 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
145 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
146 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
152 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
154 BaseScalar tol = 1.0e-14;
155 for (
size_t i=0; i<num_my_row; ++i) {
158 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
160 val.fastAccessCoeff(
j) = alpha.coeff(
j)*v + 0.12345*beta.coeff(
j)*v;
184 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
185 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
186 typedef typename Tpetra_Vector::dot_type dot_type;
189 if ( !Kokkos::is_initialized() )
190 Kokkos::initialize();
193 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
197 RCP<const Tpetra_Map> map =
198 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
200 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
201 const size_t num_my_row = myGIDs.size();
204 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
205 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
207 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
208 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
210 for (
size_t i=0; i<num_my_row; ++i) {
213 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
214 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
222 dot_type
dot = x1->dot(*x2);
226 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
229 dot_type local_val(0);
230 for (
size_t i=0; i<num_my_row; ++i) {
233 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
235 local_val += 0.12345 * v * v;
242 Teuchos::outArg(val));
248 for (
size_t i=0; i<num_my_row; ++i) {
251 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
253 local_val.fastAccessCoeff(
j) += 0.12345 * v * v;
260 Teuchos::outArg(val));
264 out <<
"dot = " << dot <<
" expected = " << val << std::endl;
266 BaseScalar tol = 1.0e-14;
286 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
287 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
290 if ( !Kokkos::is_initialized() )
291 Kokkos::initialize();
294 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
298 RCP<const Tpetra_Map> map =
299 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
301 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
302 const size_t num_my_row = myGIDs.size();
306 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
307 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
309 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
310 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
312 for (
size_t i=0; i<num_my_row; ++i) {
314 for (
size_t j=0;
j<ncol; ++
j) {
317 generate_multi_vector_coefficient<BaseScalar,size_t>(
319 val1.fastAccessCoeff(k) = v;
320 val2.fastAccessCoeff(k) = 0.12345 * v;
331 RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
332 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
338 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
340 BaseScalar tol = 1.0e-14;
341 for (
size_t i=0; i<num_my_row; ++i) {
343 for (
size_t j=0;
j<ncol; ++
j) {
345 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
347 val.fastAccessCoeff(k) = alpha.coeff(k)*v + 0.12345*beta.coeff(k)*v;
352 val.fastAccessCoeff(k), tol );
373 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
374 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
375 typedef typename Tpetra_MultiVector::dot_type dot_type;
378 if ( !Kokkos::is_initialized() )
379 Kokkos::initialize();
382 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
386 RCP<const Tpetra_Map> map =
387 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
389 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
390 const size_t num_my_row = myGIDs.size();
394 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
395 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
397 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
398 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
400 for (
size_t i=0; i<num_my_row; ++i) {
402 for (
size_t j=0;
j<ncol; ++
j) {
405 generate_multi_vector_coefficient<BaseScalar,size_t>(
407 val1.fastAccessCoeff(k) = v;
408 val2.fastAccessCoeff(k) = 0.12345 * v;
417 Array<dot_type> dots(ncol);
418 x1->dot(*x2, dots());
422 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
425 Array<dot_type> local_vals(ncol, dot_type(0));
426 for (
size_t i=0; i<num_my_row; ++i) {
428 for (
size_t j=0;
j<ncol; ++
j) {
430 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
432 local_vals[
j] += 0.12345 * v * v;
438 Array<dot_type> vals(ncol, dot_type(0));
440 local_vals.getRawPtr(), vals.getRawPtr());
445 Array<dot_type> local_vals(ncol, dot_type(
VectorSize, 0.0));
446 for (
size_t i=0; i<num_my_row; ++i) {
448 for (
size_t j=0;
j<ncol; ++
j) {
450 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
452 local_vals[
j].fastAccessCoeff(k) += 0.12345 * v * v;
458 Array<dot_type> vals(ncol, dot_type(
VectorSize, 0.0));
460 local_vals.getRawPtr(), vals.getRawPtr());
464 BaseScalar tol = 1.0e-14;
465 for (
size_t j=0;
j<ncol; ++
j) {
466 out <<
"dots(" <<
j <<
") = " << dots[
j]
467 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
488 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
489 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
490 typedef typename Tpetra_MultiVector::dot_type dot_type;
493 if ( !Kokkos::is_initialized() )
494 Kokkos::initialize();
497 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
501 RCP<const Tpetra_Map> map =
502 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
504 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
505 const size_t num_my_row = myGIDs.size();
509 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
510 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
512 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
513 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
515 for (
size_t i=0; i<num_my_row; ++i) {
517 for (
size_t j=0;
j<ncol; ++
j) {
520 generate_multi_vector_coefficient<BaseScalar,size_t>(
522 val1.fastAccessCoeff(k) = v;
523 val2.fastAccessCoeff(k) = 0.12345 * v;
534 cols[0] = 4; cols[1] = 2;
535 RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
536 RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
539 Array<dot_type> dots(ncol_sub);
540 x1_sub->dot(*x2_sub, dots());
544 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
547 Array<dot_type> local_vals(ncol_sub, dot_type(0));
548 for (
size_t i=0; i<num_my_row; ++i) {
550 for (
size_t j=0;
j<ncol_sub; ++
j) {
552 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
554 local_vals[
j] += 0.12345 * v * v;
560 Array<dot_type> vals(ncol_sub, dot_type(0));
562 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
568 Array<dot_type> local_vals(ncol_sub, dot_type(
VectorSize, 0.0));
569 for (
size_t i=0; i<num_my_row; ++i) {
571 for (
size_t j=0;
j<ncol_sub; ++
j) {
573 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
575 local_vals[
j].fastAccessCoeff(k) += 0.12345 * v * v;
581 Array<dot_type> vals(ncol_sub, dot_type(
VectorSize, 0.0));
583 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
588 BaseScalar tol = 1.0e-14;
589 for (
size_t j=0;
j<ncol_sub; ++
j) {
590 out <<
"dots(" <<
j <<
") = " << dots[
j]
591 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
612 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
613 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
614 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
615 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
618 if ( !Kokkos::is_initialized() )
619 Kokkos::initialize();
623 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
624 RCP<const Tpetra_Map> map =
625 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
627 RCP<Tpetra_CrsGraph> graph =
628 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
629 Array<GlobalOrdinal> columnIndices(2);
630 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
631 const size_t num_my_row = myGIDs.size();
632 for (
size_t i=0; i<num_my_row; ++i) {
634 columnIndices[0] = row;
637 columnIndices[1] = row+1;
640 graph->insertGlobalIndices(row, columnIndices(0,ncol));
642 graph->fillComplete();
643 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
646 Array<Scalar> vals(2);
648 for (
size_t i=0; i<num_my_row; ++i) {
650 columnIndices[0] = row;
652 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
658 columnIndices[1] = row+1;
660 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
665 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
667 matrix->fillComplete();
670 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
672 auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
673 for (
size_t i=0; i<num_my_row; ++i) {
676 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
689 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
690 matrix->apply(*x, *y);
696 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
697 BaseScalar tol = 1.0e-14;
698 for (
size_t i=0; i<num_my_row; ++i) {
701 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
703 val.fastAccessCoeff(
j) = v*v;
707 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
709 val.fastAccessCoeff(
j) += v*v;
734 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
735 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
736 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
737 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
740 if ( !Kokkos::is_initialized() )
741 Kokkos::initialize();
745 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
746 RCP<const Tpetra_Map> map =
747 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
749 RCP<Tpetra_CrsGraph> graph =
750 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
751 Array<GlobalOrdinal> columnIndices(2);
752 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
753 const size_t num_my_row = myGIDs.size();
754 for (
size_t i=0; i<num_my_row; ++i) {
756 columnIndices[0] = row;
759 columnIndices[1] = row+1;
762 graph->insertGlobalIndices(row, columnIndices(0,ncol));
764 graph->fillComplete();
765 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
768 Array<Scalar> vals(2);
770 for (
size_t i=0; i<num_my_row; ++i) {
772 columnIndices[0] = row;
774 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
780 columnIndices[1] = row+1;
782 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
787 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
789 matrix->fillComplete();
793 RCP<Tpetra_MultiVector> x = Tpetra::createMultiVector<Scalar>(map, ncol);
795 auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
796 for (
size_t i=0; i<num_my_row; ++i) {
798 for (
size_t j=0;
j<ncol; ++
j) {
801 generate_multi_vector_coefficient<BaseScalar,size_t>(
803 val.fastAccessCoeff(k) = v;
817 RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
818 matrix->apply(*x, *y);
824 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
825 BaseScalar tol = 1.0e-14;
826 for (
size_t i=0; i<num_my_row; ++i) {
828 for (
size_t j=0;
j<ncol; ++
j) {
830 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
832 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
834 val.fastAccessCoeff(k) = v1*v2;
838 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
840 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
842 val.fastAccessCoeff(k) += v1*v2;
848 val.fastAccessCoeff(k), tol );
869 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
870 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
871 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
872 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
874 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
875 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
878 if ( !Kokkos::is_initialized() )
879 Kokkos::initialize();
883 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
884 RCP<const Tpetra_Map> map =
885 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
887 RCP<Tpetra_CrsGraph> graph =
888 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
889 Array<GlobalOrdinal> columnIndices(2);
890 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
891 const size_t num_my_row = myGIDs.size();
892 for (
size_t i=0; i<num_my_row; ++i) {
894 columnIndices[0] = row;
897 columnIndices[1] = row+1;
900 graph->insertGlobalIndices(row, columnIndices(0,ncol));
902 graph->fillComplete();
903 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
906 Array<Scalar> vals(2);
908 for (
size_t i=0; i<num_my_row; ++i) {
910 columnIndices[0] = row;
912 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
918 columnIndices[1] = row+1;
920 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
925 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
927 matrix->fillComplete();
930 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
932 auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
933 for (
size_t i=0; i<num_my_row; ++i) {
936 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
943 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
944 matrix->apply(*x, *y);
947 RCP<const Tpetra_Map> flat_x_map, flat_y_map;
948 RCP<const Tpetra_CrsGraph> flat_graph =
950 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
954 RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
956 RCP<Flat_Tpetra_Vector> flat_x =
958 RCP<Flat_Tpetra_Vector> flat_y =
960 flat_matrix->apply(*flat_x, *flat_y);
967 BaseScalar tol = 1.0e-14;
968 auto y_view = y-> getLocalViewHost(Tpetra::Access::ReadOnly);
969 auto y2_view = y2->getLocalViewHost(Tpetra::Access::ReadOnly);
970 for (
size_t i=0; i<num_my_row; ++i) {
999 using DualViewType = Kokkos::DualView<Scalar*, typename Node::device_type>;
1000 using WDV = Tpetra::Details::WrappedDualView<DualViewType>;
1001 using values_view =
typename DualViewType::t_dev;
1004 if ( !Kokkos::is_initialized() )
1005 Kokkos::initialize();
1009 values_view myView(
"emptyTestView", 0);
1012 size_t use_h = wdv.getHostView(Tpetra::Access::ReadOnly).use_count();
1013 size_t use_d = wdv.getDeviceView(Tpetra::Access::ReadOnly).use_count();
1036 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1037 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1038 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1039 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1042 if ( !Kokkos::is_initialized() )
1043 Kokkos::initialize();
1047 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1048 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1049 RCP<const Tpetra_Map> map =
1050 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1052 RCP<Tpetra_CrsGraph> graph =
1053 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
1054 Array<GlobalOrdinal> columnIndices(3);
1055 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1056 const size_t num_my_row = myGIDs.size();
1057 for (
size_t i=0; i<num_my_row; ++i) {
1059 if (row == 0 || row == nrow-1) {
1060 columnIndices[0] = row;
1061 graph->insertGlobalIndices(row, columnIndices(0,1));
1064 columnIndices[0] = row-1;
1065 columnIndices[1] = row;
1066 columnIndices[2] = row+1;
1067 graph->insertGlobalIndices(row, columnIndices(0,3));
1070 graph->fillComplete();
1071 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1074 Array<Scalar> vals(3);
1077 a_val.fastAccessCoeff(
j) =
1078 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(VectorSize);
1080 for (
size_t i=0; i<num_my_row; ++i) {
1082 if (row == 0 || row == nrow-1) {
1083 columnIndices[0] = row;
1085 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1088 columnIndices[0] = row-1;
1089 columnIndices[1] = row;
1090 columnIndices[2] = row+1;
1091 vals[0] =
Scalar(-1.0) * a_val;
1092 vals[1] =
Scalar(2.0) * a_val;
1093 vals[2] =
Scalar(-1.0) * a_val;
1094 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1097 matrix->fillComplete();
1099 matrix->describe(*(Teuchos::fancyOStream(
rcp(&std::cout,
false))),
1103 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1106 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1107 b_val =
Scalar(VectorSize, BaseScalar(0.0));
1109 b_val.fastAccessCoeff(
j) =
1110 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(VectorSize);
1112 for (
size_t i=0; i<num_my_row; ++i) {
1114 if (row == 0 || row == nrow-1)
1115 b_view(i,0) =
Scalar(0.0);
1117 b_view(i,0) = -
Scalar(b_val * h * h);
1122 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1123 #if KOKKOS_VERSION >= 40799
1124 typedef KokkosKernels::ArithTraits<BaseScalar> BST;
1126 typedef Kokkos::ArithTraits<BaseScalar> BST;
1128 typedef typename BST::mag_type base_mag_type;
1129 typedef typename Tpetra_Vector::mag_type mag_type;
1130 base_mag_type btol = 1e-9;
1131 mag_type tol = btol;
1134 out.getOStream().get());
1142 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1143 Scalar
val(VectorSize, BaseScalar(0.0));
1144 for (
size_t i=0; i<num_my_row; ++i) {
1146 BaseScalar xx = row * h;
1148 val.fastAccessCoeff(
j) =
1149 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1154 Scalar v = x_view(i,0);
1157 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1159 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1168 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2)
1182 using Teuchos::getParametersFromXmlFile;
1188 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1189 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1190 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1191 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1194 if ( !Kokkos::is_initialized() )
1195 Kokkos::initialize();
1199 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1200 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1201 RCP<const Tpetra_Map> map =
1202 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1204 RCP<Tpetra_CrsGraph> graph =
1205 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
1206 Array<GlobalOrdinal> columnIndices(3);
1207 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1208 const size_t num_my_row = myGIDs.size();
1209 for (
size_t i=0; i<num_my_row; ++i) {
1211 if (row == 0 || row == nrow-1) {
1212 columnIndices[0] = row;
1213 graph->insertGlobalIndices(row, columnIndices(0,1));
1216 columnIndices[0] = row-1;
1217 columnIndices[1] = row;
1218 columnIndices[2] = row+1;
1219 graph->insertGlobalIndices(row, columnIndices(0,3));
1222 graph->fillComplete();
1223 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1226 Array<Scalar> vals(3);
1229 a_val.fastAccessCoeff(
j) =
1230 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(VectorSize);
1232 for (
size_t i=0; i<num_my_row; ++i) {
1234 if (row == 0 || row == nrow-1) {
1235 columnIndices[0] = row;
1237 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1240 columnIndices[0] = row-1;
1241 columnIndices[1] = row;
1242 columnIndices[2] = row+1;
1243 vals[0] =
Scalar(-1.0) * a_val;
1244 vals[1] =
Scalar(2.0) * a_val;
1245 vals[2] =
Scalar(-1.0) * a_val;
1246 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1249 matrix->fillComplete();
1252 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1255 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1256 b_val =
Scalar(VectorSize, BaseScalar(0.0));
1258 b_val.fastAccessCoeff(
j) =
1259 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(VectorSize);
1261 for (
size_t i=0; i<num_my_row; ++i) {
1263 if (row == 0 || row == nrow-1)
1264 b_view(i,0) =
Scalar(0.0);
1266 b_view(i,0) = -
Scalar(b_val * h * h);
1271 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1272 RCP<OP> matrix_op = matrix;
1273 RCP<ParameterList> muelu_params =
1274 getParametersFromXmlFile(
"muelu_cheby.xml");
1276 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1279 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1280 #if KOKKOS_VERSION >= 40799
1281 typedef KokkosKernels::ArithTraits<BaseScalar> BST;
1283 typedef Kokkos::ArithTraits<BaseScalar> BST;
1285 typedef typename BST::mag_type base_mag_type;
1286 typedef typename Tpetra_Vector::mag_type mag_type;
1287 base_mag_type btol = 1e-9;
1288 mag_type tol = btol;
1291 out.getOStream().get());
1299 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1300 Scalar
val(VectorSize, BaseScalar(0.0));
1301 for (
size_t i=0; i<num_my_row; ++i) {
1303 BaseScalar xx = row * h;
1305 val.fastAccessCoeff(
j) =
1306 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1311 Scalar v = x_view(i,0);
1313 if (BST::magnitude(v.coeff(
j)) < btol)
1314 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1315 if (BST::magnitude(
val.coeff(
j)) < btol)
1316 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1332 #if defined(HAVE_STOKHOS_BELOS)
1351 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1352 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1353 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1354 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1357 if ( !Kokkos::is_initialized() )
1358 Kokkos::initialize();
1362 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1363 RCP<const Tpetra_Map> map =
1364 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1366 RCP<Tpetra_CrsGraph> graph =
1367 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
1368 Array<GlobalOrdinal> columnIndices(2);
1369 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1370 const size_t num_my_row = myGIDs.size();
1371 for (
size_t i=0; i<num_my_row; ++i) {
1373 columnIndices[0] = row;
1375 if (row != nrow-1) {
1376 columnIndices[1] = row+1;
1379 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1381 graph->fillComplete();
1382 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1385 Array<Scalar> vals(2);
1386 Scalar
val(VectorSize, BaseScalar(0.0));
1387 for (
size_t i=0; i<num_my_row; ++i) {
1389 columnIndices[0] = row;
1391 val.fastAccessCoeff(
j) =
j+1;
1395 if (row != nrow-1) {
1396 columnIndices[1] = row+1;
1398 val.fastAccessCoeff(
j) =
j+1;
1402 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1404 matrix->fillComplete();
1407 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1409 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1410 for (
size_t i=0; i<num_my_row; ++i) {
1411 b_view(i,0) =
Scalar(1.0);
1417 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1418 typedef BaseScalar BelosScalar;
1420 typedef Scalar BelosScalar;
1422 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1423 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1424 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1425 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1426 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
1427 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
1428 typename ST::magnitudeType tol = 1e-9;
1429 belosParams->set(
"Flexible Gmres",
false);
1430 belosParams->set(
"Num Blocks", 100);
1431 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1432 belosParams->set(
"Maximum Iterations", 100);
1433 belosParams->set(
"Verbosity", 33);
1434 belosParams->set(
"Output Style", 1);
1435 belosParams->set(
"Output Frequency", 1);
1436 belosParams->set(
"Output Stream", out.getOStream());
1437 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1438 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1439 problem->setProblem();
1440 Belos::ReturnType ret = solver->solve();
1453 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1454 for (
size_t i=0; i<num_my_row; ++i) {
1458 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
1462 val =
Scalar(VectorSize, BaseScalar(0.0));
1466 Scalar v = x_view(i,0);
1468 if (ST::magnitude(v.coeff(
j)) < tol)
1469 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1485 using Teuchos::tuple;
1495 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1496 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1497 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1498 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1501 if ( !Kokkos::is_initialized() )
1502 Kokkos::initialize();
1506 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1507 RCP<const Tpetra_Map> map =
1508 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1510 RCP<Tpetra_CrsGraph> graph =
1511 rcp(
new Tpetra_CrsGraph(map,
size_t(1)));
1512 Array<GlobalOrdinal> columnIndices(1);
1513 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1514 const size_t num_my_row = myGIDs.size();
1515 for (
size_t i=0; i<num_my_row; ++i) {
1517 columnIndices[0] = row;
1519 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1521 graph->fillComplete();
1522 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1524 Array<Scalar> vals(1);
1525 for (
size_t i=0; i<num_my_row; ++i) {
1527 columnIndices[0] = row;
1529 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1531 matrix->fillComplete();
1540 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1542 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1543 for (
size_t i=0; i<num_my_row; ++i) {
1545 b_view(i,0) =
Scalar(0.0);
1547 if (
int(
j+2+row-VectorSize) > 0)
1548 b_view(i,0).fastAccessCoeff(
j) = BaseScalar(row+1);
1554 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1555 typedef BaseScalar BelosScalar;
1557 typedef Scalar BelosScalar;
1559 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1560 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1561 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1562 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1563 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
1564 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
1565 typename ST::magnitudeType tol = 1e-9;
1566 belosParams->set(
"Flexible Gmres",
false);
1567 belosParams->set(
"Num Blocks", 100);
1568 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1569 belosParams->set(
"Maximum Iterations", 100);
1570 belosParams->set(
"Verbosity", 33);
1571 belosParams->set(
"Output Style", 1);
1572 belosParams->set(
"Output Frequency", 1);
1573 belosParams->set(
"Output Stream", out.getOStream());
1574 belosParams->set(
"Orthogonalization",
"DGKS");
1575 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1576 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1577 problem->setProblem();
1578 Belos::ReturnType ret = solver->solve();
1581 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1582 int numItersExpected = nrow;
1583 int numIters = solver->getNumIters();
1584 out <<
"numIters = " << numIters << std::endl;
1588 std::vector<int> ensemble_iterations =
1590 out <<
"Ensemble iterations = ";
1591 for (
auto ensemble_iteration : ensemble_iterations)
1592 out << ensemble_iteration <<
" ";
1596 if (
int(
j+1+nrow-VectorSize) > 0) {
1612 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1613 for (
size_t i=0; i<num_my_row; ++i) {
1615 Scalar v = x_view(i,0);
1618 if (ST::magnitude(v.coeff(
j)) < tol)
1619 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1625 if (
j+2+row-VectorSize > 0)
1626 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1641 using Teuchos::tuple;
1651 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1652 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1653 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1654 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1657 if ( !Kokkos::is_initialized() )
1658 Kokkos::initialize();
1662 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1663 RCP<const Tpetra_Map> map =
1664 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1666 RCP<Tpetra_CrsGraph> graph =
1667 rcp(
new Tpetra_CrsGraph(map,
size_t(1)));
1668 Array<GlobalOrdinal> columnIndices(1);
1669 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1670 const size_t num_my_row = myGIDs.size();
1671 for (
size_t i=0; i<num_my_row; ++i) {
1673 columnIndices[0] = row;
1675 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1677 graph->fillComplete();
1678 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1680 Array<Scalar> vals(1);
1681 for (
size_t i=0; i<num_my_row; ++i) {
1683 columnIndices[0] = row;
1685 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1687 matrix->fillComplete();
1696 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1698 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1699 for (
size_t i=0; i<num_my_row; ++i) {
1701 b_view(i,0) =
Scalar(0.0);
1703 if (
int(
j+2+row-VectorSize) > 0)
1704 b_view(i,0).fastAccessCoeff(
j) = BaseScalar(row+1);
1710 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1711 typedef BaseScalar BelosScalar;
1713 typedef Scalar BelosScalar;
1715 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1716 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1717 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1718 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1719 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
1720 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
1721 typename ST::magnitudeType tol = 1e-9;
1722 belosParams->set(
"Flexible Gmres",
false);
1723 belosParams->set(
"Num Blocks", 100);
1724 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1725 belosParams->set(
"Maximum Iterations", 100);
1726 belosParams->set(
"Verbosity", 33);
1727 belosParams->set(
"Output Style", 1);
1728 belosParams->set(
"Output Frequency", 1);
1729 belosParams->set(
"Output Stream", out.getOStream());
1730 belosParams->set(
"Orthogonalization",
"ICGS");
1731 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1732 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1733 problem->setProblem();
1734 Belos::ReturnType ret = solver->solve();
1737 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1738 int numItersExpected = nrow;
1739 int numIters = solver->getNumIters();
1740 out <<
"numIters = " << numIters << std::endl;
1744 std::vector<int> ensemble_iterations =
1746 out <<
"Ensemble iterations = ";
1747 for (
auto ensemble_iteration : ensemble_iterations)
1748 out << ensemble_iteration <<
" ";
1752 if (
int(
j+1+nrow-VectorSize) > 0) {
1768 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1769 for (
size_t i=0; i<num_my_row; ++i) {
1771 Scalar v = x_view(i,0);
1774 if (ST::magnitude(v.coeff(
j)) < tol)
1775 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1778 Scalar val =
Scalar(0.0);
1781 if (
j+2+row-VectorSize > 0)
1782 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1797 using Teuchos::tuple;
1807 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1808 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1809 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1810 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1813 if ( !Kokkos::is_initialized() )
1814 Kokkos::initialize();
1818 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1819 RCP<const Tpetra_Map> map =
1820 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1822 RCP<Tpetra_CrsGraph> graph =
1823 rcp(
new Tpetra_CrsGraph(map,
size_t(1)));
1824 Array<GlobalOrdinal> columnIndices(1);
1825 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1826 const size_t num_my_row = myGIDs.size();
1827 for (
size_t i=0; i<num_my_row; ++i) {
1829 columnIndices[0] = row;
1831 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1833 graph->fillComplete();
1834 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
1836 Array<Scalar> vals(1);
1837 for (
size_t i=0; i<num_my_row; ++i) {
1839 columnIndices[0] = row;
1841 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1843 matrix->fillComplete();
1852 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1854 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1855 for (
size_t i=0; i<num_my_row; ++i) {
1857 b_view(i,0) =
Scalar(0.0);
1859 if (
int(
j+2+row-VectorSize) > 0)
1860 b_view(i,0).fastAccessCoeff(
j) = BaseScalar(row+1);
1866 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1867 typedef BaseScalar BelosScalar;
1869 typedef Scalar BelosScalar;
1871 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1872 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1873 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1874 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1875 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
1876 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
1877 typename ST::magnitudeType tol = 1e-9;
1878 belosParams->set(
"Flexible Gmres",
false);
1879 belosParams->set(
"Num Blocks", 100);
1880 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1881 belosParams->set(
"Maximum Iterations", 100);
1882 belosParams->set(
"Verbosity", 33);
1883 belosParams->set(
"Output Style", 1);
1884 belosParams->set(
"Output Frequency", 1);
1885 belosParams->set(
"Output Stream", out.getOStream());
1886 belosParams->set(
"Orthogonalization",
"IMGS");
1887 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1888 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1889 problem->setProblem();
1890 Belos::ReturnType ret = solver->solve();
1893 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1894 int numItersExpected = nrow;
1895 int numIters = solver->getNumIters();
1896 out <<
"numIters = " << numIters << std::endl;
1900 std::vector<int> ensemble_iterations =
1902 out <<
"Ensemble iterations = ";
1903 for (
auto ensemble_iteration : ensemble_iterations)
1904 out << ensemble_iteration <<
" ";
1908 if (
int(
j+1+nrow-VectorSize) > 0) {
1924 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1925 for (
size_t i=0; i<num_my_row; ++i) {
1927 Scalar v = x_view(i,0);
1930 if (ST::magnitude(v.coeff(
j)) < tol)
1931 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1934 Scalar val =
Scalar(0.0);
1937 if (
j+2+row-VectorSize > 0)
1938 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1965 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2)
1985 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1986 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1987 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1988 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1991 if ( !Kokkos::is_initialized() )
1992 Kokkos::initialize();
1996 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1997 RCP<const Tpetra_Map> map =
1998 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2000 RCP<Tpetra_CrsGraph> graph =
2001 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
2002 Array<GlobalOrdinal> columnIndices(2);
2003 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2004 const size_t num_my_row = myGIDs.size();
2005 for (
size_t i=0; i<num_my_row; ++i) {
2007 columnIndices[0] = row;
2009 if (row != nrow-1) {
2010 columnIndices[1] = row+1;
2013 graph->insertGlobalIndices(row, columnIndices(0,ncol));
2015 graph->fillComplete();
2016 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
2019 Array<Scalar> vals(2);
2020 Scalar
val(VectorSize, BaseScalar(0.0));
2021 for (
size_t i=0; i<num_my_row; ++i) {
2023 columnIndices[0] = row;
2025 val.fastAccessCoeff(
j) =
j+1;
2029 if (row != nrow-1) {
2030 columnIndices[1] = row+1;
2032 val.fastAccessCoeff(
j) =
j+1;
2036 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
2038 matrix->fillComplete();
2041 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2043 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
2044 for (
size_t i=0; i<num_my_row; ++i) {
2045 b_view(i,0) =
Scalar(1.0);
2050 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
2051 Ifpack2::Factory factory;
2052 RCP<Prec> M = factory.create<Tpetra_CrsMatrix>(
"RILUK", matrix);
2058 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
2059 typedef BaseScalar BelosScalar;
2061 typedef Scalar BelosScalar;
2063 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2064 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2065 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2066 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2067 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
2068 problem->setRightPrec(M);
2069 problem->setProblem();
2070 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
2071 typename ST::magnitudeType tol = 1e-9;
2072 belosParams->set(
"Flexible Gmres",
false);
2073 belosParams->set(
"Num Blocks", 100);
2074 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
2075 belosParams->set(
"Maximum Iterations", 100);
2076 belosParams->set(
"Verbosity", 33);
2077 belosParams->set(
"Output Style", 1);
2078 belosParams->set(
"Output Frequency", 1);
2079 belosParams->set(
"Output Stream", out.getOStream());
2081 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2082 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2083 Belos::ReturnType ret = solver->solve();
2096 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
2097 for (
size_t i=0; i<num_my_row; ++i) {
2101 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
2105 val =
Scalar(VectorSize, BaseScalar(0.0));
2109 Scalar v = x_view(i,0);
2111 if (ST::magnitude(v.coeff(
j)) < tol)
2112 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2128 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2)
2142 using Teuchos::getParametersFromXmlFile;
2148 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2149 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2150 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2151 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2154 if ( !Kokkos::is_initialized() )
2155 Kokkos::initialize();
2159 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
2160 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2161 RCP<const Tpetra_Map> map =
2162 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2164 RCP<Tpetra_CrsGraph> graph =
2165 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
2166 Array<GlobalOrdinal> columnIndices(3);
2167 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2168 const size_t num_my_row = myGIDs.size();
2169 for (
size_t i=0; i<num_my_row; ++i) {
2171 if (row == 0 || row == nrow-1) {
2172 columnIndices[0] = row;
2173 graph->insertGlobalIndices(row, columnIndices(0,1));
2176 columnIndices[0] = row-1;
2177 columnIndices[1] = row;
2178 columnIndices[2] = row+1;
2179 graph->insertGlobalIndices(row, columnIndices(0,3));
2182 graph->fillComplete();
2183 RCP<Tpetra_CrsMatrix> matrix =
rcp(
new Tpetra_CrsMatrix(graph));
2186 Array<Scalar> vals(3);
2187 Scalar a_val(VectorSize, BaseScalar(0.0));
2189 a_val.fastAccessCoeff(
j) =
2190 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(VectorSize);
2192 for (
size_t i=0; i<num_my_row; ++i) {
2194 if (row == 0 || row == nrow-1) {
2195 columnIndices[0] = row;
2197 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2200 columnIndices[0] = row-1;
2201 columnIndices[1] = row;
2202 columnIndices[2] = row+1;
2203 vals[0] =
Scalar(-1.0) * a_val;
2204 vals[1] =
Scalar(2.0) * a_val;
2205 vals[2] =
Scalar(-1.0) * a_val;
2206 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2209 matrix->fillComplete();
2212 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2215 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
2216 b_val =
Scalar(VectorSize, BaseScalar(0.0));
2218 b_val.fastAccessCoeff(
j) =
2219 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(VectorSize);
2221 for (
size_t i=0; i<num_my_row; ++i) {
2223 if (row == 0 || row == nrow-1)
2224 b_view(i,0) =
Scalar(0.0);
2226 b_view(i,0) = -
Scalar(b_val * h * h);
2231 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2232 RCP<ParameterList> muelu_params =
2233 getParametersFromXmlFile(
"muelu_cheby.xml");
2234 RCP<OP> matrix_op = matrix;
2236 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
2240 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
2241 typedef BaseScalar BelosScalar;
2243 typedef Scalar BelosScalar;
2245 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2246 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2247 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
2248 RCP< BLinProb > problem =
rcp(
new BLinProb(matrix, x, b));
2249 problem->setRightPrec(M);
2250 problem->setProblem();
2251 RCP<ParameterList> belosParams =
rcp(
new ParameterList);
2252 typename ST::magnitudeType tol = 1e-9;
2253 belosParams->set(
"Num Blocks", 100);
2254 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
2255 belosParams->set(
"Maximum Iterations", 100);
2256 belosParams->set(
"Verbosity", 33);
2257 belosParams->set(
"Output Style", 1);
2258 belosParams->set(
"Output Frequency", 1);
2259 belosParams->set(
"Output Stream", out.getOStream());
2262 belosParams->set(
"Implicit Residual Scaling",
"None");
2264 RCP<Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true> > solver =
2265 rcp(
new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true>(problem, belosParams));
2266 Belos::ReturnType ret = solver->solve();
2269 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
2271 std::vector<int> ensemble_iterations =
2272 solver->getResidualStatusTest()->getEnsembleIterations();
2273 out <<
"Ensemble iterations = ";
2275 out << ensemble_iterations[i] <<
" ";
2284 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
2285 Scalar
val(VectorSize, BaseScalar(0.0));
2286 for (
size_t i=0; i<num_my_row; ++i) {
2288 BaseScalar xx = row * h;
2290 val.fastAccessCoeff(
j) =
2291 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
2296 Scalar v = x_view(i,0);
2298 if (ST::magnitude(v.coeff(
j)) < tol)
2299 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2300 if (ST::magnitude(val.coeff(
j)) < tol)
2301 val.fastAccessCoeff(
j) = BaseScalar(0.0);
2318 #if defined(HAVE_STOKHOS_AMESOS2)
2323 #define TEST_AMESOS2_SOLVER(SolverName) \
2324 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, Amesos2_##SolverName, Storage, LocalOrdinal, GlobalOrdinal, Node) \
2326 using Teuchos::RCP; \
2327 using Teuchos::rcp; \
2328 using Teuchos::ArrayView; \
2329 using Teuchos::Array; \
2330 using Teuchos::ArrayRCP; \
2331 using Teuchos::ParameterList; \
2333 typedef typename Storage::value_type BaseScalar; \
2334 typedef Sacado::MP::Vector<Storage> Scalar; \
2336 typedef Teuchos::Comm<int> Tpetra_Comm; \
2337 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map; \
2338 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector; \
2339 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector; \
2340 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix; \
2341 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph; \
2344 if ( !Kokkos::is_initialized() ) \
2345 Kokkos::initialize(); \
2348 GlobalOrdinal nrow = 50; \
2349 BaseScalar h = 1.0 / static_cast<BaseScalar>(nrow-1); \
2350 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm(); \
2351 RCP<const Tpetra_Map> map = \
2352 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>( \
2354 RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map, size_t(3)); \
2355 Array<GlobalOrdinal> columnIndices(3); \
2356 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList(); \
2357 const size_t num_my_row = myGIDs.size(); \
2358 for (size_t i=0; i<num_my_row; ++i) { \
2359 const GlobalOrdinal row = myGIDs[i]; \
2360 if (row == 0 || row == nrow-1) { \
2361 columnIndices[0] = row; \
2362 graph->insertGlobalIndices(row, columnIndices(0,1)); \
2365 columnIndices[0] = row-1; \
2366 columnIndices[1] = row; \
2367 columnIndices[2] = row+1; \
2368 graph->insertGlobalIndices(row, columnIndices(0,3)); \
2371 graph->fillComplete(); \
2372 RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph)); \
2375 Array<Scalar> vals(3); \
2376 Scalar a_val(VectorSize, BaseScalar(0.0)); \
2377 for (LocalOrdinal j=0; j<VectorSize; ++j) { \
2378 a_val.fastAccessCoeff(j) = \
2379 BaseScalar(1.0) + BaseScalar(j) / BaseScalar(VectorSize); \
2381 for (size_t i=0; i<num_my_row; ++i) { \
2382 const GlobalOrdinal row = myGIDs[i]; \
2383 if (row == 0 || row == nrow-1) { \
2384 columnIndices[0] = row; \
2385 vals[0] = Scalar(1.0); \
2386 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1)); \
2389 columnIndices[0] = row-1; \
2390 columnIndices[1] = row; \
2391 columnIndices[2] = row+1; \
2392 vals[0] = Scalar(-1.0) * a_val; \
2393 vals[1] = Scalar(2.0) * a_val; \
2394 vals[2] = Scalar(-1.0) * a_val; \
2395 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3)); \
2398 matrix->fillComplete();\
2401 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map); \
2404 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll); \
2405 b_val = Scalar(VectorSize, BaseScalar(0.0)); \
2406 for (LocalOrdinal j=0; j<VectorSize; ++j) { \
2407 b_val.fastAccessCoeff(j) = \
2408 BaseScalar(-1.0) + BaseScalar(j) / BaseScalar(VectorSize); \
2410 for (size_t i=0; i<num_my_row; ++i) { \
2411 const GlobalOrdinal row = myGIDs[i]; \
2412 if (row == 0 || row == nrow-1) \
2413 b_view(i,0) = Scalar(0.0); \
2415 b_view(i,0) = -Scalar(b_val * h * h); \
2420 typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver; \
2421 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map); \
2422 std::string solver_name(#SolverName); \
2423 out << "Solving linear system with " << solver_name << std::endl; \
2424 RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>( \
2425 solver_name, matrix, x, b); \
2429 solver = Teuchos::null; \
2430 typedef Teuchos::ScalarTraits<BaseScalar> ST; \
2431 typename ST::magnitudeType tol = 1e-9; \
2432 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly); \
2433 Scalar val(VectorSize, BaseScalar(0.0)); \
2434 for (size_t i=0; i<num_my_row; ++i) { \
2435 const GlobalOrdinal row = myGIDs[i]; \
2436 BaseScalar xx = row * h; \
2437 for (LocalOrdinal j=0; j<VectorSize; ++j) { \
2438 val.fastAccessCoeff(j) = \
2439 BaseScalar(0.5) * (b_val.coeff(j)/a_val.coeff(j)) * xx * (xx - BaseScalar(1.0)); \
2441 TEST_EQUALITY( x_view(i,0).size(), VectorSize ); \
2444 Scalar v = x_view(i,0); \
2445 for (LocalOrdinal j=0; j<VectorSize; ++j) { \
2446 if (ST::magnitude(v.coeff(j)) < tol) \
2447 v.fastAccessCoeff(j) = BaseScalar(0.0); \
2448 if (ST::magnitude(val.coeff(j)) < tol) \
2449 val.fastAccessCoeff(j) = BaseScalar(0.0); \
2452 for (LocalOrdinal j=0; j<VectorSize; ++j) \
2453 TEST_FLOATING_EQUALITY(v.coeff(j), val.coeff(j), tol); \
2457 #endif // defined(HAVE_STOKHOS_AMESOS2)
2459 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_BASKER)
2460 TEST_AMESOS2_SOLVER(basker)
2465 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_KLU2)
2466 TEST_AMESOS2_SOLVER(klu2)
2471 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_SUPERLUDIST)
2472 TEST_AMESOS2_SOLVER(superlu_dist)
2477 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_SUPERLUMT)
2478 TEST_AMESOS2_SOLVER(superlu_mt)
2483 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_SUPERLU)
2484 TEST_AMESOS2_SOLVER(superlu)
2489 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_PARDISO_MKL)
2490 TEST_AMESOS2_SOLVER(pardisomkl)
2495 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_LAPACK)
2496 TEST_AMESOS2_SOLVER(lapack)
2501 #if defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
2502 TEST_AMESOS2_SOLVER(cholmod)
2507 #define CRSMATRIX_MP_VECTOR_TESTS_SLGN(S, LO, GO, N) \
2508 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorAdd, S, LO, GO, N ) \
2509 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorDot, S, LO, GO, N ) \
2510 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorAdd, S, LO, GO, N ) \
2511 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDot, S, LO, GO, N ) \
2512 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDotSub, S, LO, GO, N ) \
2513 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixVectorMultiply, S, LO, GO, N ) \
2514 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, S, LO, GO, N ) \
2515 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, WrappedDualView, S, LO, GO, N ) \
2516 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Flatten, S, LO, GO, N ) \
2517 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimpleCG, S, LO, GO, N ) \
2518 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimplePCG_Muelu, S, LO, GO, N ) \
2519 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES, S, LO, GO, N ) \
2520 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, S, LO, GO, N ) \
2521 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, S, LO, GO, N ) \
2522 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, S, LO, GO, N ) \
2523 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, S, LO, GO, N ) \
2524 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosCG_Muelu, S, LO, GO, N ) \
2525 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_basker, S, LO, GO, N ) \
2526 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_klu2, S, LO, GO, N ) \
2527 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_superlu_dist, S, LO, GO, N ) \
2528 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_superlu_mt, S, LO, GO, N ) \
2529 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_superlu, S, LO, GO, N ) \
2530 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_pardisomkl, S, LO, GO, N ) \
2531 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_lapack, S, LO, GO, N ) \
2532 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2_cholmod, S, LO, GO, N )
2534 #define CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N) \
2535 typedef Stokhos::DeviceForNode<N>::type Device; \
2536 typedef Stokhos::StaticFixedStorage<int,double,VectorSize,Device::execution_space> SFS; \
2537 using default_global_ordinal_type = ::Tpetra::Map<>::global_ordinal_type; \
2538 using default_local_ordinal_type = ::Tpetra::Map<>::local_ordinal_type; \
2539 CRSMATRIX_MP_VECTOR_TESTS_SLGN(SFS, default_local_ordinal_type, default_global_ordinal_type, N)
2541 #define CRSMATRIX_MP_VECTOR_TESTS_N(N) \
2542 CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Node > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &cijk_graph, const CijkType &cijk_dev)
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)
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
#define TEST_EQUALITY_CONST(v1, v2)
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
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_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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)
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > create_flat_mp_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > &graph, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_range_map, const LocalOrdinal block_size)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j)-expr2.val(j)
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
Tpetra::KokkosClassic::DefaultNode::DefaultNodeType Node
#define TEST_EQUALITY(v1, v2)
const unsigned VectorSize