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" 87 template <
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;
100 template <
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;
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->getNodeElementList();
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);
159 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
160 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
162 for (
size_t i=0; i<num_my_row; ++i) {
165 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
166 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
171 x1_view = Teuchos::null;
172 x2_view = Teuchos::null;
177 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
178 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
184 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
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].size(),
VectorSize );
208 using Teuchos::ArrayView;
209 using Teuchos::Array;
210 using Teuchos::ArrayRCP;
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->getNodeElementList();
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);
238 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
239 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
241 for (
size_t i=0; i<num_my_row; ++i) {
244 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
245 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
250 x1_view = Teuchos::null;
251 x2_view = Teuchos::null;
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;
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->getNodeElementList();
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);
340 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
341 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
343 for (
size_t i=0; i<num_my_row; ++i) {
345 for (
size_t j=0;
j<ncol; ++
j) {
348 generate_multi_vector_coefficient<BaseScalar,size_t>(
350 val1.fastAccessCoeff(k) = v;
351 val2.fastAccessCoeff(k) = 0.12345 * v;
353 x1_view[
j][i] = val1;
354 x2_view[
j][i] = val2;
357 x1_view = Teuchos::null;
358 x2_view = Teuchos::null;
363 RCP<Tpetra_MultiVector>
y = Tpetra::createMultiVector<Scalar>(map, ncol);
364 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
370 ArrayRCP< ArrayRCP<Scalar> > y_view =
y->get2dViewNonConst();
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;
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->getNodeElementList();
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);
428 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
429 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
431 for (
size_t i=0; i<num_my_row; ++i) {
433 for (
size_t j=0;
j<ncol; ++
j) {
436 generate_multi_vector_coefficient<BaseScalar,size_t>(
438 val1.fastAccessCoeff(k) = v;
439 val2.fastAccessCoeff(k) = 0.12345 * v;
441 x1_view[
j][i] = val1;
442 x2_view[
j][i] = val2;
445 x1_view = Teuchos::null;
446 x2_view = Teuchos::null;
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;
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->getNodeElementList();
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);
543 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
544 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
546 for (
size_t i=0; i<num_my_row; ++i) {
548 for (
size_t j=0;
j<ncol; ++
j) {
551 generate_multi_vector_coefficient<BaseScalar,size_t>(
553 val1.fastAccessCoeff(k) = v;
554 val2.fastAccessCoeff(k) = 0.12345 * v;
556 x1_view[
j][i] = val1;
557 x2_view[
j][i] = val2;
560 x1_view = Teuchos::null;
561 x2_view = Teuchos::null;
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;
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), Tpetra::StaticProfile));
661 Array<GlobalOrdinal> columnIndices(2);
662 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
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);
703 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
704 for (
size_t i=0; i<num_my_row; ++i) {
707 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
711 x_view = Teuchos::null;
720 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
721 matrix->apply(*
x, *
y);
727 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
728 BaseScalar tol = 1.0e-14;
729 for (
size_t i=0; i<num_my_row; ++i) {
732 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
734 val.fastAccessCoeff(
j) = v*v;
738 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
740 val.fastAccessCoeff(
j) += v*v;
743 TEST_EQUALITY( y_view[i].size(),
VectorSize );
757 using Teuchos::ArrayView;
758 using Teuchos::Array;
759 using Teuchos::ArrayRCP;
764 typedef Teuchos::Comm<int> Tpetra_Comm;
765 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
766 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
767 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
768 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
771 if ( !Kokkos::is_initialized() )
772 Kokkos::initialize();
776 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
777 RCP<const Tpetra_Map> map =
778 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
780 RCP<Tpetra_CrsGraph> graph =
781 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
782 Array<GlobalOrdinal> columnIndices(2);
783 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
784 const size_t num_my_row = myGIDs.size();
785 for (
size_t i=0; i<num_my_row; ++i) {
787 columnIndices[0] = row;
790 columnIndices[1] = row+1;
793 graph->insertGlobalIndices(row, columnIndices(0,ncol));
795 graph->fillComplete();
796 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
799 Array<Scalar> vals(2);
801 for (
size_t i=0; i<num_my_row; ++i) {
803 columnIndices[0] = row;
805 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
811 columnIndices[1] = row+1;
813 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
818 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
820 matrix->fillComplete();
824 RCP<Tpetra_MultiVector>
x = Tpetra::createMultiVector<Scalar>(map, ncol);
825 ArrayRCP< ArrayRCP<Scalar> > x_view =
x->get2dViewNonConst();
826 for (
size_t i=0; i<num_my_row; ++i) {
828 for (
size_t j=0;
j<ncol; ++
j) {
831 generate_multi_vector_coefficient<BaseScalar,size_t>(
833 val.fastAccessCoeff(k) = v;
838 x_view = Teuchos::null;
847 RCP<Tpetra_MultiVector>
y = Tpetra::createMultiVector<Scalar>(map, ncol);
848 matrix->apply(*
x, *
y);
854 ArrayRCP< ArrayRCP<Scalar> > y_view =
y->get2dViewNonConst();
855 BaseScalar tol = 1.0e-14;
856 for (
size_t i=0; i<num_my_row; ++i) {
858 for (
size_t j=0;
j<ncol; ++
j) {
860 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
862 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
864 val.fastAccessCoeff(k) = v1*v2;
868 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
870 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
872 val.fastAccessCoeff(k) += v1*v2;
878 val.fastAccessCoeff(k), tol );
891 using Teuchos::ArrayView;
892 using Teuchos::Array;
893 using Teuchos::ArrayRCP;
898 typedef Teuchos::Comm<int> Tpetra_Comm;
899 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
900 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
901 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
902 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
904 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
905 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
908 if ( !Kokkos::is_initialized() )
909 Kokkos::initialize();
913 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
914 RCP<const Tpetra_Map> map =
915 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
917 RCP<Tpetra_CrsGraph> graph =
918 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
919 Array<GlobalOrdinal> columnIndices(2);
920 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
921 const size_t num_my_row = myGIDs.size();
922 for (
size_t i=0; i<num_my_row; ++i) {
924 columnIndices[0] = row;
927 columnIndices[1] = row+1;
930 graph->insertGlobalIndices(row, columnIndices(0,ncol));
932 graph->fillComplete();
933 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
936 Array<Scalar> vals(2);
938 for (
size_t i=0; i<num_my_row; ++i) {
940 columnIndices[0] = row;
942 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
948 columnIndices[1] = row+1;
950 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
955 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
957 matrix->fillComplete();
960 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
961 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
962 for (
size_t i=0; i<num_my_row; ++i) {
965 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
971 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
972 matrix->apply(*
x, *
y);
975 RCP<const Tpetra_Map> flat_x_map, flat_y_map;
976 RCP<const Tpetra_CrsGraph> flat_graph =
978 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
982 RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
983 RCP<Flat_Tpetra_Vector> flat_x =
985 RCP<Flat_Tpetra_Vector> flat_y =
987 flat_matrix->apply(*flat_x, *flat_y);
993 BaseScalar tol = 1.0e-14;
994 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
995 ArrayRCP<Scalar> y2_view = y2->get1dViewNonConst();
996 for (
size_t i=0; i<num_my_row; ++i) {
997 TEST_EQUALITY( y_view[i].size(),
VectorSize );
998 TEST_EQUALITY( y2_view[i].size(),
VectorSize );
1018 using Teuchos::ArrayView;
1019 using Teuchos::Array;
1020 using Teuchos::ArrayRCP;
1025 using DualViewType = Kokkos::DualView<Scalar*, typename Node::device_type>;
1026 using WDV = Tpetra::Details::WrappedDualView<DualViewType>;
1027 using values_view =
typename DualViewType::t_dev;
1030 if ( !Kokkos::is_initialized() )
1031 Kokkos::initialize();
1035 values_view myView(
"emptyTestView", 0);
1038 size_t use_h = wdv.getHostView(Tpetra::Access::ReadOnly).use_count();
1039 size_t use_d = wdv.getDeviceView(Tpetra::Access::ReadOnly).use_count();
1042 TEST_EQUALITY(use_h, use_d);
1053 using Teuchos::ArrayView;
1054 using Teuchos::Array;
1055 using Teuchos::ArrayRCP;
1056 using Teuchos::ParameterList;
1061 typedef Teuchos::Comm<int> Tpetra_Comm;
1062 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1063 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1064 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1065 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1068 if ( !Kokkos::is_initialized() )
1069 Kokkos::initialize();
1073 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1074 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1075 RCP<const Tpetra_Map> map =
1076 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1078 RCP<Tpetra_CrsGraph> graph =
1079 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1080 Array<GlobalOrdinal> columnIndices(3);
1081 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1082 const size_t num_my_row = myGIDs.size();
1083 for (
size_t i=0; i<num_my_row; ++i) {
1085 if (row == 0 || row == nrow-1) {
1086 columnIndices[0] = row;
1087 graph->insertGlobalIndices(row, columnIndices(0,1));
1090 columnIndices[0] = row-1;
1091 columnIndices[1] = row;
1092 columnIndices[2] = row+1;
1093 graph->insertGlobalIndices(row, columnIndices(0,3));
1096 graph->fillComplete();
1097 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1100 Array<Scalar> vals(3);
1103 a_val.fastAccessCoeff(
j) =
1104 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1106 for (
size_t i=0; i<num_my_row; ++i) {
1108 if (row == 0 || row == nrow-1) {
1109 columnIndices[0] = row;
1111 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1114 columnIndices[0] = row-1;
1115 columnIndices[1] = row;
1116 columnIndices[2] = row+1;
1117 vals[0] =
Scalar(-1.0) * a_val;
1118 vals[1] =
Scalar(2.0) * a_val;
1119 vals[2] =
Scalar(-1.0) * a_val;
1120 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1123 matrix->fillComplete();
1125 matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
1126 Teuchos::VERB_EXTREME);
1129 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1130 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1133 b_val.fastAccessCoeff(
j) =
1134 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1136 for (
size_t i=0; i<num_my_row; ++i) {
1138 if (row == 0 || row == nrow-1)
1141 b_view[i] = -
Scalar(b_val * h * h);
1145 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1146 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1147 typedef typename BST::mag_type base_mag_type;
1148 typedef typename Tpetra_Vector::mag_type mag_type;
1149 base_mag_type btol = 1e-9;
1150 mag_type tol = btol;
1153 out.getOStream().get());
1154 TEST_EQUALITY_CONST( solved,
true );
1161 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1163 for (
size_t i=0; i<num_my_row; ++i) {
1165 BaseScalar xx = row * h;
1167 val.fastAccessCoeff(
j) =
1168 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1170 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1176 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1178 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1182 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), btol);
1187 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2) 1197 using Teuchos::ArrayView;
1198 using Teuchos::Array;
1199 using Teuchos::ArrayRCP;
1200 using Teuchos::ParameterList;
1201 using Teuchos::getParametersFromXmlFile;
1206 typedef Teuchos::Comm<int> Tpetra_Comm;
1207 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1208 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1209 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1210 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1213 if ( !Kokkos::is_initialized() )
1214 Kokkos::initialize();
1218 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1219 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1220 RCP<const Tpetra_Map> map =
1221 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1223 RCP<Tpetra_CrsGraph> graph =
1224 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1225 Array<GlobalOrdinal> columnIndices(3);
1226 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1227 const size_t num_my_row = myGIDs.size();
1228 for (
size_t i=0; i<num_my_row; ++i) {
1230 if (row == 0 || row == nrow-1) {
1231 columnIndices[0] = row;
1232 graph->insertGlobalIndices(row, columnIndices(0,1));
1235 columnIndices[0] = row-1;
1236 columnIndices[1] = row;
1237 columnIndices[2] = row+1;
1238 graph->insertGlobalIndices(row, columnIndices(0,3));
1241 graph->fillComplete();
1242 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1245 Array<Scalar> vals(3);
1248 a_val.fastAccessCoeff(
j) =
1249 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1251 for (
size_t i=0; i<num_my_row; ++i) {
1253 if (row == 0 || row == nrow-1) {
1254 columnIndices[0] = row;
1256 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1259 columnIndices[0] = row-1;
1260 columnIndices[1] = row;
1261 columnIndices[2] = row+1;
1262 vals[0] =
Scalar(-1.0) * a_val;
1263 vals[1] =
Scalar(2.0) * a_val;
1264 vals[2] =
Scalar(-1.0) * a_val;
1265 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1268 matrix->fillComplete();
1271 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1272 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1275 b_val.fastAccessCoeff(
j) =
1276 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1278 for (
size_t i=0; i<num_my_row; ++i) {
1280 if (row == 0 || row == nrow-1)
1283 b_view[i] = -
Scalar(b_val * h * h);
1287 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1288 RCP<OP> matrix_op = matrix;
1289 RCP<ParameterList> muelu_params =
1290 getParametersFromXmlFile(
"muelu_cheby.xml");
1292 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1295 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1296 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1297 typedef typename BST::mag_type base_mag_type;
1298 typedef typename Tpetra_Vector::mag_type mag_type;
1299 base_mag_type btol = 1e-9;
1300 mag_type tol = btol;
1303 out.getOStream().get());
1304 TEST_EQUALITY_CONST( solved,
true );
1311 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1313 for (
size_t i=0; i<num_my_row; ++i) {
1315 BaseScalar xx = row * h;
1317 val.fastAccessCoeff(
j) =
1318 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1320 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1325 if (BST::magnitude(v.coeff(
j)) < btol)
1326 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1327 if (BST::magnitude(
val.coeff(
j)) < btol)
1328 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1332 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), btol);
1344 #if defined(HAVE_STOKHOS_BELOS) 1354 using Teuchos::ArrayView;
1355 using Teuchos::Array;
1356 using Teuchos::ArrayRCP;
1357 using Teuchos::ParameterList;
1362 typedef Teuchos::Comm<int> Tpetra_Comm;
1363 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1364 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1365 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1366 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1369 if ( !Kokkos::is_initialized() )
1370 Kokkos::initialize();
1374 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1375 RCP<const Tpetra_Map> map =
1376 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1378 RCP<Tpetra_CrsGraph> graph =
1379 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1380 Array<GlobalOrdinal> columnIndices(2);
1381 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1382 const size_t num_my_row = myGIDs.size();
1383 for (
size_t i=0; i<num_my_row; ++i) {
1385 columnIndices[0] = row;
1387 if (row != nrow-1) {
1388 columnIndices[1] = row+1;
1391 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1393 graph->fillComplete();
1394 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1397 Array<Scalar> vals(2);
1399 for (
size_t i=0; i<num_my_row; ++i) {
1401 columnIndices[0] = row;
1403 val.fastAccessCoeff(
j) =
j+1;
1407 if (row != nrow-1) {
1408 columnIndices[1] = row+1;
1410 val.fastAccessCoeff(
j) =
j+1;
1414 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1416 matrix->fillComplete();
1419 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1420 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1421 for (
size_t i=0; i<num_my_row; ++i) {
1426 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1427 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1428 typedef BaseScalar BelosScalar;
1430 typedef Scalar BelosScalar;
1432 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1433 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1434 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1435 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1436 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1437 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1438 typename ST::magnitudeType tol = 1e-9;
1439 belosParams->set(
"Flexible Gmres",
false);
1440 belosParams->set(
"Num Blocks", 100);
1441 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1442 belosParams->set(
"Maximum Iterations", 100);
1443 belosParams->set(
"Verbosity", 33);
1444 belosParams->set(
"Output Style", 1);
1445 belosParams->set(
"Output Frequency", 1);
1446 belosParams->set(
"Output Stream", out.getOStream());
1447 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1448 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1449 problem->setProblem();
1450 Belos::ReturnType ret = solver->solve();
1451 TEST_EQUALITY_CONST( ret, Belos::Converged );
1463 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1464 for (
size_t i=0; i<num_my_row; ++i) {
1468 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
1473 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1478 if (ST::magnitude(v.coeff(
j)) < tol)
1479 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1483 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1495 using Teuchos::tuple;
1496 using Teuchos::ArrayView;
1497 using Teuchos::Array;
1498 using Teuchos::ArrayRCP;
1499 using Teuchos::ParameterList;
1504 typedef Teuchos::Comm<int> Tpetra_Comm;
1505 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1506 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1507 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1508 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1511 if ( !Kokkos::is_initialized() )
1512 Kokkos::initialize();
1516 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1517 RCP<const Tpetra_Map> map =
1518 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1520 RCP<Tpetra_CrsGraph> graph =
1521 rcp(
new Tpetra_CrsGraph(map,
size_t(1), Tpetra::StaticProfile));
1522 Array<GlobalOrdinal> columnIndices(1);
1523 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1524 const size_t num_my_row = myGIDs.size();
1525 for (
size_t i=0; i<num_my_row; ++i) {
1527 columnIndices[0] = row;
1529 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1531 graph->fillComplete();
1532 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1534 Array<Scalar> vals(1);
1535 for (
size_t i=0; i<num_my_row; ++i) {
1537 columnIndices[0] = row;
1539 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1541 matrix->fillComplete();
1550 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1551 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1552 for (
size_t i=0; i<num_my_row; ++i) {
1557 b_view[i].fastAccessCoeff(
j) = BaseScalar(row+1);
1561 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1562 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1563 typedef BaseScalar BelosScalar;
1565 typedef Scalar BelosScalar;
1567 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1568 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1569 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1570 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1571 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1572 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1573 typename ST::magnitudeType tol = 1e-9;
1574 belosParams->set(
"Flexible Gmres",
false);
1575 belosParams->set(
"Num Blocks", 100);
1576 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1577 belosParams->set(
"Maximum Iterations", 100);
1578 belosParams->set(
"Verbosity", 33);
1579 belosParams->set(
"Output Style", 1);
1580 belosParams->set(
"Output Frequency", 1);
1581 belosParams->set(
"Output Stream", out.getOStream());
1582 belosParams->set(
"Orthogonalization",
"DGKS");
1583 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1584 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1585 problem->setProblem();
1586 Belos::ReturnType ret = solver->solve();
1587 TEST_EQUALITY_CONST( ret, Belos::Converged );
1589 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT 1590 int numItersExpected = nrow;
1591 int numIters = solver->getNumIters();
1592 out <<
"numIters = " << numIters << std::endl;
1593 TEST_EQUALITY( numIters, numItersExpected);
1596 std::vector<int> ensemble_iterations =
1598 out <<
"Ensemble iterations = ";
1599 for (
auto ensemble_iteration : ensemble_iterations)
1600 out << ensemble_iteration <<
" ";
1605 TEST_EQUALITY(
int(
j+1+nrow-
VectorSize), ensemble_iterations[
j]);
1608 TEST_EQUALITY(
int(0), ensemble_iterations[
j]);
1620 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1621 for (
size_t i=0; i<num_my_row; ++i) {
1626 if (ST::magnitude(v.coeff(
j)) < tol)
1627 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1634 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1637 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1649 using Teuchos::tuple;
1650 using Teuchos::ArrayView;
1651 using Teuchos::Array;
1652 using Teuchos::ArrayRCP;
1653 using Teuchos::ParameterList;
1658 typedef Teuchos::Comm<int> Tpetra_Comm;
1659 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1660 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1661 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1662 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1665 if ( !Kokkos::is_initialized() )
1666 Kokkos::initialize();
1670 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1671 RCP<const Tpetra_Map> map =
1672 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1674 RCP<Tpetra_CrsGraph> graph =
1675 rcp(
new Tpetra_CrsGraph(map,
size_t(1), Tpetra::StaticProfile));
1676 Array<GlobalOrdinal> columnIndices(1);
1677 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1678 const size_t num_my_row = myGIDs.size();
1679 for (
size_t i=0; i<num_my_row; ++i) {
1681 columnIndices[0] = row;
1683 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1685 graph->fillComplete();
1686 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1688 Array<Scalar> vals(1);
1689 for (
size_t i=0; i<num_my_row; ++i) {
1691 columnIndices[0] = row;
1693 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1695 matrix->fillComplete();
1704 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1705 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1706 for (
size_t i=0; i<num_my_row; ++i) {
1711 b_view[i].fastAccessCoeff(
j) = BaseScalar(row+1);
1715 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1716 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1717 typedef BaseScalar BelosScalar;
1719 typedef Scalar BelosScalar;
1721 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1722 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1723 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1724 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1725 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1726 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1727 typename ST::magnitudeType tol = 1e-9;
1728 belosParams->set(
"Flexible Gmres",
false);
1729 belosParams->set(
"Num Blocks", 100);
1730 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1731 belosParams->set(
"Maximum Iterations", 100);
1732 belosParams->set(
"Verbosity", 33);
1733 belosParams->set(
"Output Style", 1);
1734 belosParams->set(
"Output Frequency", 1);
1735 belosParams->set(
"Output Stream", out.getOStream());
1736 belosParams->set(
"Orthogonalization",
"ICGS");
1737 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1738 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1739 problem->setProblem();
1740 Belos::ReturnType ret = solver->solve();
1741 TEST_EQUALITY_CONST( ret, Belos::Converged );
1743 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT 1744 int numItersExpected = nrow;
1745 int numIters = solver->getNumIters();
1746 out <<
"numIters = " << numIters << std::endl;
1747 TEST_EQUALITY( numIters, numItersExpected);
1750 std::vector<int> ensemble_iterations =
1752 out <<
"Ensemble iterations = ";
1753 for (
auto ensemble_iteration : ensemble_iterations)
1754 out << ensemble_iteration <<
" ";
1759 TEST_EQUALITY(
int(
j+1+nrow-
VectorSize), ensemble_iterations[
j]);
1762 TEST_EQUALITY(
int(0), ensemble_iterations[
j]);
1774 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1775 for (
size_t i=0; i<num_my_row; ++i) {
1780 if (ST::magnitude(v.coeff(
j)) < tol)
1781 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1788 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1791 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1803 using Teuchos::tuple;
1804 using Teuchos::ArrayView;
1805 using Teuchos::Array;
1806 using Teuchos::ArrayRCP;
1807 using Teuchos::ParameterList;
1812 typedef Teuchos::Comm<int> Tpetra_Comm;
1813 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1814 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1815 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1816 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1819 if ( !Kokkos::is_initialized() )
1820 Kokkos::initialize();
1824 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1825 RCP<const Tpetra_Map> map =
1826 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1828 RCP<Tpetra_CrsGraph> graph =
1829 rcp(
new Tpetra_CrsGraph(map,
size_t(1), Tpetra::StaticProfile));
1830 Array<GlobalOrdinal> columnIndices(1);
1831 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1832 const size_t num_my_row = myGIDs.size();
1833 for (
size_t i=0; i<num_my_row; ++i) {
1835 columnIndices[0] = row;
1837 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1839 graph->fillComplete();
1840 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1842 Array<Scalar> vals(1);
1843 for (
size_t i=0; i<num_my_row; ++i) {
1845 columnIndices[0] = row;
1847 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1849 matrix->fillComplete();
1858 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1859 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1860 for (
size_t i=0; i<num_my_row; ++i) {
1865 b_view[i].fastAccessCoeff(
j) = BaseScalar(row+1);
1869 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1870 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1871 typedef BaseScalar BelosScalar;
1873 typedef Scalar BelosScalar;
1875 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1876 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1877 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1878 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1879 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1880 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1881 typename ST::magnitudeType tol = 1e-9;
1882 belosParams->set(
"Flexible Gmres",
false);
1883 belosParams->set(
"Num Blocks", 100);
1884 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1885 belosParams->set(
"Maximum Iterations", 100);
1886 belosParams->set(
"Verbosity", 33);
1887 belosParams->set(
"Output Style", 1);
1888 belosParams->set(
"Output Frequency", 1);
1889 belosParams->set(
"Output Stream", out.getOStream());
1890 belosParams->set(
"Orthogonalization",
"IMGS");
1891 RCP<Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP> > solver =
1892 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1893 problem->setProblem();
1894 Belos::ReturnType ret = solver->solve();
1895 TEST_EQUALITY_CONST( ret, Belos::Converged );
1897 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT 1898 int numItersExpected = nrow;
1899 int numIters = solver->getNumIters();
1900 out <<
"numIters = " << numIters << std::endl;
1901 TEST_EQUALITY( numIters, numItersExpected);
1904 std::vector<int> ensemble_iterations =
1906 out <<
"Ensemble iterations = ";
1907 for (
auto ensemble_iteration : ensemble_iterations)
1908 out << ensemble_iteration <<
" ";
1913 TEST_EQUALITY(
int(
j+1+nrow-
VectorSize), ensemble_iterations[
j]);
1916 TEST_EQUALITY(
int(0), ensemble_iterations[
j]);
1928 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1929 for (
size_t i=0; i<num_my_row; ++i) {
1934 if (ST::magnitude(v.coeff(
j)) < tol)
1935 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1942 val.fastAccessCoeff(
j) = BaseScalar(1.0);
1945 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1969 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) 1980 using Teuchos::ArrayView;
1981 using Teuchos::Array;
1982 using Teuchos::ArrayRCP;
1983 using Teuchos::ParameterList;
1988 typedef Teuchos::Comm<int> Tpetra_Comm;
1989 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1990 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1991 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1992 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1995 if ( !Kokkos::is_initialized() )
1996 Kokkos::initialize();
2000 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2001 RCP<const Tpetra_Map> map =
2002 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2004 RCP<Tpetra_CrsGraph> graph =
2005 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
2006 Array<GlobalOrdinal> columnIndices(2);
2007 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
2008 const size_t num_my_row = myGIDs.size();
2009 for (
size_t i=0; i<num_my_row; ++i) {
2011 columnIndices[0] = row;
2013 if (row != nrow-1) {
2014 columnIndices[1] = row+1;
2017 graph->insertGlobalIndices(row, columnIndices(0,ncol));
2019 graph->fillComplete();
2020 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
2023 Array<Scalar> vals(2);
2025 for (
size_t i=0; i<num_my_row; ++i) {
2027 columnIndices[0] = row;
2029 val.fastAccessCoeff(
j) =
j+1;
2033 if (row != nrow-1) {
2034 columnIndices[1] = row+1;
2036 val.fastAccessCoeff(
j) =
j+1;
2040 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
2042 matrix->fillComplete();
2045 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2046 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2047 for (
size_t i=0; i<num_my_row; ++i) {
2052 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
2053 Ifpack2::Factory factory;
2054 RCP<Prec> M = factory.create<Tpetra_CrsMatrix>(
"RILUK", matrix);
2059 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2060 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 2061 typedef BaseScalar BelosScalar;
2063 typedef Scalar BelosScalar;
2065 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2066 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2067 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2068 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2069 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
2070 problem->setRightPrec(M);
2071 problem->setProblem();
2072 RCP<ParameterList> belosParams = rcp(
new ParameterList);
2073 typename ST::magnitudeType tol = 1e-9;
2074 belosParams->set(
"Flexible Gmres",
false);
2075 belosParams->set(
"Num Blocks", 100);
2076 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
2077 belosParams->set(
"Maximum Iterations", 100);
2078 belosParams->set(
"Verbosity", 33);
2079 belosParams->set(
"Output Style", 1);
2080 belosParams->set(
"Output Frequency", 1);
2081 belosParams->set(
"Output Stream", out.getOStream());
2083 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2084 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2085 Belos::ReturnType ret = solver->solve();
2086 TEST_EQUALITY_CONST( ret, Belos::Converged );
2098 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
2099 for (
size_t i=0; i<num_my_row; ++i) {
2103 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
2108 TEST_EQUALITY( x_view[i].size(),
VectorSize );
2113 if (ST::magnitude(v.coeff(
j)) < tol)
2114 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2118 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
2130 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) 2140 using Teuchos::ArrayView;
2141 using Teuchos::Array;
2142 using Teuchos::ArrayRCP;
2143 using Teuchos::ParameterList;
2144 using Teuchos::getParametersFromXmlFile;
2149 typedef Teuchos::Comm<int> Tpetra_Comm;
2150 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2151 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2152 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2153 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2156 if ( !Kokkos::is_initialized() )
2157 Kokkos::initialize();
2161 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
2162 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2163 RCP<const Tpetra_Map> map =
2164 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2166 RCP<Tpetra_CrsGraph> graph =
2167 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
2168 Array<GlobalOrdinal> columnIndices(3);
2169 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
2170 const size_t num_my_row = myGIDs.size();
2171 for (
size_t i=0; i<num_my_row; ++i) {
2173 if (row == 0 || row == nrow-1) {
2174 columnIndices[0] = row;
2175 graph->insertGlobalIndices(row, columnIndices(0,1));
2178 columnIndices[0] = row-1;
2179 columnIndices[1] = row;
2180 columnIndices[2] = row+1;
2181 graph->insertGlobalIndices(row, columnIndices(0,3));
2184 graph->fillComplete();
2185 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
2188 Array<Scalar> vals(3);
2191 a_val.fastAccessCoeff(
j) =
2192 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
2194 for (
size_t i=0; i<num_my_row; ++i) {
2196 if (row == 0 || row == nrow-1) {
2197 columnIndices[0] = row;
2199 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2202 columnIndices[0] = row-1;
2203 columnIndices[1] = row;
2204 columnIndices[2] = row+1;
2205 vals[0] =
Scalar(-1.0) * a_val;
2206 vals[1] =
Scalar(2.0) * a_val;
2207 vals[2] =
Scalar(-1.0) * a_val;
2208 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2211 matrix->fillComplete();
2214 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2215 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
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)
2226 b_view[i] = -
Scalar(b_val * h * h);
2230 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2231 RCP<ParameterList> muelu_params =
2232 getParametersFromXmlFile(
"muelu_cheby.xml");
2233 RCP<OP> matrix_op = matrix;
2235 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
2238 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2239 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 2240 typedef BaseScalar BelosScalar;
2242 typedef Scalar BelosScalar;
2244 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2245 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2246 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2247 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
2248 problem->setRightPrec(M);
2249 problem->setProblem();
2250 RCP<ParameterList> belosParams = rcp(
new ParameterList);
2251 typename ST::magnitudeType tol = 1e-9;
2252 belosParams->set(
"Num Blocks", 100);
2253 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
2254 belosParams->set(
"Maximum Iterations", 100);
2255 belosParams->set(
"Verbosity", 33);
2256 belosParams->set(
"Output Style", 1);
2257 belosParams->set(
"Output Frequency", 1);
2258 belosParams->set(
"Output Stream", out.getOStream());
2261 belosParams->set(
"Implicit Residual Scaling",
"None");
2263 RCP<Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true> > solver =
2264 rcp(
new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true>(problem, belosParams));
2265 Belos::ReturnType ret = solver->solve();
2266 TEST_EQUALITY_CONST( ret, Belos::Converged );
2268 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT 2270 std::vector<int> ensemble_iterations =
2271 solver->getResidualStatusTest()->getEnsembleIterations();
2272 out <<
"Ensemble iterations = ";
2274 out << ensemble_iterations[i] <<
" ";
2283 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
2285 for (
size_t i=0; i<num_my_row; ++i) {
2287 BaseScalar xx = row * h;
2289 val.fastAccessCoeff(
j) =
2290 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
2292 TEST_EQUALITY( x_view[i].size(),
VectorSize );
2297 if (ST::magnitude(v.coeff(
j)) < tol)
2298 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2299 if (ST::magnitude(
val.coeff(
j)) < tol)
2300 val.fastAccessCoeff(
j) = BaseScalar(0.0);
2304 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
2317 #if defined(HAVE_STOKHOS_AMESOS2) 2327 using Teuchos::ArrayView;
2328 using Teuchos::Array;
2329 using Teuchos::ArrayRCP;
2330 using Teuchos::ParameterList;
2335 typedef Teuchos::Comm<int> Tpetra_Comm;
2336 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2337 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2338 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
2339 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2340 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2343 if ( !Kokkos::is_initialized() )
2344 Kokkos::initialize();
2348 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
2349 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2350 RCP<const Tpetra_Map> map =
2351 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2353 RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map,
size_t(3));
2354 Array<GlobalOrdinal> columnIndices(3);
2355 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
2356 const size_t num_my_row = myGIDs.size();
2357 for (
size_t i=0; i<num_my_row; ++i) {
2359 if (row == 0 || row == nrow-1) {
2360 columnIndices[0] = row;
2361 graph->insertGlobalIndices(row, columnIndices(0,1));
2364 columnIndices[0] = row-1;
2365 columnIndices[1] = row;
2366 columnIndices[2] = row+1;
2367 graph->insertGlobalIndices(row, columnIndices(0,3));
2370 graph->fillComplete();
2371 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
2374 Array<Scalar> vals(3);
2377 a_val.fastAccessCoeff(
j) =
2378 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
2380 for (
size_t i=0; i<num_my_row; ++i) {
2382 if (row == 0 || row == nrow-1) {
2383 columnIndices[0] = row;
2385 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2388 columnIndices[0] = row-1;
2389 columnIndices[1] = row;
2390 columnIndices[2] = row+1;
2391 vals[0] =
Scalar(-1.0) * a_val;
2392 vals[1] =
Scalar(2.0) * a_val;
2393 vals[2] =
Scalar(-1.0) * a_val;
2394 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2397 matrix->fillComplete();
2400 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2401 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2404 b_val.fastAccessCoeff(
j) =
2405 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
2407 for (
size_t i=0; i<num_my_row; ++i) {
2409 if (row == 0 || row == nrow-1)
2412 b_view[i] = -
Scalar(b_val * h * h);
2416 typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
2417 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2418 std::string solver_name;
2419 #if defined(HAVE_AMESOS2_BASKER) 2420 solver_name =
"basker";
2421 #elif defined(HAVE_AMESOS2_KLU2) 2422 solver_name =
"klu2";
2423 #elif defined(HAVE_AMESOS2_SUPERLUDIST) 2424 solver_name =
"superlu_dist";
2425 #elif defined(HAVE_AMESOS2_SUPERLUMT) 2426 solver_name =
"superlu_mt";
2427 #elif defined(HAVE_AMESOS2_SUPERLU) 2428 solver_name =
"superlu";
2429 #elif defined(HAVE_AMESOS2_PARDISO_MKL) 2430 solver_name =
"pardisomkl";
2431 #elif defined(HAVE_AMESOS2_LAPACK) 2432 solver_name =
"lapack";
2433 #elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL) 2434 solver_name =
"lapack";
2440 out <<
"Solving linear system with " << solver_name << std::endl;
2441 RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2442 solver_name, matrix,
x, b);
2449 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2450 typename ST::magnitudeType tol = 1e-9;
2451 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
2453 for (
size_t i=0; i<num_my_row; ++i) {
2455 BaseScalar xx = row * h;
2457 val.fastAccessCoeff(
j) =
2458 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
2460 TEST_EQUALITY( x_view[i].size(),
VectorSize );
2465 if (ST::magnitude(v.coeff(
j)) < tol)
2466 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2467 if (ST::magnitude(
val.coeff(
j)) < tol)
2468 val.fastAccessCoeff(
j) = BaseScalar(0.0);
2472 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
2484 #define CRSMATRIX_MP_VECTOR_TESTS_SLGN(S, LO, GO, N) \ 2485 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorAdd, S, LO, GO, N ) \ 2486 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorDot, S, LO, GO, N ) \ 2487 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorAdd, S, LO, GO, N ) \ 2488 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDot, S, LO, GO, N ) \ 2489 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDotSub, S, LO, GO, N ) \ 2490 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixVectorMultiply, S, LO, GO, N ) \ 2491 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, S, LO, GO, N ) \ 2492 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, WrappedDualView, S, LO, GO, N ) \ 2493 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Flatten, S, LO, GO, N ) \ 2494 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimpleCG, S, LO, GO, N ) \ 2495 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimplePCG_Muelu, S, LO, GO, N ) \ 2496 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES, S, LO, GO, N ) \ 2497 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_DGKS, S, LO, GO, N ) \ 2498 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_ICGS, S, LO, GO, N ) \ 2499 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_IMGS, S, LO, GO, N ) \ 2500 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, S, LO, GO, N ) \ 2501 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosCG_Muelu, S, LO, GO, N ) \ 2502 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2, S, LO, GO, N ) 2504 #define CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N) \ 2505 typedef Stokhos::DeviceForNode<N>::type Device; \ 2506 typedef Stokhos::StaticFixedStorage<int,double,VectorSize,Device::execution_space> SFS; \ 2507 using default_global_ordinal_type = ::Tpetra::Map<>::global_ordinal_type; \ 2508 using default_local_ordinal_type = ::Tpetra::Map<>::local_ordinal_type; \ 2509 CRSMATRIX_MP_VECTOR_TESTS_SLGN(SFS, default_local_ordinal_type, default_global_ordinal_type, N) 2511 #define CRSMATRIX_MP_VECTOR_TESTS_N(N) \ 2512 CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N) 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...
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j) - expr2.val(j)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec_const, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
bool PCG_Solve(const Matrix &A, Vector &x, const Vector &b, const Prec &M, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
KokkosClassic::DefaultNode::DefaultNodeType Node
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
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)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &flat_graph, const LocalOrdinal block_size)
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)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)