42 #include "Teuchos_UnitTestHelpers.hpp" 46 #include "Teuchos_XMLParameterListCoreHelpers.hpp" 51 #include "Tpetra_ConfigDefs.hpp" 52 #include "Tpetra_DefaultPlatform.hpp" 53 #include "Tpetra_Map.hpp" 54 #include "Tpetra_MultiVector.hpp" 55 #include "Tpetra_Vector.hpp" 56 #include "Tpetra_CrsGraph.hpp" 57 #include "Tpetra_CrsMatrix.hpp" 61 #ifdef HAVE_STOKHOS_BELOS 63 #include "BelosLinearProblem.hpp" 64 #include "BelosPseudoBlockGmresSolMgr.hpp" 65 #include "BelosPseudoBlockCGSolMgr.hpp" 69 #ifdef HAVE_STOKHOS_IFPACK2 71 #include "Ifpack2_Factory.hpp" 75 #ifdef HAVE_STOKHOS_MUELU 77 #include "MueLu_CreateTpetraPreconditioner.hpp" 81 #ifdef HAVE_STOKHOS_AMESOS2 83 #include "Amesos2_Factory.hpp" 86 template <
typename scalar,
typename ordinal>
90 const ordinal iColFEM,
91 const ordinal iStoch )
93 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
94 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
95 return X_fem + X_stoch;
99 template <
typename scalar,
typename ordinal>
103 const ordinal nStoch,
104 const ordinal iColFEM,
106 const ordinal iStoch)
108 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
109 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
110 const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
111 return X_fem + X_stoch + X_col;
120 #if defined(__CUDACC__) 134 using Teuchos::ArrayView;
135 using Teuchos::Array;
136 using Teuchos::ArrayRCP;
142 typedef Teuchos::Comm<int> Tpetra_Comm;
143 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
147 if (!Kokkos::HostSpace::execution_space::is_initialized())
148 Kokkos::HostSpace::execution_space::initialize();
149 if (!Device::is_initialized())
150 Device::initialize();
153 RCP<const Tpetra_Comm> comm =
154 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
158 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
159 RCP<const Tpetra_Map> map =
160 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
162 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
163 const size_t num_my_row = myGIDs.size();
166 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
167 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
168 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
169 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
171 for (
size_t i=0; i<num_my_row; ++i) {
174 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
175 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
180 x1_view = Teuchos::null;
181 x2_view = Teuchos::null;
186 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
187 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
193 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
195 BaseScalar tol = 1.0e-14;
196 for (
size_t i=0; i<num_my_row; ++i) {
199 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
201 val.fastAccessCoeff(
j) = alpha.coeff(
j)*v + 0.12345*beta.coeff(
j)*v;
203 TEST_EQUALITY( y_view[i].size(),
VectorSize );
205 TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(
j),
val.fastAccessCoeff(
j), tol );
217 using Teuchos::ArrayView;
218 using Teuchos::Array;
219 using Teuchos::ArrayRCP;
225 typedef Teuchos::Comm<int> Tpetra_Comm;
226 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
228 typedef typename Tpetra_Vector::dot_type dot_type;
231 if (!Kokkos::HostSpace::execution_space::is_initialized())
232 Kokkos::HostSpace::execution_space::initialize();
233 if (!Device::is_initialized())
234 Device::initialize();
237 RCP<const Tpetra_Comm> comm =
238 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
242 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
243 RCP<const Tpetra_Map> map =
244 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
246 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
247 const size_t num_my_row = myGIDs.size();
250 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
251 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
252 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
253 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
255 for (
size_t i=0; i<num_my_row; ++i) {
258 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
259 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
264 x1_view = Teuchos::null;
265 x2_view = Teuchos::null;
268 dot_type
dot = x1->dot(*x2);
272 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 275 dot_type local_val(0);
276 for (
size_t i=0; i<num_my_row; ++i) {
279 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
281 local_val += 0.12345 * v * v;
287 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
288 Teuchos::outArg(
val));
294 for (
size_t i=0; i<num_my_row; ++i) {
297 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
299 local_val.fastAccessCoeff(
j) += 0.12345 * v * v;
305 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
306 Teuchos::outArg(
val));
310 out <<
"dot = " <<
dot <<
" expected = " <<
val << std::endl;
312 BaseScalar tol = 1.0e-14;
313 TEST_FLOATING_EQUALITY(
dot,
val, tol );
324 using Teuchos::ArrayView;
325 using Teuchos::Array;
326 using Teuchos::ArrayRCP;
332 typedef Teuchos::Comm<int> Tpetra_Comm;
333 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
337 if (!Kokkos::HostSpace::execution_space::is_initialized())
338 Kokkos::HostSpace::execution_space::initialize();
339 if (!Device::is_initialized())
340 Device::initialize();
343 RCP<const Tpetra_Comm> comm =
344 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
348 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
349 RCP<const Tpetra_Map> map =
350 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
352 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
353 const size_t num_my_row = myGIDs.size();
357 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
358 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
359 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
360 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
362 for (
size_t i=0; i<num_my_row; ++i) {
364 for (
size_t j=0;
j<ncol; ++
j) {
367 generate_multi_vector_coefficient<BaseScalar,size_t>(
369 val1.fastAccessCoeff(k) = v;
370 val2.fastAccessCoeff(k) = 0.12345 * v;
372 x1_view[
j][i] = val1;
373 x2_view[
j][i] = val2;
376 x1_view = Teuchos::null;
377 x2_view = Teuchos::null;
382 RCP<Tpetra_MultiVector>
y = Tpetra::createMultiVector<Scalar>(map, ncol);
383 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
389 ArrayRCP< ArrayRCP<Scalar> > y_view =
y->get2dViewNonConst();
391 BaseScalar tol = 1.0e-14;
392 for (
size_t i=0; i<num_my_row; ++i) {
394 for (
size_t j=0;
j<ncol; ++
j) {
396 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
398 val.fastAccessCoeff(k) = alpha.coeff(k)*v + 0.12345*beta.coeff(k)*v;
402 TEST_FLOATING_EQUALITY( y_view[
j][i].fastAccessCoeff(k),
403 val.fastAccessCoeff(k), tol );
416 using Teuchos::ArrayView;
417 using Teuchos::Array;
418 using Teuchos::ArrayRCP;
424 typedef Teuchos::Comm<int> Tpetra_Comm;
425 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
427 typedef typename Tpetra_MultiVector::dot_type dot_type;
430 if (!Kokkos::HostSpace::execution_space::is_initialized())
431 Kokkos::HostSpace::execution_space::initialize();
432 if (!Device::is_initialized())
433 Device::initialize();
436 RCP<const Tpetra_Comm> comm =
437 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
441 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
442 RCP<const Tpetra_Map> map =
443 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
445 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
446 const size_t num_my_row = myGIDs.size();
450 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
451 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
452 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
453 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
455 for (
size_t i=0; i<num_my_row; ++i) {
457 for (
size_t j=0;
j<ncol; ++
j) {
460 generate_multi_vector_coefficient<BaseScalar,size_t>(
462 val1.fastAccessCoeff(k) = v;
463 val2.fastAccessCoeff(k) = 0.12345 * v;
465 x1_view[
j][i] = val1;
466 x2_view[
j][i] = val2;
469 x1_view = Teuchos::null;
470 x2_view = Teuchos::null;
473 Array<dot_type> dots(ncol);
474 x1->dot(*x2, dots());
478 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 481 Array<dot_type> local_vals(ncol, dot_type(0));
482 for (
size_t i=0; i<num_my_row; ++i) {
484 for (
size_t j=0;
j<ncol; ++
j) {
486 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
488 local_vals[
j] += 0.12345 * v * v;
494 Array<dot_type> vals(ncol, dot_type(0));
495 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
496 local_vals.getRawPtr(), vals.getRawPtr());
501 Array<dot_type> local_vals(ncol, dot_type(
VectorSize, 0.0));
502 for (
size_t i=0; i<num_my_row; ++i) {
504 for (
size_t j=0;
j<ncol; ++
j) {
506 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
508 local_vals[
j].fastAccessCoeff(k) += 0.12345 * v * v;
514 Array<dot_type> vals(ncol, dot_type(
VectorSize, 0.0));
515 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
516 local_vals.getRawPtr(), vals.getRawPtr());
520 BaseScalar tol = 1.0e-14;
521 for (
size_t j=0;
j<ncol; ++
j) {
522 out <<
"dots(" <<
j <<
") = " << dots[
j]
523 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
524 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
536 using Teuchos::ArrayView;
537 using Teuchos::Array;
538 using Teuchos::ArrayRCP;
544 typedef Teuchos::Comm<int> Tpetra_Comm;
545 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
547 typedef typename Tpetra_MultiVector::dot_type dot_type;
550 if (!Kokkos::HostSpace::execution_space::is_initialized())
551 Kokkos::HostSpace::execution_space::initialize();
552 if (!Device::is_initialized())
553 Device::initialize();
556 RCP<const Tpetra_Comm> comm =
557 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
561 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
562 RCP<const Tpetra_Map> map =
563 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
565 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
566 const size_t num_my_row = myGIDs.size();
570 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
571 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
572 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
573 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
575 for (
size_t i=0; i<num_my_row; ++i) {
577 for (
size_t j=0;
j<ncol; ++
j) {
580 generate_multi_vector_coefficient<BaseScalar,size_t>(
582 val1.fastAccessCoeff(k) = v;
583 val2.fastAccessCoeff(k) = 0.12345 * v;
585 x1_view[
j][i] = val1;
586 x2_view[
j][i] = val2;
589 x1_view = Teuchos::null;
590 x2_view = Teuchos::null;
594 Teuchos::Array<size_t> cols(ncol_sub);
595 cols[0] = 4; cols[1] = 2;
596 RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
597 RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
600 Array<dot_type> dots(ncol_sub);
601 x1_sub->dot(*x2_sub, dots());
605 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 608 Array<dot_type> local_vals(ncol_sub, dot_type(0));
609 for (
size_t i=0; i<num_my_row; ++i) {
611 for (
size_t j=0;
j<ncol_sub; ++
j) {
613 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
615 local_vals[
j] += 0.12345 * v * v;
621 Array<dot_type> vals(ncol_sub, dot_type(0));
622 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
623 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
629 Array<dot_type> local_vals(ncol_sub, dot_type(
VectorSize, 0.0));
630 for (
size_t i=0; i<num_my_row; ++i) {
632 for (
size_t j=0;
j<ncol_sub; ++
j) {
634 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
636 local_vals[
j].fastAccessCoeff(k) += 0.12345 * v * v;
642 Array<dot_type> vals(ncol_sub, dot_type(
VectorSize, 0.0));
643 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
644 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
649 BaseScalar tol = 1.0e-14;
650 for (
size_t j=0;
j<ncol_sub; ++
j) {
651 out <<
"dots(" <<
j <<
") = " << dots[
j]
652 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
653 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
665 using Teuchos::ArrayView;
666 using Teuchos::Array;
667 using Teuchos::ArrayRCP;
673 typedef Teuchos::Comm<int> Tpetra_Comm;
674 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
676 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
677 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
680 if (!Kokkos::HostSpace::execution_space::is_initialized())
681 Kokkos::HostSpace::execution_space::initialize();
682 if (!Device::is_initialized())
683 Device::initialize();
687 RCP<const Tpetra_Comm> comm =
688 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
689 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
690 RCP<const Tpetra_Map> map =
691 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
693 RCP<Tpetra_CrsGraph> graph =
694 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
695 Array<GlobalOrdinal> columnIndices(2);
696 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
697 const size_t num_my_row = myGIDs.size();
698 for (
size_t i=0; i<num_my_row; ++i) {
700 columnIndices[0] = row;
703 columnIndices[1] = row+1;
706 graph->insertGlobalIndices(row, columnIndices(0,ncol));
708 graph->fillComplete();
709 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
712 Array<Scalar> vals(2);
714 for (
size_t i=0; i<num_my_row; ++i) {
716 columnIndices[0] = row;
718 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
724 columnIndices[1] = row+1;
726 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
731 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
733 matrix->fillComplete();
736 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
737 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
738 for (
size_t i=0; i<num_my_row; ++i) {
741 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
745 x_view = Teuchos::null;
754 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
755 matrix->apply(*
x, *
y);
761 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
762 BaseScalar tol = 1.0e-14;
763 for (
size_t i=0; i<num_my_row; ++i) {
766 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
768 val.fastAccessCoeff(
j) = v*v;
772 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
774 val.fastAccessCoeff(
j) += v*v;
777 TEST_EQUALITY( y_view[i].size(),
VectorSize );
779 TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(
j),
val.fastAccessCoeff(
j), tol );
791 using Teuchos::ArrayView;
792 using Teuchos::Array;
793 using Teuchos::ArrayRCP;
799 typedef Teuchos::Comm<int> Tpetra_Comm;
800 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
802 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
803 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
806 if (!Kokkos::HostSpace::execution_space::is_initialized())
807 Kokkos::HostSpace::execution_space::initialize();
808 if (!Device::is_initialized())
809 Device::initialize();
813 RCP<const Tpetra_Comm> comm =
814 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
815 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
816 RCP<const Tpetra_Map> map =
817 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
819 RCP<Tpetra_CrsGraph> graph =
820 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
821 Array<GlobalOrdinal> columnIndices(2);
822 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
823 const size_t num_my_row = myGIDs.size();
824 for (
size_t i=0; i<num_my_row; ++i) {
826 columnIndices[0] = row;
829 columnIndices[1] = row+1;
832 graph->insertGlobalIndices(row, columnIndices(0,ncol));
834 graph->fillComplete();
835 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
838 Array<Scalar> vals(2);
840 for (
size_t i=0; i<num_my_row; ++i) {
842 columnIndices[0] = row;
844 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
850 columnIndices[1] = row+1;
852 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
857 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
859 matrix->fillComplete();
863 RCP<Tpetra_MultiVector>
x = Tpetra::createMultiVector<Scalar>(map, ncol);
864 ArrayRCP< ArrayRCP<Scalar> > x_view =
x->get2dViewNonConst();
865 for (
size_t i=0; i<num_my_row; ++i) {
867 for (
size_t j=0;
j<ncol; ++
j) {
870 generate_multi_vector_coefficient<BaseScalar,size_t>(
872 val.fastAccessCoeff(k) = v;
877 x_view = Teuchos::null;
886 RCP<Tpetra_MultiVector>
y = Tpetra::createMultiVector<Scalar>(map, ncol);
887 matrix->apply(*
x, *
y);
893 ArrayRCP< ArrayRCP<Scalar> > y_view =
y->get2dViewNonConst();
894 BaseScalar tol = 1.0e-14;
895 for (
size_t i=0; i<num_my_row; ++i) {
897 for (
size_t j=0;
j<ncol; ++
j) {
899 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
901 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
903 val.fastAccessCoeff(k) = v1*v2;
907 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
909 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
911 val.fastAccessCoeff(k) += v1*v2;
916 TEST_FLOATING_EQUALITY( y_view[
j][i].fastAccessCoeff(k),
917 val.fastAccessCoeff(k), tol );
930 using Teuchos::ArrayView;
931 using Teuchos::Array;
932 using Teuchos::ArrayRCP;
938 typedef Teuchos::Comm<int> Tpetra_Comm;
939 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
941 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
942 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
945 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
948 if (!Kokkos::HostSpace::execution_space::is_initialized())
949 Kokkos::HostSpace::execution_space::initialize();
950 if (!Device::is_initialized())
951 Device::initialize();
955 RCP<const Tpetra_Comm> comm =
956 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
957 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
958 RCP<const Tpetra_Map> map =
959 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
961 RCP<Tpetra_CrsGraph> graph =
962 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
963 Array<GlobalOrdinal> columnIndices(2);
964 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
965 const size_t num_my_row = myGIDs.size();
966 for (
size_t i=0; i<num_my_row; ++i) {
968 columnIndices[0] = row;
971 columnIndices[1] = row+1;
974 graph->insertGlobalIndices(row, columnIndices(0,ncol));
976 graph->fillComplete();
977 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
980 Array<Scalar> vals(2);
982 for (
size_t i=0; i<num_my_row; ++i) {
984 columnIndices[0] = row;
986 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
992 columnIndices[1] = row+1;
994 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
999 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1001 matrix->fillComplete();
1004 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1005 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1006 for (
size_t i=0; i<num_my_row; ++i) {
1009 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
1015 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
1016 matrix->apply(*
x, *
y);
1019 RCP<const Tpetra_Map> flat_x_map, flat_y_map;
1020 RCP<const Tpetra_CrsGraph> flat_graph =
1022 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1026 RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
1027 RCP<Flat_Tpetra_Vector> flat_x =
1029 RCP<Flat_Tpetra_Vector> flat_y =
1031 flat_matrix->apply(*flat_x, *flat_y);
1037 BaseScalar tol = 1.0e-14;
1038 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
1039 ArrayRCP<Scalar> y2_view = y2->get1dViewNonConst();
1040 for (
size_t i=0; i<num_my_row; ++i) {
1041 TEST_EQUALITY( y_view[i].size(),
VectorSize );
1042 TEST_EQUALITY( y2_view[i].size(),
VectorSize );
1044 TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(
j),
1045 y2_view[i].fastAccessCoeff(
j), tol );
1057 using Teuchos::ArrayView;
1058 using Teuchos::Array;
1059 using Teuchos::ArrayRCP;
1060 using Teuchos::ParameterList;
1066 typedef Teuchos::Comm<int> Tpetra_Comm;
1067 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1069 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1070 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1073 if (!Kokkos::HostSpace::execution_space::is_initialized())
1074 Kokkos::HostSpace::execution_space::initialize();
1075 if (!Device::is_initialized())
1076 Device::initialize();
1080 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1081 RCP<const Tpetra_Comm> comm =
1082 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1083 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1084 RCP<const Tpetra_Map> map =
1085 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1087 RCP<Tpetra_CrsGraph> graph =
1088 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1089 Array<GlobalOrdinal> columnIndices(3);
1090 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1091 const size_t num_my_row = myGIDs.size();
1092 for (
size_t i=0; i<num_my_row; ++i) {
1094 if (row == 0 || row == nrow-1) {
1095 columnIndices[0] = row;
1096 graph->insertGlobalIndices(row, columnIndices(0,1));
1099 columnIndices[0] = row-1;
1100 columnIndices[1] = row;
1101 columnIndices[2] = row+1;
1102 graph->insertGlobalIndices(row, columnIndices(0,3));
1105 graph->fillComplete();
1106 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1109 Array<Scalar> vals(3);
1112 a_val.fastAccessCoeff(
j) =
1113 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1115 for (
size_t i=0; i<num_my_row; ++i) {
1117 if (row == 0 || row == nrow-1) {
1118 columnIndices[0] = row;
1120 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1123 columnIndices[0] = row-1;
1124 columnIndices[1] = row;
1125 columnIndices[2] = row+1;
1126 vals[0] =
Scalar(-1.0) * a_val;
1127 vals[1] =
Scalar(2.0) * a_val;
1128 vals[2] =
Scalar(-1.0) * a_val;
1129 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1132 matrix->fillComplete();
1134 matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
1135 Teuchos::VERB_EXTREME);
1138 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1139 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1142 b_val.fastAccessCoeff(
j) =
1143 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1145 for (
size_t i=0; i<num_my_row; ++i) {
1147 if (row == 0 || row == nrow-1)
1150 b_view[i] = -
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 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
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].size(),
VectorSize );
1185 v.fastAccessCoeff(
j) = BaseScalar(0.0);
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;
1216 typedef Teuchos::Comm<int> Tpetra_Comm;
1217 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1219 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1220 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1223 if (!Kokkos::HostSpace::execution_space::is_initialized())
1224 Kokkos::HostSpace::execution_space::initialize();
1225 if (!Device::is_initialized())
1226 Device::initialize();
1230 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1231 RCP<const Tpetra_Comm> comm =
1232 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1233 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1234 RCP<const Tpetra_Map> map =
1235 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1237 RCP<Tpetra_CrsGraph> graph =
1238 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1239 Array<GlobalOrdinal> columnIndices(3);
1240 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1241 const size_t num_my_row = myGIDs.size();
1242 for (
size_t i=0; i<num_my_row; ++i) {
1244 if (row == 0 || row == nrow-1) {
1245 columnIndices[0] = row;
1246 graph->insertGlobalIndices(row, columnIndices(0,1));
1249 columnIndices[0] = row-1;
1250 columnIndices[1] = row;
1251 columnIndices[2] = row+1;
1252 graph->insertGlobalIndices(row, columnIndices(0,3));
1255 graph->fillComplete();
1256 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1259 Array<Scalar> vals(3);
1262 a_val.fastAccessCoeff(
j) =
1263 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1265 for (
size_t i=0; i<num_my_row; ++i) {
1267 if (row == 0 || row == nrow-1) {
1268 columnIndices[0] = row;
1270 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1273 columnIndices[0] = row-1;
1274 columnIndices[1] = row;
1275 columnIndices[2] = row+1;
1276 vals[0] =
Scalar(-1.0) * a_val;
1277 vals[1] =
Scalar(2.0) * a_val;
1278 vals[2] =
Scalar(-1.0) * a_val;
1279 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1282 matrix->fillComplete();
1285 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1286 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1289 b_val.fastAccessCoeff(
j) =
1290 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1292 for (
size_t i=0; i<num_my_row; ++i) {
1294 if (row == 0 || row == nrow-1)
1297 b_view[i] = -
Scalar(b_val * h * h);
1301 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1302 RCP<OP> matrix_op = matrix;
1303 RCP<ParameterList> muelu_params =
1304 getParametersFromXmlFile(
"muelu_cheby.xml");
1306 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1309 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1310 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1311 typedef typename BST::mag_type base_mag_type;
1312 typedef typename Tpetra_Vector::mag_type mag_type;
1313 base_mag_type btol = 1e-9;
1314 mag_type tol = btol;
1317 out.getOStream().get());
1318 TEST_EQUALITY_CONST( solved,
true );
1325 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1327 for (
size_t i=0; i<num_my_row; ++i) {
1329 BaseScalar xx = row * h;
1331 val.fastAccessCoeff(
j) =
1332 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1334 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1339 if (BST::magnitude(v.coeff(
j)) < btol)
1340 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1341 if (BST::magnitude(
val.coeff(
j)) < btol)
1342 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1346 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), btol);
1358 #if defined(HAVE_STOKHOS_BELOS) 1368 using Teuchos::ArrayView;
1369 using Teuchos::Array;
1370 using Teuchos::ArrayRCP;
1371 using Teuchos::ParameterList;
1377 typedef Teuchos::Comm<int> Tpetra_Comm;
1378 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1380 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1381 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1384 if (!Kokkos::HostSpace::execution_space::is_initialized())
1385 Kokkos::HostSpace::execution_space::initialize();
1386 if (!Device::is_initialized())
1387 Device::initialize();
1391 RCP<const Tpetra_Comm> comm =
1392 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1393 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1394 RCP<const Tpetra_Map> map =
1395 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1397 RCP<Tpetra_CrsGraph> graph =
1398 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1399 Array<GlobalOrdinal> columnIndices(2);
1400 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1401 const size_t num_my_row = myGIDs.size();
1402 for (
size_t i=0; i<num_my_row; ++i) {
1404 columnIndices[0] = row;
1406 if (row != nrow-1) {
1407 columnIndices[1] = row+1;
1410 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1412 graph->fillComplete();
1413 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1416 Array<Scalar> vals(2);
1418 for (
size_t i=0; i<num_my_row; ++i) {
1420 columnIndices[0] = row;
1422 val.fastAccessCoeff(
j) =
j+1;
1426 if (row != nrow-1) {
1427 columnIndices[1] = row+1;
1429 val.fastAccessCoeff(
j) =
j+1;
1433 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1435 matrix->fillComplete();
1438 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1439 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1440 for (
size_t i=0; i<num_my_row; ++i) {
1445 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1446 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1447 typedef BaseScalar BelosScalar;
1449 typedef Scalar BelosScalar;
1452 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1453 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1454 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1455 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1456 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1457 typename ST::magnitudeType tol = 1e-9;
1458 belosParams->set(
"Flexible Gmres",
false);
1459 belosParams->set(
"Num Blocks", 100);
1460 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1461 belosParams->set(
"Maximum Iterations", 100);
1462 belosParams->set(
"Verbosity", 33);
1463 belosParams->set(
"Output Style", 1);
1464 belosParams->set(
"Output Frequency", 1);
1465 belosParams->set(
"Output Stream", out.getOStream());
1466 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1467 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1468 problem->setProblem();
1469 Belos::ReturnType ret = solver->solve();
1470 TEST_EQUALITY_CONST( ret, Belos::Converged );
1482 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1483 for (
size_t i=0; i<num_my_row; ++i) {
1487 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
1492 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1497 if (ST::magnitude(v.coeff(
j)) < tol)
1498 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1502 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1514 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) 1525 using Teuchos::ArrayView;
1526 using Teuchos::Array;
1527 using Teuchos::ArrayRCP;
1528 using Teuchos::ParameterList;
1534 typedef Teuchos::Comm<int> Tpetra_Comm;
1535 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1537 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1538 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1541 if (!Kokkos::HostSpace::execution_space::is_initialized())
1542 Kokkos::HostSpace::execution_space::initialize();
1543 if (!Device::is_initialized())
1544 Device::initialize();
1548 RCP<const Tpetra_Comm> comm =
1549 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1550 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1551 RCP<const Tpetra_Map> map =
1552 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1554 RCP<Tpetra_CrsGraph> graph =
1555 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1556 Array<GlobalOrdinal> columnIndices(2);
1557 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1558 const size_t num_my_row = myGIDs.size();
1559 for (
size_t i=0; i<num_my_row; ++i) {
1561 columnIndices[0] = row;
1563 if (row != nrow-1) {
1564 columnIndices[1] = row+1;
1567 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1569 graph->fillComplete();
1570 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1573 Array<Scalar> vals(2);
1575 for (
size_t i=0; i<num_my_row; ++i) {
1577 columnIndices[0] = row;
1579 val.fastAccessCoeff(
j) =
j+1;
1583 if (row != nrow-1) {
1584 columnIndices[1] = row+1;
1586 val.fastAccessCoeff(
j) =
j+1;
1590 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1592 matrix->fillComplete();
1595 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1596 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1597 for (
size_t i=0; i<num_my_row; ++i) {
1602 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
1603 Ifpack2::Factory factory;
1604 RCP<Prec> M = factory.create<Tpetra_CrsMatrix>(
"RILUK", matrix);
1609 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1610 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1611 typedef BaseScalar BelosScalar;
1613 typedef Scalar BelosScalar;
1616 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1617 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1618 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1619 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1620 problem->setRightPrec(M);
1621 problem->setProblem();
1622 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1623 typename ST::magnitudeType tol = 1e-9;
1624 belosParams->set(
"Flexible Gmres",
false);
1625 belosParams->set(
"Num Blocks", 100);
1626 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1627 belosParams->set(
"Maximum Iterations", 100);
1628 belosParams->set(
"Verbosity", 33);
1629 belosParams->set(
"Output Style", 1);
1630 belosParams->set(
"Output Frequency", 1);
1631 belosParams->set(
"Output Stream", out.getOStream());
1633 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1634 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1635 Belos::ReturnType ret = solver->solve();
1636 TEST_EQUALITY_CONST( ret, Belos::Converged );
1648 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1649 for (
size_t i=0; i<num_my_row; ++i) {
1653 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
1658 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1663 if (ST::magnitude(v.coeff(
j)) < tol)
1664 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1668 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1680 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) 1690 using Teuchos::ArrayView;
1691 using Teuchos::Array;
1692 using Teuchos::ArrayRCP;
1693 using Teuchos::ParameterList;
1694 using Teuchos::getParametersFromXmlFile;
1700 typedef Teuchos::Comm<int> Tpetra_Comm;
1701 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1703 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1704 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1707 if (!Kokkos::HostSpace::execution_space::is_initialized())
1708 Kokkos::HostSpace::execution_space::initialize();
1709 if (!Device::is_initialized())
1710 Device::initialize();
1714 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1715 RCP<const Tpetra_Comm> comm =
1716 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1717 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1718 RCP<const Tpetra_Map> map =
1719 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1721 RCP<Tpetra_CrsGraph> graph =
1722 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1723 Array<GlobalOrdinal> columnIndices(3);
1724 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1725 const size_t num_my_row = myGIDs.size();
1726 for (
size_t i=0; i<num_my_row; ++i) {
1728 if (row == 0 || row == nrow-1) {
1729 columnIndices[0] = row;
1730 graph->insertGlobalIndices(row, columnIndices(0,1));
1733 columnIndices[0] = row-1;
1734 columnIndices[1] = row;
1735 columnIndices[2] = row+1;
1736 graph->insertGlobalIndices(row, columnIndices(0,3));
1739 graph->fillComplete();
1740 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1743 Array<Scalar> vals(3);
1746 a_val.fastAccessCoeff(
j) =
1747 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1749 for (
size_t i=0; i<num_my_row; ++i) {
1751 if (row == 0 || row == nrow-1) {
1752 columnIndices[0] = row;
1754 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1757 columnIndices[0] = row-1;
1758 columnIndices[1] = row;
1759 columnIndices[2] = row+1;
1760 vals[0] =
Scalar(-1.0) * a_val;
1761 vals[1] =
Scalar(2.0) * a_val;
1762 vals[2] =
Scalar(-1.0) * a_val;
1763 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1766 matrix->fillComplete();
1769 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1770 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1773 b_val.fastAccessCoeff(
j) =
1774 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1776 for (
size_t i=0; i<num_my_row; ++i) {
1778 if (row == 0 || row == nrow-1)
1781 b_view[i] = -
Scalar(b_val * h * h);
1785 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1786 RCP<ParameterList> muelu_params =
1787 getParametersFromXmlFile(
"muelu_cheby.xml");
1788 RCP<OP> matrix_op = matrix;
1790 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1793 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1794 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1795 typedef BaseScalar BelosScalar;
1797 typedef Scalar BelosScalar;
1800 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1801 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1802 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1803 problem->setRightPrec(M);
1804 problem->setProblem();
1805 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1806 typename ST::magnitudeType tol = 1e-9;
1807 belosParams->set(
"Num Blocks", 100);
1808 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1809 belosParams->set(
"Maximum Iterations", 100);
1810 belosParams->set(
"Verbosity", 33);
1811 belosParams->set(
"Output Style", 1);
1812 belosParams->set(
"Output Frequency", 1);
1813 belosParams->set(
"Output Stream", out.getOStream());
1816 belosParams->set(
"Implicit Residual Scaling",
"None");
1818 RCP<Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true> > solver =
1819 rcp(
new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true>(problem, belosParams));
1820 Belos::ReturnType ret = solver->solve();
1821 TEST_EQUALITY_CONST( ret, Belos::Converged );
1823 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT 1825 std::vector<int> ensemble_iterations =
1826 solver->getResidualStatusTest()->getEnsembleIterations();
1827 out <<
"Ensemble iterations = ";
1829 out << ensemble_iterations[i] <<
" ";
1838 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1840 for (
size_t i=0; i<num_my_row; ++i) {
1842 BaseScalar xx = row * h;
1844 val.fastAccessCoeff(
j) =
1845 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1847 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1852 if (ST::magnitude(v.coeff(
j)) < tol)
1853 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1854 if (ST::magnitude(
val.coeff(
j)) < tol)
1855 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1859 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1872 #if defined(HAVE_STOKHOS_AMESOS2) 1882 using Teuchos::ArrayView;
1883 using Teuchos::Array;
1884 using Teuchos::ArrayRCP;
1885 using Teuchos::ParameterList;
1891 typedef Teuchos::Comm<int> Tpetra_Comm;
1892 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1895 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1896 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1899 if (!Kokkos::HostSpace::execution_space::is_initialized())
1900 Kokkos::HostSpace::execution_space::initialize();
1901 if (!Device::is_initialized())
1902 Device::initialize();
1906 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1907 RCP<const Tpetra_Comm> comm =
1908 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1909 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1910 RCP<const Tpetra_Map> map =
1911 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1913 RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map,
size_t(3));
1914 Array<GlobalOrdinal> columnIndices(3);
1915 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1916 const size_t num_my_row = myGIDs.size();
1917 for (
size_t i=0; i<num_my_row; ++i) {
1919 if (row == 0 || row == nrow-1) {
1920 columnIndices[0] = row;
1921 graph->insertGlobalIndices(row, columnIndices(0,1));
1924 columnIndices[0] = row-1;
1925 columnIndices[1] = row;
1926 columnIndices[2] = row+1;
1927 graph->insertGlobalIndices(row, columnIndices(0,3));
1930 graph->fillComplete();
1931 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1934 Array<Scalar> vals(3);
1937 a_val.fastAccessCoeff(
j) =
1938 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1940 for (
size_t i=0; i<num_my_row; ++i) {
1942 if (row == 0 || row == nrow-1) {
1943 columnIndices[0] = row;
1945 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1948 columnIndices[0] = row-1;
1949 columnIndices[1] = row;
1950 columnIndices[2] = row+1;
1951 vals[0] =
Scalar(-1.0) * a_val;
1952 vals[1] =
Scalar(2.0) * a_val;
1953 vals[2] =
Scalar(-1.0) * a_val;
1954 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1957 matrix->fillComplete();
1960 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1961 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1964 b_val.fastAccessCoeff(
j) =
1965 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1967 for (
size_t i=0; i<num_my_row; ++i) {
1969 if (row == 0 || row == nrow-1)
1972 b_view[i] = -
Scalar(b_val * h * h);
1976 typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
1977 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1978 std::string solver_name;
1979 #if defined(HAVE_AMESOS2_BASKER) 1980 solver_name =
"basker";
1981 #elif defined(HAVE_AMESOS2_KLU2) 1982 solver_name =
"klu2";
1983 #elif defined(HAVE_AMESOS2_SUPERLUDIST) 1984 solver_name =
"superlu_dist";
1985 #elif defined(HAVE_AMESOS2_SUPERLUMT) 1986 solver_name =
"superlu_mt";
1987 #elif defined(HAVE_AMESOS2_SUPERLU) 1988 solver_name =
"superlu";
1989 #elif defined(HAVE_AMESOS2_PARDISO_MKL) 1990 solver_name =
"pardisomkl";
1991 #elif defined(HAVE_AMESOS2_LAPACK) 1992 solver_name =
"lapack";
1993 #elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL) 1994 solver_name =
"lapack";
2000 out <<
"Solving linear system with " << solver_name << std::endl;
2001 RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2002 solver_name, matrix,
x, b);
2009 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2010 typename ST::magnitudeType tol = 1e-9;
2011 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
2013 for (
size_t i=0; i<num_my_row; ++i) {
2015 BaseScalar xx = row * h;
2017 val.fastAccessCoeff(
j) =
2018 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
2020 TEST_EQUALITY( x_view[i].size(),
VectorSize );
2025 if (ST::magnitude(v.coeff(
j)) < tol)
2026 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2027 if (ST::magnitude(
val.coeff(
j)) < tol)
2028 val.fastAccessCoeff(
j) = BaseScalar(0.0);
2032 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
2044 #define CRSMATRIX_MP_VECTOR_TESTS_SLGN(S, LO, GO, N) \ 2045 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorAdd, S, LO, GO, N ) \ 2046 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorDot, S, LO, GO, N ) \ 2047 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorAdd, S, LO, GO, N ) \ 2048 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDot, S, LO, GO, N ) \ 2049 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDotSub, S, LO, GO, N ) \ 2050 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixVectorMultiply, S, LO, GO, N ) \ 2051 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, S, LO, GO, N ) \ 2052 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Flatten, S, LO, GO, N ) \ 2053 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimpleCG, S, LO, GO, N ) \ 2054 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimplePCG_Muelu, S, LO, GO, N ) \ 2055 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES, S, LO, GO, N ) \ 2056 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, S, LO, GO, N ) \ 2057 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosCG_Muelu, S, LO, GO, N ) \ 2058 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2, S, LO, GO, N ) 2060 #define CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N) \ 2061 typedef Stokhos::DeviceForNode<N>::type Device; \ 2062 typedef Stokhos::StaticFixedStorage<int,double,VectorSize,Device::execution_space> SFS; \ 2063 CRSMATRIX_MP_VECTOR_TESTS_SLGN(SFS, int, int, N) 2065 #define CRSMATRIX_MP_VECTOR_TESTS_N(N) \ 2066 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)
Kokkos::DefaultExecutionSpace execution_space
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)