42#include "Teuchos_UnitTestHelpers.hpp"
44#include "Stokhos_Ensemble_Sizes.hpp"
47#include "Teuchos_XMLParameterListCoreHelpers.hpp"
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"
58#include "Tpetra_Details_WrappedDualView.hpp"
62#ifdef HAVE_STOKHOS_BELOS
64#include "BelosLinearProblem.hpp"
65#include "BelosPseudoBlockGmresSolMgr.hpp"
66#include "BelosPseudoBlockCGSolMgr.hpp"
70#ifdef HAVE_STOKHOS_IFPACK2
72#include "Ifpack2_Factory.hpp"
76#ifdef HAVE_STOKHOS_MUELU
78#include "MueLu_CreateTpetraPreconditioner.hpp"
82#ifdef HAVE_STOKHOS_AMESOS2
84#include "Amesos2_Factory.hpp"
87template <
typename scalar,
typename ordinal>
91 const ordinal iColFEM,
92 const ordinal iStoch )
94 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
95 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
96 return X_fem + X_stoch;
100template <
typename scalar,
typename ordinal>
104 const ordinal nStoch,
105 const ordinal iColFEM,
107 const ordinal iStoch)
109 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
110 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
111 const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
112 return X_fem + X_stoch + X_col;
130 using Teuchos::ArrayView;
131 using Teuchos::Array;
132 using Teuchos::ArrayRCP;
134 typedef typename Storage::value_type BaseScalar;
137 typedef Teuchos::Comm<int> Tpetra_Comm;
138 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
139 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
142 if ( !Kokkos::is_initialized() )
143 Kokkos::initialize();
146 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
150 RCP<const Tpetra_Map> map =
151 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
153 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
154 const size_t num_my_row = myGIDs.size();
157 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
158 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
160 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
161 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
163 for (
size_t i=0; i<num_my_row; ++i) {
166 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
167 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
177 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
178 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
184 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
186 BaseScalar tol = 1.0e-14;
187 for (
size_t i=0; i<num_my_row; ++i) {
190 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
192 val.fastAccessCoeff(
j) = alpha.coeff(
j)*v + 0.12345*beta.coeff(
j)*v;
194 TEST_EQUALITY( y_view(i,0).size(),
VectorSize );
208 using Teuchos::ArrayView;
209 using Teuchos::Array;
210 using Teuchos::ArrayRCP;
212 typedef typename Storage::value_type BaseScalar;
215 typedef Teuchos::Comm<int> Tpetra_Comm;
216 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
217 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
218 typedef typename Tpetra_Vector::dot_type dot_type;
221 if ( !Kokkos::is_initialized() )
222 Kokkos::initialize();
225 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
229 RCP<const Tpetra_Map> map =
230 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
232 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
233 const size_t num_my_row = myGIDs.size();
236 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
237 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
239 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
240 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
242 for (
size_t i=0; i<num_my_row; ++i) {
245 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
246 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
254 dot_type dot = x1->dot(*x2);
258#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
261 dot_type local_val(0);
262 for (
size_t i=0; i<num_my_row; ++i) {
265 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
267 local_val += 0.12345 * v * v;
273 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
274 Teuchos::outArg(
val));
280 for (
size_t i=0; i<num_my_row; ++i) {
283 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
285 local_val.fastAccessCoeff(
j) += 0.12345 * v * v;
291 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
292 Teuchos::outArg(
val));
296 out <<
"dot = " << dot <<
" expected = " <<
val << std::endl;
298 BaseScalar tol = 1.0e-14;
299 TEST_FLOATING_EQUALITY( dot,
val, tol );
310 using Teuchos::ArrayView;
311 using Teuchos::Array;
312 using Teuchos::ArrayRCP;
314 typedef typename Storage::value_type BaseScalar;
317 typedef Teuchos::Comm<int> Tpetra_Comm;
318 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
319 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
322 if ( !Kokkos::is_initialized() )
323 Kokkos::initialize();
326 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
330 RCP<const Tpetra_Map> map =
331 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
333 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
334 const size_t num_my_row = myGIDs.size();
338 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
339 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
341 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
342 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
344 for (
size_t i=0; i<num_my_row; ++i) {
346 for (
size_t j=0;
j<ncol; ++
j) {
349 generate_multi_vector_coefficient<BaseScalar,size_t>(
351 val1.fastAccessCoeff(k) = v;
352 val2.fastAccessCoeff(k) = 0.12345 * v;
363 RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
364 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
370 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
372 BaseScalar tol = 1.0e-14;
373 for (
size_t i=0; i<num_my_row; ++i) {
375 for (
size_t j=0;
j<ncol; ++
j) {
377 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
379 val.fastAccessCoeff(k) = alpha.coeff(k)*v + 0.12345*beta.coeff(k)*v;
384 val.fastAccessCoeff(k), tol );
397 using Teuchos::ArrayView;
398 using Teuchos::Array;
399 using Teuchos::ArrayRCP;
401 typedef typename Storage::value_type BaseScalar;
404 typedef Teuchos::Comm<int> Tpetra_Comm;
405 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
406 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
407 typedef typename Tpetra_MultiVector::dot_type dot_type;
410 if ( !Kokkos::is_initialized() )
411 Kokkos::initialize();
414 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
418 RCP<const Tpetra_Map> map =
419 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
421 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
422 const size_t num_my_row = myGIDs.size();
426 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
427 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
429 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
430 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
432 for (
size_t i=0; i<num_my_row; ++i) {
434 for (
size_t j=0;
j<ncol; ++
j) {
437 generate_multi_vector_coefficient<BaseScalar,size_t>(
439 val1.fastAccessCoeff(k) = v;
440 val2.fastAccessCoeff(k) = 0.12345 * v;
449 Array<dot_type> dots(ncol);
450 x1->dot(*x2, dots());
454#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
457 Array<dot_type> local_vals(ncol, dot_type(0));
458 for (
size_t i=0; i<num_my_row; ++i) {
460 for (
size_t j=0;
j<ncol; ++
j) {
462 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
464 local_vals[
j] += 0.12345 * v * v;
470 Array<dot_type> vals(ncol, dot_type(0));
471 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
472 local_vals.getRawPtr(), vals.getRawPtr());
477 Array<dot_type> local_vals(ncol, dot_type(
VectorSize, 0.0));
478 for (
size_t i=0; i<num_my_row; ++i) {
480 for (
size_t j=0;
j<ncol; ++
j) {
482 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
484 local_vals[
j].fastAccessCoeff(k) += 0.12345 * v * v;
490 Array<dot_type> vals(ncol, dot_type(
VectorSize, 0.0));
491 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
492 local_vals.getRawPtr(), vals.getRawPtr());
496 BaseScalar tol = 1.0e-14;
497 for (
size_t j=0;
j<ncol; ++
j) {
498 out <<
"dots(" <<
j <<
") = " << dots[
j]
499 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
500 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
512 using Teuchos::ArrayView;
513 using Teuchos::Array;
514 using Teuchos::ArrayRCP;
516 typedef typename Storage::value_type BaseScalar;
519 typedef Teuchos::Comm<int> Tpetra_Comm;
520 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
521 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
522 typedef typename Tpetra_MultiVector::dot_type dot_type;
525 if ( !Kokkos::is_initialized() )
526 Kokkos::initialize();
529 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
533 RCP<const Tpetra_Map> map =
534 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
536 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
537 const size_t num_my_row = myGIDs.size();
541 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
542 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
544 auto x1_view = x1->getLocalViewHost(Tpetra::Access::OverwriteAll);
545 auto x2_view = x2->getLocalViewHost(Tpetra::Access::OverwriteAll);
547 for (
size_t i=0; i<num_my_row; ++i) {
549 for (
size_t j=0;
j<ncol; ++
j) {
552 generate_multi_vector_coefficient<BaseScalar,size_t>(
554 val1.fastAccessCoeff(k) = v;
555 val2.fastAccessCoeff(k) = 0.12345 * v;
565 Teuchos::Array<size_t> cols(ncol_sub);
566 cols[0] = 4; cols[1] = 2;
567 RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
568 RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
571 Array<dot_type> dots(ncol_sub);
572 x1_sub->dot(*x2_sub, dots());
576#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
579 Array<dot_type> local_vals(ncol_sub, dot_type(0));
580 for (
size_t i=0; i<num_my_row; ++i) {
582 for (
size_t j=0;
j<ncol_sub; ++
j) {
584 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
586 local_vals[
j] += 0.12345 * v * v;
592 Array<dot_type> vals(ncol_sub, dot_type(0));
593 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
594 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
600 Array<dot_type> local_vals(ncol_sub, dot_type(
VectorSize, 0.0));
601 for (
size_t i=0; i<num_my_row; ++i) {
603 for (
size_t j=0;
j<ncol_sub; ++
j) {
605 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
607 local_vals[
j].fastAccessCoeff(k) += 0.12345 * v * v;
613 Array<dot_type> vals(ncol_sub, dot_type(
VectorSize, 0.0));
614 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
615 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
620 BaseScalar tol = 1.0e-14;
621 for (
size_t j=0;
j<ncol_sub; ++
j) {
622 out <<
"dots(" <<
j <<
") = " << dots[
j]
623 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
624 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
636 using Teuchos::ArrayView;
637 using Teuchos::Array;
638 using Teuchos::ArrayRCP;
640 typedef typename Storage::value_type BaseScalar;
643 typedef Teuchos::Comm<int> Tpetra_Comm;
644 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
645 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
646 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
647 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
650 if ( !Kokkos::is_initialized() )
651 Kokkos::initialize();
655 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
656 RCP<const Tpetra_Map> map =
657 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
659 RCP<Tpetra_CrsGraph> graph =
660 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
661 Array<GlobalOrdinal> columnIndices(2);
662 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
663 const size_t num_my_row = myGIDs.size();
664 for (
size_t i=0; i<num_my_row; ++i) {
666 columnIndices[0] = row;
669 columnIndices[1] = row+1;
672 graph->insertGlobalIndices(row, columnIndices(0,ncol));
674 graph->fillComplete();
675 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
678 Array<Scalar> vals(2);
680 for (
size_t i=0; i<num_my_row; ++i) {
682 columnIndices[0] = row;
684 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
690 columnIndices[1] = row+1;
692 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
697 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
699 matrix->fillComplete();
702 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
704 auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
705 for (
size_t i=0; i<num_my_row; ++i) {
708 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
721 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
722 matrix->apply(*x, *y);
728 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
729 BaseScalar tol = 1.0e-14;
730 for (
size_t i=0; i<num_my_row; ++i) {
733 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
735 val.fastAccessCoeff(
j) = v*v;
739 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
741 val.fastAccessCoeff(
j) += v*v;
744 TEST_EQUALITY( y_view(i,0).size(),
VectorSize );
758 using Teuchos::ArrayView;
759 using Teuchos::Array;
760 using Teuchos::ArrayRCP;
762 typedef typename Storage::value_type BaseScalar;
765 typedef Teuchos::Comm<int> Tpetra_Comm;
766 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
767 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
768 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
769 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
772 if ( !Kokkos::is_initialized() )
773 Kokkos::initialize();
777 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
778 RCP<const Tpetra_Map> map =
779 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
781 RCP<Tpetra_CrsGraph> graph =
782 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
783 Array<GlobalOrdinal> columnIndices(2);
784 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
785 const size_t num_my_row = myGIDs.size();
786 for (
size_t i=0; i<num_my_row; ++i) {
788 columnIndices[0] = row;
791 columnIndices[1] = row+1;
794 graph->insertGlobalIndices(row, columnIndices(0,ncol));
796 graph->fillComplete();
797 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
800 Array<Scalar> vals(2);
802 for (
size_t i=0; i<num_my_row; ++i) {
804 columnIndices[0] = row;
806 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
812 columnIndices[1] = row+1;
814 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
819 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
821 matrix->fillComplete();
825 RCP<Tpetra_MultiVector> x = Tpetra::createMultiVector<Scalar>(map, ncol);
827 auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
828 for (
size_t i=0; i<num_my_row; ++i) {
830 for (
size_t j=0;
j<ncol; ++
j) {
833 generate_multi_vector_coefficient<BaseScalar,size_t>(
835 val.fastAccessCoeff(k) = v;
849 RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
850 matrix->apply(*x, *y);
856 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
857 BaseScalar tol = 1.0e-14;
858 for (
size_t i=0; i<num_my_row; ++i) {
860 for (
size_t j=0;
j<ncol; ++
j) {
862 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
864 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
866 val.fastAccessCoeff(k) = v1*v2;
870 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
872 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
874 val.fastAccessCoeff(k) += v1*v2;
880 val.fastAccessCoeff(k), tol );
893 using Teuchos::ArrayView;
894 using Teuchos::Array;
895 using Teuchos::ArrayRCP;
897 typedef typename Storage::value_type BaseScalar;
900 typedef Teuchos::Comm<int> Tpetra_Comm;
901 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
902 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
903 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
904 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
906 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
907 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
910 if ( !Kokkos::is_initialized() )
911 Kokkos::initialize();
915 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
916 RCP<const Tpetra_Map> map =
917 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
919 RCP<Tpetra_CrsGraph> graph =
920 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
921 Array<GlobalOrdinal> columnIndices(2);
922 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
923 const size_t num_my_row = myGIDs.size();
924 for (
size_t i=0; i<num_my_row; ++i) {
926 columnIndices[0] = row;
929 columnIndices[1] = row+1;
932 graph->insertGlobalIndices(row, columnIndices(0,ncol));
934 graph->fillComplete();
935 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
938 Array<Scalar> vals(2);
940 for (
size_t i=0; i<num_my_row; ++i) {
942 columnIndices[0] = row;
944 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
950 columnIndices[1] = row+1;
952 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
957 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
959 matrix->fillComplete();
962 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
964 auto x_view = x->getLocalViewHost(Tpetra::Access::OverwriteAll);
965 for (
size_t i=0; i<num_my_row; ++i) {
968 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
975 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
976 matrix->apply(*x, *y);
979 RCP<const Tpetra_Map> flat_x_map, flat_y_map;
980 RCP<const Tpetra_CrsGraph> flat_graph =
982 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
986 RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
988 RCP<Flat_Tpetra_Vector> flat_x =
990 RCP<Flat_Tpetra_Vector> flat_y =
992 flat_matrix->apply(*flat_x, *flat_y);
999 BaseScalar tol = 1.0e-14;
1000 auto y_view = y-> getLocalViewHost(Tpetra::Access::ReadOnly);
1001 auto y2_view = y2->getLocalViewHost(Tpetra::Access::ReadOnly);
1002 for (
size_t i=0; i<num_my_row; ++i) {
1003 TEST_EQUALITY( y_view( i,0).size(),
VectorSize );
1004 TEST_EQUALITY( y2_view(i,0).size(),
VectorSize );
1024 using Teuchos::ArrayView;
1025 using Teuchos::Array;
1026 using Teuchos::ArrayRCP;
1031 using DualViewType = Kokkos::DualView<Scalar*, typename Node::device_type>;
1032 using WDV = Tpetra::Details::WrappedDualView<DualViewType>;
1033 using values_view =
typename DualViewType::t_dev;
1036 if ( !Kokkos::is_initialized() )
1037 Kokkos::initialize();
1041 values_view myView(
"emptyTestView", 0);
1044 size_t use_h = wdv.getHostView(Tpetra::Access::ReadOnly).use_count();
1045 size_t use_d = wdv.getDeviceView(Tpetra::Access::ReadOnly).use_count();
1048 TEST_EQUALITY(use_h, use_d);
1059 using Teuchos::ArrayView;
1060 using Teuchos::Array;
1061 using Teuchos::ArrayRCP;
1062 using Teuchos::ParameterList;
1064 typedef typename Storage::value_type BaseScalar;
1067 typedef Teuchos::Comm<int> Tpetra_Comm;
1068 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1069 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1070 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1071 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1074 if ( !Kokkos::is_initialized() )
1075 Kokkos::initialize();
1079 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1080 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1081 RCP<const Tpetra_Map> map =
1082 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1084 RCP<Tpetra_CrsGraph> graph =
1085 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
1086 Array<GlobalOrdinal> columnIndices(3);
1087 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1088 const size_t num_my_row = myGIDs.size();
1089 for (
size_t i=0; i<num_my_row; ++i) {
1091 if (row == 0 || row == nrow-1) {
1092 columnIndices[0] = row;
1093 graph->insertGlobalIndices(row, columnIndices(0,1));
1096 columnIndices[0] = row-1;
1097 columnIndices[1] = row;
1098 columnIndices[2] = row+1;
1099 graph->insertGlobalIndices(row, columnIndices(0,3));
1102 graph->fillComplete();
1103 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1106 Array<Scalar> vals(3);
1109 a_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 columnIndices[0] = row;
1117 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1120 columnIndices[0] = row-1;
1121 columnIndices[1] = row;
1122 columnIndices[2] = row+1;
1123 vals[0] =
Scalar(-1.0) * a_val;
1124 vals[1] =
Scalar(2.0) * a_val;
1125 vals[2] =
Scalar(-1.0) * a_val;
1126 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1129 matrix->fillComplete();
1131 matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
1132 Teuchos::VERB_EXTREME);
1135 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1138 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1141 b_val.fastAccessCoeff(
j) =
1142 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1144 for (
size_t i=0; i<num_my_row; ++i) {
1146 if (row == 0 || row == nrow-1)
1147 b_view(i,0) =
Scalar(0.0);
1149 b_view(i,0) = -
Scalar(b_val * h * h);
1154 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1155 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1156 typedef typename BST::mag_type base_mag_type;
1157 typedef typename Tpetra_Vector::mag_type mag_type;
1158 base_mag_type btol = 1e-9;
1159 mag_type tol = btol;
1162 out.getOStream().get());
1163 TEST_EQUALITY_CONST( solved,
true );
1170 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
1172 for (
size_t i=0; i<num_my_row; ++i) {
1174 BaseScalar xx = row * h;
1176 val.fastAccessCoeff(
j) =
1177 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1179 TEST_EQUALITY( x_view(i,0).size(),
VectorSize );
1184 if (BST::abs(v.coeff(
j)) < btol)
1185 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1186 if (BST::abs(
val.coeff(
j)) < btol)
1187 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1191 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), btol);
1196#if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2)
1206 using Teuchos::ArrayView;
1207 using Teuchos::Array;
1208 using Teuchos::ArrayRCP;
1209 using Teuchos::ParameterList;
1210 using Teuchos::getParametersFromXmlFile;
1212 typedef typename Storage::value_type BaseScalar;
1215 typedef Teuchos::Comm<int> Tpetra_Comm;
1216 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1217 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1218 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1219 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1222 if ( !Kokkos::is_initialized() )
1223 Kokkos::initialize();
1227 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1228 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1229 RCP<const Tpetra_Map> map =
1230 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1232 RCP<Tpetra_CrsGraph> graph =
1233 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
1234 Array<GlobalOrdinal> columnIndices(3);
1235 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1236 const size_t num_my_row = myGIDs.size();
1237 for (
size_t i=0; i<num_my_row; ++i) {
1239 if (row == 0 || row == nrow-1) {
1240 columnIndices[0] = row;
1241 graph->insertGlobalIndices(row, columnIndices(0,1));
1244 columnIndices[0] = row-1;
1245 columnIndices[1] = row;
1246 columnIndices[2] = row+1;
1247 graph->insertGlobalIndices(row, columnIndices(0,3));
1250 graph->fillComplete();
1251 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1254 Array<Scalar> vals(3);
1257 a_val.fastAccessCoeff(
j) =
1258 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1260 for (
size_t i=0; i<num_my_row; ++i) {
1262 if (row == 0 || row == nrow-1) {
1263 columnIndices[0] = row;
1265 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1268 columnIndices[0] = row-1;
1269 columnIndices[1] = row;
1270 columnIndices[2] = row+1;
1271 vals[0] =
Scalar(-1.0) * a_val;
1272 vals[1] =
Scalar(2.0) * a_val;
1273 vals[2] =
Scalar(-1.0) * a_val;
1274 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1277 matrix->fillComplete();
1280 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1283 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1286 b_val.fastAccessCoeff(
j) =
1287 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1289 for (
size_t i=0; i<num_my_row; ++i) {
1291 if (row == 0 || row == nrow-1)
1292 b_view(i,0) =
Scalar(0.0);
1294 b_view(i,0) = -
Scalar(b_val * h * h);
1299 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1300 RCP<OP> matrix_op = matrix;
1301 RCP<ParameterList> muelu_params =
1302 getParametersFromXmlFile(
"muelu_cheby.xml");
1304 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1307 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1308 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1309 typedef typename BST::mag_type base_mag_type;
1310 typedef typename Tpetra_Vector::mag_type mag_type;
1311 base_mag_type btol = 1e-9;
1312 mag_type tol = btol;
1315 out.getOStream().get());
1316 TEST_EQUALITY_CONST( solved,
true );
1323 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
1325 for (
size_t i=0; i<num_my_row; ++i) {
1327 BaseScalar xx = row * h;
1329 val.fastAccessCoeff(
j) =
1330 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1332 TEST_EQUALITY( x_view(i,0).size(),
VectorSize );
1337 if (BST::magnitude(v.coeff(
j)) < btol)
1338 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1339 if (BST::magnitude(
val.coeff(
j)) < btol)
1340 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1344 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), btol);
1356#if defined(HAVE_STOKHOS_BELOS)
1366 using Teuchos::ArrayView;
1367 using Teuchos::Array;
1368 using Teuchos::ArrayRCP;
1369 using Teuchos::ParameterList;
1371 typedef typename Storage::value_type BaseScalar;
1374 typedef Teuchos::Comm<int> Tpetra_Comm;
1375 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1376 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1377 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1378 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1381 if ( !Kokkos::is_initialized() )
1382 Kokkos::initialize();
1386 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1387 RCP<const Tpetra_Map> map =
1388 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1390 RCP<Tpetra_CrsGraph> graph =
1391 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
1392 Array<GlobalOrdinal> columnIndices(2);
1393 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1394 const size_t num_my_row = myGIDs.size();
1395 for (
size_t i=0; i<num_my_row; ++i) {
1397 columnIndices[0] = row;
1399 if (row != nrow-1) {
1400 columnIndices[1] = row+1;
1403 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1405 graph->fillComplete();
1406 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1409 Array<Scalar> vals(2);
1411 for (
size_t i=0; i<num_my_row; ++i) {
1413 columnIndices[0] = row;
1415 val.fastAccessCoeff(
j) =
j+1;
1419 if (row != nrow-1) {
1420 columnIndices[1] = row+1;
1422 val.fastAccessCoeff(
j) =
j+1;
1426 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1428 matrix->fillComplete();
1431 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1433 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1434 for (
size_t i=0; i<num_my_row; ++i) {
1435 b_view(i,0) =
Scalar(1.0);
1440 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1441#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1442 typedef BaseScalar BelosScalar;
1444 typedef Scalar BelosScalar;
1446 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1447 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1448 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1449 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1450 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
1451 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1452 typename ST::magnitudeType tol = 1e-9;
1453 belosParams->set(
"Flexible Gmres",
false);
1454 belosParams->set(
"Num Blocks", 100);
1455 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1456 belosParams->set(
"Maximum Iterations", 100);
1457 belosParams->set(
"Verbosity", 33);
1458 belosParams->set(
"Output Style", 1);
1459 belosParams->set(
"Output Frequency", 1);
1460 belosParams->set(
"Output Stream", out.getOStream());
1461 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1462 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1463 problem->setProblem();
1464 Belos::ReturnType ret = solver->solve();
1465 TEST_EQUALITY_CONST( ret, Belos::Converged );
1477 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
1478 for (
size_t i=0; i<num_my_row; ++i) {
1482 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
1487 TEST_EQUALITY( x_view(i,0).size(),
VectorSize );
1492 if (ST::magnitude(v.coeff(
j)) < tol)
1493 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1497 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1509 using Teuchos::tuple;
1510 using Teuchos::ArrayView;
1511 using Teuchos::Array;
1512 using Teuchos::ArrayRCP;
1513 using Teuchos::ParameterList;
1515 typedef typename Storage::value_type BaseScalar;
1518 typedef Teuchos::Comm<int> Tpetra_Comm;
1519 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1520 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1521 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1522 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1525 if ( !Kokkos::is_initialized() )
1526 Kokkos::initialize();
1530 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1531 RCP<const Tpetra_Map> map =
1532 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1534 RCP<Tpetra_CrsGraph> graph =
1535 rcp(
new Tpetra_CrsGraph(map,
size_t(1)));
1536 Array<GlobalOrdinal> columnIndices(1);
1537 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1538 const size_t num_my_row = myGIDs.size();
1539 for (
size_t i=0; i<num_my_row; ++i) {
1541 columnIndices[0] = row;
1543 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1545 graph->fillComplete();
1546 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1548 Array<Scalar> vals(1);
1549 for (
size_t i=0; i<num_my_row; ++i) {
1551 columnIndices[0] = row;
1553 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1555 matrix->fillComplete();
1564 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1566 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1567 for (
size_t i=0; i<num_my_row; ++i) {
1569 b_view(i,0) =
Scalar(0.0);
1572 b_view(i,0).fastAccessCoeff(
j) = BaseScalar(row+1);
1577 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1578#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1579 typedef BaseScalar BelosScalar;
1581 typedef Scalar BelosScalar;
1583 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1584 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1585 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1586 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1587 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
1588 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1589 typename ST::magnitudeType tol = 1e-9;
1590 belosParams->set(
"Flexible Gmres",
false);
1591 belosParams->set(
"Num Blocks", 100);
1592 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1593 belosParams->set(
"Maximum Iterations", 100);
1594 belosParams->set(
"Verbosity", 33);
1595 belosParams->set(
"Output Style", 1);
1596 belosParams->set(
"Output Frequency", 1);
1597 belosParams->set(
"Output Stream", out.getOStream());
1598 belosParams->set(
"Orthogonalization",
"DGKS");
1599 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1600 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1601 problem->setProblem();
1602 Belos::ReturnType ret = solver->solve();
1603 TEST_EQUALITY_CONST( ret, Belos::Converged );
1605#ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1606 int numItersExpected = nrow;
1607 int numIters = solver->getNumIters();
1608 out <<
"numIters = " << numIters << std::endl;
1609 TEST_EQUALITY( numIters, numItersExpected);
1612 std::vector<int> ensemble_iterations =
1614 out <<
"Ensemble iterations = ";
1615 for (
auto ensemble_iteration : ensemble_iterations)
1616 out << ensemble_iteration <<
" ";
1621 TEST_EQUALITY(
int(
j+1+nrow-
VectorSize), ensemble_iterations[
j]);
1624 TEST_EQUALITY(
int(0), ensemble_iterations[
j]);
1636 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
1637 for (
size_t i=0; i<num_my_row; ++i) {
1642 if (ST::magnitude(v.coeff(
j)) < tol)
1643 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1650 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1653 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1665 using Teuchos::tuple;
1666 using Teuchos::ArrayView;
1667 using Teuchos::Array;
1668 using Teuchos::ArrayRCP;
1669 using Teuchos::ParameterList;
1671 typedef typename Storage::value_type BaseScalar;
1674 typedef Teuchos::Comm<int> Tpetra_Comm;
1675 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1676 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1677 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1678 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1681 if ( !Kokkos::is_initialized() )
1682 Kokkos::initialize();
1686 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1687 RCP<const Tpetra_Map> map =
1688 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1690 RCP<Tpetra_CrsGraph> graph =
1691 rcp(
new Tpetra_CrsGraph(map,
size_t(1)));
1692 Array<GlobalOrdinal> columnIndices(1);
1693 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1694 const size_t num_my_row = myGIDs.size();
1695 for (
size_t i=0; i<num_my_row; ++i) {
1697 columnIndices[0] = row;
1699 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1701 graph->fillComplete();
1702 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1704 Array<Scalar> vals(1);
1705 for (
size_t i=0; i<num_my_row; ++i) {
1707 columnIndices[0] = row;
1709 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1711 matrix->fillComplete();
1720 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1722 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1723 for (
size_t i=0; i<num_my_row; ++i) {
1725 b_view(i,0) =
Scalar(0.0);
1728 b_view(i,0).fastAccessCoeff(
j) = BaseScalar(row+1);
1733 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1734#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1735 typedef BaseScalar BelosScalar;
1737 typedef Scalar BelosScalar;
1739 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1740 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1741 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1742 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1743 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
1744 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1745 typename ST::magnitudeType tol = 1e-9;
1746 belosParams->set(
"Flexible Gmres",
false);
1747 belosParams->set(
"Num Blocks", 100);
1748 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1749 belosParams->set(
"Maximum Iterations", 100);
1750 belosParams->set(
"Verbosity", 33);
1751 belosParams->set(
"Output Style", 1);
1752 belosParams->set(
"Output Frequency", 1);
1753 belosParams->set(
"Output Stream", out.getOStream());
1754 belosParams->set(
"Orthogonalization",
"ICGS");
1755 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1756 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1757 problem->setProblem();
1758 Belos::ReturnType ret = solver->solve();
1759 TEST_EQUALITY_CONST( ret, Belos::Converged );
1761#ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1762 int numItersExpected = nrow;
1763 int numIters = solver->getNumIters();
1764 out <<
"numIters = " << numIters << std::endl;
1765 TEST_EQUALITY( numIters, numItersExpected);
1768 std::vector<int> ensemble_iterations =
1770 out <<
"Ensemble iterations = ";
1771 for (
auto ensemble_iteration : ensemble_iterations)
1772 out << ensemble_iteration <<
" ";
1777 TEST_EQUALITY(
int(
j+1+nrow-
VectorSize), ensemble_iterations[
j]);
1780 TEST_EQUALITY(
int(0), ensemble_iterations[
j]);
1792 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
1793 for (
size_t i=0; i<num_my_row; ++i) {
1798 if (ST::magnitude(v.coeff(
j)) < tol)
1799 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1806 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1809 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1821 using Teuchos::tuple;
1822 using Teuchos::ArrayView;
1823 using Teuchos::Array;
1824 using Teuchos::ArrayRCP;
1825 using Teuchos::ParameterList;
1827 typedef typename Storage::value_type BaseScalar;
1830 typedef Teuchos::Comm<int> Tpetra_Comm;
1831 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1832 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1833 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1834 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1837 if ( !Kokkos::is_initialized() )
1838 Kokkos::initialize();
1842 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1843 RCP<const Tpetra_Map> map =
1844 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1846 RCP<Tpetra_CrsGraph> graph =
1847 rcp(
new Tpetra_CrsGraph(map,
size_t(1)));
1848 Array<GlobalOrdinal> columnIndices(1);
1849 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1850 const size_t num_my_row = myGIDs.size();
1851 for (
size_t i=0; i<num_my_row; ++i) {
1853 columnIndices[0] = row;
1855 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1857 graph->fillComplete();
1858 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1860 Array<Scalar> vals(1);
1861 for (
size_t i=0; i<num_my_row; ++i) {
1863 columnIndices[0] = row;
1865 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1867 matrix->fillComplete();
1876 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1878 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
1879 for (
size_t i=0; i<num_my_row; ++i) {
1881 b_view(i,0) =
Scalar(0.0);
1884 b_view(i,0).fastAccessCoeff(
j) = BaseScalar(row+1);
1889 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1890#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
1891 typedef BaseScalar BelosScalar;
1893 typedef Scalar BelosScalar;
1895 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1896 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1897 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1898 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1899 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
1900 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1901 typename ST::magnitudeType tol = 1e-9;
1902 belosParams->set(
"Flexible Gmres",
false);
1903 belosParams->set(
"Num Blocks", 100);
1904 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1905 belosParams->set(
"Maximum Iterations", 100);
1906 belosParams->set(
"Verbosity", 33);
1907 belosParams->set(
"Output Style", 1);
1908 belosParams->set(
"Output Frequency", 1);
1909 belosParams->set(
"Output Stream", out.getOStream());
1910 belosParams->set(
"Orthogonalization",
"IMGS");
1911 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1912 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1913 problem->setProblem();
1914 Belos::ReturnType ret = solver->solve();
1915 TEST_EQUALITY_CONST( ret, Belos::Converged );
1917#ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
1918 int numItersExpected = nrow;
1919 int numIters = solver->getNumIters();
1920 out <<
"numIters = " << numIters << std::endl;
1921 TEST_EQUALITY( numIters, numItersExpected);
1924 std::vector<int> ensemble_iterations =
1926 out <<
"Ensemble iterations = ";
1927 for (
auto ensemble_iteration : ensemble_iterations)
1928 out << ensemble_iteration <<
" ";
1933 TEST_EQUALITY(
int(
j+1+nrow-
VectorSize), ensemble_iterations[
j]);
1936 TEST_EQUALITY(
int(0), ensemble_iterations[
j]);
1948 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
1949 for (
size_t i=0; i<num_my_row; ++i) {
1954 if (ST::magnitude(v.coeff(
j)) < tol)
1955 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1962 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1965 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1989#if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2)
2000 using Teuchos::ArrayView;
2001 using Teuchos::Array;
2002 using Teuchos::ArrayRCP;
2003 using Teuchos::ParameterList;
2005 typedef typename Storage::value_type BaseScalar;
2008 typedef Teuchos::Comm<int> Tpetra_Comm;
2009 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2010 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2011 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2012 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2015 if ( !Kokkos::is_initialized() )
2016 Kokkos::initialize();
2020 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2021 RCP<const Tpetra_Map> map =
2022 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2024 RCP<Tpetra_CrsGraph> graph =
2025 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
2026 Array<GlobalOrdinal> columnIndices(2);
2027 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2028 const size_t num_my_row = myGIDs.size();
2029 for (
size_t i=0; i<num_my_row; ++i) {
2031 columnIndices[0] = row;
2033 if (row != nrow-1) {
2034 columnIndices[1] = row+1;
2037 graph->insertGlobalIndices(row, columnIndices(0,ncol));
2039 graph->fillComplete();
2040 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
2043 Array<Scalar> vals(2);
2045 for (
size_t i=0; i<num_my_row; ++i) {
2047 columnIndices[0] = row;
2049 val.fastAccessCoeff(
j) =
j+1;
2053 if (row != nrow-1) {
2054 columnIndices[1] = row+1;
2056 val.fastAccessCoeff(
j) =
j+1;
2060 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
2062 matrix->fillComplete();
2065 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2067 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
2068 for (
size_t i=0; i<num_my_row; ++i) {
2069 b_view(i,0) =
Scalar(1.0);
2074 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
2075 Ifpack2::Factory factory;
2076 RCP<Prec> M = factory.create<Tpetra_CrsMatrix>(
"RILUK", matrix);
2081 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2082#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
2083 typedef BaseScalar BelosScalar;
2085 typedef Scalar BelosScalar;
2087 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2088 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2089 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2090 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2091 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
2092 problem->setRightPrec(M);
2093 problem->setProblem();
2094 RCP<ParameterList> belosParams = rcp(
new ParameterList);
2095 typename ST::magnitudeType tol = 1e-9;
2096 belosParams->set(
"Flexible Gmres",
false);
2097 belosParams->set(
"Num Blocks", 100);
2098 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
2099 belosParams->set(
"Maximum Iterations", 100);
2100 belosParams->set(
"Verbosity", 33);
2101 belosParams->set(
"Output Style", 1);
2102 belosParams->set(
"Output Frequency", 1);
2103 belosParams->set(
"Output Stream", out.getOStream());
2105 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2106 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2107 Belos::ReturnType ret = solver->solve();
2108 TEST_EQUALITY_CONST( ret, Belos::Converged );
2120 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
2121 for (
size_t i=0; i<num_my_row; ++i) {
2125 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
2130 TEST_EQUALITY( x_view(i,0).size(),
VectorSize );
2135 if (ST::magnitude(v.coeff(
j)) < tol)
2136 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2140 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
2152#if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2)
2162 using Teuchos::ArrayView;
2163 using Teuchos::Array;
2164 using Teuchos::ArrayRCP;
2165 using Teuchos::ParameterList;
2166 using Teuchos::getParametersFromXmlFile;
2168 typedef typename Storage::value_type BaseScalar;
2171 typedef Teuchos::Comm<int> Tpetra_Comm;
2172 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2173 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2174 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2175 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2178 if ( !Kokkos::is_initialized() )
2179 Kokkos::initialize();
2183 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
2184 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2185 RCP<const Tpetra_Map> map =
2186 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2188 RCP<Tpetra_CrsGraph> graph =
2189 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
2190 Array<GlobalOrdinal> columnIndices(3);
2191 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2192 const size_t num_my_row = myGIDs.size();
2193 for (
size_t i=0; i<num_my_row; ++i) {
2195 if (row == 0 || row == nrow-1) {
2196 columnIndices[0] = row;
2197 graph->insertGlobalIndices(row, columnIndices(0,1));
2200 columnIndices[0] = row-1;
2201 columnIndices[1] = row;
2202 columnIndices[2] = row+1;
2203 graph->insertGlobalIndices(row, columnIndices(0,3));
2206 graph->fillComplete();
2207 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
2210 Array<Scalar> vals(3);
2213 a_val.fastAccessCoeff(
j) =
2214 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
2216 for (
size_t i=0; i<num_my_row; ++i) {
2218 if (row == 0 || row == nrow-1) {
2219 columnIndices[0] = row;
2221 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2224 columnIndices[0] = row-1;
2225 columnIndices[1] = row;
2226 columnIndices[2] = row+1;
2227 vals[0] =
Scalar(-1.0) * a_val;
2228 vals[1] =
Scalar(2.0) * a_val;
2229 vals[2] =
Scalar(-1.0) * a_val;
2230 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2233 matrix->fillComplete();
2236 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2239 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
2242 b_val.fastAccessCoeff(
j) =
2243 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
2245 for (
size_t i=0; i<num_my_row; ++i) {
2247 if (row == 0 || row == nrow-1)
2248 b_view(i,0) =
Scalar(0.0);
2250 b_view(i,0) = -
Scalar(b_val * h * h);
2255 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2256 RCP<ParameterList> muelu_params =
2257 getParametersFromXmlFile(
"muelu_cheby.xml");
2258 RCP<OP> matrix_op = matrix;
2260 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
2263 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2264#ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT
2265 typedef BaseScalar BelosScalar;
2267 typedef Scalar BelosScalar;
2269 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2270 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2271 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2272 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
2273 problem->setRightPrec(M);
2274 problem->setProblem();
2275 RCP<ParameterList> belosParams = rcp(
new ParameterList);
2276 typename ST::magnitudeType tol = 1e-9;
2277 belosParams->set(
"Num Blocks", 100);
2278 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
2279 belosParams->set(
"Maximum Iterations", 100);
2280 belosParams->set(
"Verbosity", 33);
2281 belosParams->set(
"Output Style", 1);
2282 belosParams->set(
"Output Frequency", 1);
2283 belosParams->set(
"Output Stream", out.getOStream());
2286 belosParams->set(
"Implicit Residual Scaling",
"None");
2288 RCP<Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true> > solver =
2289 rcp(
new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true>(problem, belosParams));
2290 Belos::ReturnType ret = solver->solve();
2291 TEST_EQUALITY_CONST( ret, Belos::Converged );
2293#ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT
2295 std::vector<int> ensemble_iterations =
2296 solver->getResidualStatusTest()->getEnsembleIterations();
2297 out <<
"Ensemble iterations = ";
2299 out << ensemble_iterations[i] <<
" ";
2308 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
2310 for (
size_t i=0; i<num_my_row; ++i) {
2312 BaseScalar xx = row * h;
2314 val.fastAccessCoeff(
j) =
2315 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
2317 TEST_EQUALITY( x_view(i,0).size(),
VectorSize );
2322 if (ST::magnitude(v.coeff(
j)) < tol)
2323 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2324 if (ST::magnitude(
val.coeff(
j)) < tol)
2325 val.fastAccessCoeff(
j) = BaseScalar(0.0);
2329 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
2342#if defined(HAVE_STOKHOS_AMESOS2)
2352 using Teuchos::ArrayView;
2353 using Teuchos::Array;
2354 using Teuchos::ArrayRCP;
2355 using Teuchos::ParameterList;
2357 typedef typename Storage::value_type BaseScalar;
2360 typedef Teuchos::Comm<int> Tpetra_Comm;
2361 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2362 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2363 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
2364 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2365 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2368 if ( !Kokkos::is_initialized() )
2369 Kokkos::initialize();
2373 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
2374 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2375 RCP<const Tpetra_Map> map =
2376 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2378 RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map,
size_t(3));
2379 Array<GlobalOrdinal> columnIndices(3);
2380 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2381 const size_t num_my_row = myGIDs.size();
2382 for (
size_t i=0; i<num_my_row; ++i) {
2384 if (row == 0 || row == nrow-1) {
2385 columnIndices[0] = row;
2386 graph->insertGlobalIndices(row, columnIndices(0,1));
2389 columnIndices[0] = row-1;
2390 columnIndices[1] = row;
2391 columnIndices[2] = row+1;
2392 graph->insertGlobalIndices(row, columnIndices(0,3));
2395 graph->fillComplete();
2396 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
2399 Array<Scalar> vals(3);
2402 a_val.fastAccessCoeff(
j) =
2403 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
2405 for (
size_t i=0; i<num_my_row; ++i) {
2407 if (row == 0 || row == nrow-1) {
2408 columnIndices[0] = row;
2410 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2413 columnIndices[0] = row-1;
2414 columnIndices[1] = row;
2415 columnIndices[2] = row+1;
2416 vals[0] =
Scalar(-1.0) * a_val;
2417 vals[1] =
Scalar(2.0) * a_val;
2418 vals[2] =
Scalar(-1.0) * a_val;
2419 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2422 matrix->fillComplete();
2425 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2428 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
2431 b_val.fastAccessCoeff(
j) =
2432 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
2434 for (
size_t i=0; i<num_my_row; ++i) {
2436 if (row == 0 || row == nrow-1)
2437 b_view(i,0) =
Scalar(0.0);
2439 b_view(i,0) = -
Scalar(b_val * h * h);
2444 typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
2445 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2446 std::string solver_name;
2447#if defined(HAVE_AMESOS2_BASKER)
2448 solver_name =
"basker";
2449#elif defined(HAVE_AMESOS2_KLU2)
2450 solver_name =
"klu2";
2451#elif defined(HAVE_AMESOS2_SUPERLUDIST)
2452 solver_name =
"superlu_dist";
2453#elif defined(HAVE_AMESOS2_SUPERLUMT)
2454 solver_name =
"superlu_mt";
2455#elif defined(HAVE_AMESOS2_SUPERLU)
2456 solver_name =
"superlu";
2457#elif defined(HAVE_AMESOS2_PARDISO_MKL)
2458 solver_name =
"pardisomkl";
2459#elif defined(HAVE_AMESOS2_LAPACK)
2460 solver_name =
"lapack";
2461#elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
2462 solver_name =
"lapack";
2468 out <<
"Solving linear system with " << solver_name << std::endl;
2469 RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2470 solver_name, matrix, x, b);
2477 solver = Teuchos::null;
2478 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2479 typename ST::magnitudeType tol = 1e-9;
2480 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
2482 for (
size_t i=0; i<num_my_row; ++i) {
2484 BaseScalar xx = row * h;
2486 val.fastAccessCoeff(
j) =
2487 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
2489 TEST_EQUALITY( x_view(i,0).size(),
VectorSize );
2494 if (ST::magnitude(v.coeff(
j)) < tol)
2495 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2496 if (ST::magnitude(
val.coeff(
j)) < tol)
2497 val.fastAccessCoeff(
j) = BaseScalar(0.0);
2501 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
2513#define CRSMATRIX_MP_VECTOR_TESTS_SLGN(S, LO, GO, N) \
2514 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorAdd, S, LO, GO, N ) \
2515 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorDot, S, LO, GO, N ) \
2516 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorAdd, S, LO, GO, N ) \
2517 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDot, S, LO, GO, N ) \
2518 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDotSub, S, LO, GO, N ) \
2519 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixVectorMultiply, S, LO, GO, N ) \
2520 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, S, LO, GO, N ) \
2521 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, WrappedDualView, S, LO, GO, N ) \
2522 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Flatten, S, LO, GO, N ) \
2523 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimpleCG, S, LO, GO, N ) \
2524 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimplePCG_Muelu, S, LO, GO, N ) \
2525 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES, S, LO, GO, N ) \
2526 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, S, LO, GO, N ) \
2527 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, S, LO, GO, N ) \
2528 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, S, LO, GO, N ) \
2529 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, S, LO, GO, N ) \
2530 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosCG_Muelu, S, LO, GO, N ) \
2531 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2, S, LO, GO, N )
2533#define CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N) \
2534 typedef Stokhos::DeviceForNode<N>::type Device; \
2535 typedef Stokhos::StaticFixedStorage<int,double,VectorSize,Device::execution_space> SFS; \
2536 using default_global_ordinal_type = ::Tpetra::Map<>::global_ordinal_type; \
2537 using default_local_ordinal_type = ::Tpetra::Map<>::local_ordinal_type; \
2538 CRSMATRIX_MP_VECTOR_TESTS_SLGN(SFS, default_local_ordinal_type, default_global_ordinal_type, N)
2540#define CRSMATRIX_MP_VECTOR_TESTS_N(N) \
2541 CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N)
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j) - expr2.val(j)
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node)
scalar generate_multi_vector_coefficient(const ordinal nFEM, const ordinal nVec, const ordinal nStoch, const ordinal iColFEM, const ordinal iVec, const ordinal iStoch)
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
Convergence test using the implicit residual norm(s), with an explicit residual norm(s) check for los...
KokkosClassic::DefaultNode::DefaultNodeType Node
bool CG_Solve(const Matrix &A, Vector &x, const Vector &b, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
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< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_map)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &cijk_graph, const CijkType &cijk_dev)
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)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x