42 #include "Teuchos_UnitTestHelpers.hpp" 46 #include "Teuchos_XMLParameterListCoreHelpers.hpp" 51 #include "Tpetra_Core.hpp" 52 #include "Tpetra_Map.hpp" 53 #include "Tpetra_MultiVector.hpp" 54 #include "Tpetra_Vector.hpp" 55 #include "Tpetra_CrsGraph.hpp" 56 #include "Tpetra_CrsMatrix.hpp" 60 #ifdef HAVE_STOKHOS_BELOS 62 #include "BelosLinearProblem.hpp" 63 #include "BelosPseudoBlockGmresSolMgr.hpp" 64 #include "BelosPseudoBlockCGSolMgr.hpp" 68 #ifdef HAVE_STOKHOS_IFPACK2 70 #include "Ifpack2_Factory.hpp" 74 #ifdef HAVE_STOKHOS_MUELU 76 #include "MueLu_CreateTpetraPreconditioner.hpp" 80 #ifdef HAVE_STOKHOS_AMESOS2 82 #include "Amesos2_Factory.hpp" 85 template <
typename scalar,
typename ordinal>
89 const ordinal iColFEM,
90 const ordinal iStoch )
92 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
93 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
94 return X_fem + X_stoch;
98 template <
typename scalar,
typename ordinal>
102 const ordinal nStoch,
103 const ordinal iColFEM,
105 const ordinal iStoch)
107 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
108 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
109 const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
110 return X_fem + X_stoch + X_col;
119 #if defined(__CUDACC__) 133 using Teuchos::ArrayView;
134 using Teuchos::Array;
135 using Teuchos::ArrayRCP;
140 typedef Teuchos::Comm<int> Tpetra_Comm;
141 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
142 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
145 if ( !Kokkos::is_initialized() )
146 Kokkos::initialize();
149 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
153 RCP<const Tpetra_Map> map =
154 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
156 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
157 const size_t num_my_row = myGIDs.size();
160 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
161 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
162 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
163 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
165 for (
size_t i=0; i<num_my_row; ++i) {
168 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
169 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
174 x1_view = Teuchos::null;
175 x2_view = Teuchos::null;
180 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
181 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
187 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
189 BaseScalar tol = 1.0e-14;
190 for (
size_t i=0; i<num_my_row; ++i) {
193 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
195 val.fastAccessCoeff(
j) = alpha.coeff(
j)*v + 0.12345*beta.coeff(
j)*v;
197 TEST_EQUALITY( y_view[i].size(),
VectorSize );
211 using Teuchos::ArrayView;
212 using Teuchos::Array;
213 using Teuchos::ArrayRCP;
218 typedef Teuchos::Comm<int> Tpetra_Comm;
219 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
220 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
221 typedef typename Tpetra_Vector::dot_type dot_type;
224 if ( !Kokkos::is_initialized() )
225 Kokkos::initialize();
228 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
232 RCP<const Tpetra_Map> map =
233 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
235 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
236 const size_t num_my_row = myGIDs.size();
239 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
240 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
241 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
242 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
244 for (
size_t i=0; i<num_my_row; ++i) {
247 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
248 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow,
VectorSize, row,
j);
253 x1_view = Teuchos::null;
254 x2_view = Teuchos::null;
257 dot_type
dot = x1->dot(*x2);
261 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 264 dot_type local_val(0);
265 for (
size_t i=0; i<num_my_row; ++i) {
268 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
270 local_val += 0.12345 * v * v;
276 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
277 Teuchos::outArg(
val));
283 for (
size_t i=0; i<num_my_row; ++i) {
286 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
288 local_val.fastAccessCoeff(
j) += 0.12345 * v * v;
294 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
295 Teuchos::outArg(
val));
299 out <<
"dot = " <<
dot <<
" expected = " <<
val << std::endl;
301 BaseScalar tol = 1.0e-14;
302 TEST_FLOATING_EQUALITY(
dot,
val, tol );
313 using Teuchos::ArrayView;
314 using Teuchos::Array;
315 using Teuchos::ArrayRCP;
320 typedef Teuchos::Comm<int> Tpetra_Comm;
321 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
322 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
325 if ( !Kokkos::is_initialized() )
326 Kokkos::initialize();
329 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
333 RCP<const Tpetra_Map> map =
334 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
336 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
337 const size_t num_my_row = myGIDs.size();
341 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
342 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
343 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
344 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
346 for (
size_t i=0; i<num_my_row; ++i) {
348 for (
size_t j=0;
j<ncol; ++
j) {
351 generate_multi_vector_coefficient<BaseScalar,size_t>(
353 val1.fastAccessCoeff(k) = v;
354 val2.fastAccessCoeff(k) = 0.12345 * v;
356 x1_view[
j][i] = val1;
357 x2_view[
j][i] = val2;
360 x1_view = Teuchos::null;
361 x2_view = Teuchos::null;
366 RCP<Tpetra_MultiVector>
y = Tpetra::createMultiVector<Scalar>(map, ncol);
367 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
373 ArrayRCP< ArrayRCP<Scalar> > y_view =
y->get2dViewNonConst();
375 BaseScalar tol = 1.0e-14;
376 for (
size_t i=0; i<num_my_row; ++i) {
378 for (
size_t j=0;
j<ncol; ++
j) {
380 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
382 val.fastAccessCoeff(k) = alpha.coeff(k)*v + 0.12345*beta.coeff(k)*v;
387 val.fastAccessCoeff(k), tol );
400 using Teuchos::ArrayView;
401 using Teuchos::Array;
402 using Teuchos::ArrayRCP;
407 typedef Teuchos::Comm<int> Tpetra_Comm;
408 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
409 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
410 typedef typename Tpetra_MultiVector::dot_type dot_type;
413 if ( !Kokkos::is_initialized() )
414 Kokkos::initialize();
417 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
421 RCP<const Tpetra_Map> map =
422 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
424 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
425 const size_t num_my_row = myGIDs.size();
429 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
430 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
431 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
432 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
434 for (
size_t i=0; i<num_my_row; ++i) {
436 for (
size_t j=0;
j<ncol; ++
j) {
439 generate_multi_vector_coefficient<BaseScalar,size_t>(
441 val1.fastAccessCoeff(k) = v;
442 val2.fastAccessCoeff(k) = 0.12345 * v;
444 x1_view[
j][i] = val1;
445 x2_view[
j][i] = val2;
448 x1_view = Teuchos::null;
449 x2_view = Teuchos::null;
452 Array<dot_type> dots(ncol);
453 x1->dot(*x2, dots());
457 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 460 Array<dot_type> local_vals(ncol, dot_type(0));
461 for (
size_t i=0; i<num_my_row; ++i) {
463 for (
size_t j=0;
j<ncol; ++
j) {
465 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
467 local_vals[
j] += 0.12345 * v * v;
473 Array<dot_type> vals(ncol, dot_type(0));
474 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
475 local_vals.getRawPtr(), vals.getRawPtr());
480 Array<dot_type> local_vals(ncol, dot_type(
VectorSize, 0.0));
481 for (
size_t i=0; i<num_my_row; ++i) {
483 for (
size_t j=0;
j<ncol; ++
j) {
485 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
487 local_vals[
j].fastAccessCoeff(k) += 0.12345 * v * v;
493 Array<dot_type> vals(ncol, dot_type(
VectorSize, 0.0));
494 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
495 local_vals.getRawPtr(), vals.getRawPtr());
499 BaseScalar tol = 1.0e-14;
500 for (
size_t j=0;
j<ncol; ++
j) {
501 out <<
"dots(" <<
j <<
") = " << dots[
j]
502 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
503 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
515 using Teuchos::ArrayView;
516 using Teuchos::Array;
517 using Teuchos::ArrayRCP;
522 typedef Teuchos::Comm<int> Tpetra_Comm;
523 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
524 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
525 typedef typename Tpetra_MultiVector::dot_type dot_type;
528 if ( !Kokkos::is_initialized() )
529 Kokkos::initialize();
532 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
536 RCP<const Tpetra_Map> map =
537 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
539 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
540 const size_t num_my_row = myGIDs.size();
544 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
545 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
546 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
547 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
549 for (
size_t i=0; i<num_my_row; ++i) {
551 for (
size_t j=0;
j<ncol; ++
j) {
554 generate_multi_vector_coefficient<BaseScalar,size_t>(
556 val1.fastAccessCoeff(k) = v;
557 val2.fastAccessCoeff(k) = 0.12345 * v;
559 x1_view[
j][i] = val1;
560 x2_view[
j][i] = val2;
563 x1_view = Teuchos::null;
564 x2_view = Teuchos::null;
568 Teuchos::Array<size_t> cols(ncol_sub);
569 cols[0] = 4; cols[1] = 2;
570 RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
571 RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
574 Array<dot_type> dots(ncol_sub);
575 x1_sub->dot(*x2_sub, dots());
579 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 582 Array<dot_type> local_vals(ncol_sub, dot_type(0));
583 for (
size_t i=0; i<num_my_row; ++i) {
585 for (
size_t j=0;
j<ncol_sub; ++
j) {
587 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
589 local_vals[
j] += 0.12345 * v * v;
595 Array<dot_type> vals(ncol_sub, dot_type(0));
596 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
597 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
603 Array<dot_type> local_vals(ncol_sub, dot_type(
VectorSize, 0.0));
604 for (
size_t i=0; i<num_my_row; ++i) {
606 for (
size_t j=0;
j<ncol_sub; ++
j) {
608 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
610 local_vals[
j].fastAccessCoeff(k) += 0.12345 * v * v;
616 Array<dot_type> vals(ncol_sub, dot_type(
VectorSize, 0.0));
617 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
618 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
623 BaseScalar tol = 1.0e-14;
624 for (
size_t j=0;
j<ncol_sub; ++
j) {
625 out <<
"dots(" <<
j <<
") = " << dots[
j]
626 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
627 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
639 using Teuchos::ArrayView;
640 using Teuchos::Array;
641 using Teuchos::ArrayRCP;
646 typedef Teuchos::Comm<int> Tpetra_Comm;
647 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
648 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
649 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
650 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
653 if ( !Kokkos::is_initialized() )
654 Kokkos::initialize();
658 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
659 RCP<const Tpetra_Map> map =
660 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
662 RCP<Tpetra_CrsGraph> graph =
663 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
664 Array<GlobalOrdinal> columnIndices(2);
665 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
666 const size_t num_my_row = myGIDs.size();
667 for (
size_t i=0; i<num_my_row; ++i) {
669 columnIndices[0] = row;
672 columnIndices[1] = row+1;
675 graph->insertGlobalIndices(row, columnIndices(0,ncol));
677 graph->fillComplete();
678 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
681 Array<Scalar> vals(2);
683 for (
size_t i=0; i<num_my_row; ++i) {
685 columnIndices[0] = row;
687 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
693 columnIndices[1] = row+1;
695 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
700 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
702 matrix->fillComplete();
705 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
706 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
707 for (
size_t i=0; i<num_my_row; ++i) {
710 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
714 x_view = Teuchos::null;
723 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
724 matrix->apply(*
x, *
y);
730 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
731 BaseScalar tol = 1.0e-14;
732 for (
size_t i=0; i<num_my_row; ++i) {
735 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
737 val.fastAccessCoeff(
j) = v*v;
741 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
743 val.fastAccessCoeff(
j) += v*v;
746 TEST_EQUALITY( y_view[i].size(),
VectorSize );
760 using Teuchos::ArrayView;
761 using Teuchos::Array;
762 using Teuchos::ArrayRCP;
767 typedef Teuchos::Comm<int> Tpetra_Comm;
768 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
769 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
770 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
771 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
774 if ( !Kokkos::is_initialized() )
775 Kokkos::initialize();
779 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
780 RCP<const Tpetra_Map> map =
781 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
783 RCP<Tpetra_CrsGraph> graph =
784 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
785 Array<GlobalOrdinal> columnIndices(2);
786 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
787 const size_t num_my_row = myGIDs.size();
788 for (
size_t i=0; i<num_my_row; ++i) {
790 columnIndices[0] = row;
793 columnIndices[1] = row+1;
796 graph->insertGlobalIndices(row, columnIndices(0,ncol));
798 graph->fillComplete();
799 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
802 Array<Scalar> vals(2);
804 for (
size_t i=0; i<num_my_row; ++i) {
806 columnIndices[0] = row;
808 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
814 columnIndices[1] = row+1;
816 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
821 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
823 matrix->fillComplete();
827 RCP<Tpetra_MultiVector>
x = Tpetra::createMultiVector<Scalar>(map, ncol);
828 ArrayRCP< ArrayRCP<Scalar> > x_view =
x->get2dViewNonConst();
829 for (
size_t i=0; i<num_my_row; ++i) {
831 for (
size_t j=0;
j<ncol; ++
j) {
834 generate_multi_vector_coefficient<BaseScalar,size_t>(
836 val.fastAccessCoeff(k) = v;
841 x_view = Teuchos::null;
850 RCP<Tpetra_MultiVector>
y = Tpetra::createMultiVector<Scalar>(map, ncol);
851 matrix->apply(*
x, *
y);
857 ArrayRCP< ArrayRCP<Scalar> > y_view =
y->get2dViewNonConst();
858 BaseScalar tol = 1.0e-14;
859 for (
size_t i=0; i<num_my_row; ++i) {
861 for (
size_t j=0;
j<ncol; ++
j) {
863 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
865 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
867 val.fastAccessCoeff(k) = v1*v2;
871 BaseScalar v1 = generate_vector_coefficient<BaseScalar,size_t>(
873 BaseScalar v2 = generate_multi_vector_coefficient<BaseScalar,size_t>(
875 val.fastAccessCoeff(k) += v1*v2;
881 val.fastAccessCoeff(k), tol );
894 using Teuchos::ArrayView;
895 using Teuchos::Array;
896 using Teuchos::ArrayRCP;
901 typedef Teuchos::Comm<int> Tpetra_Comm;
902 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
903 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
904 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
905 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
907 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
908 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
911 if ( !Kokkos::is_initialized() )
912 Kokkos::initialize();
916 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
917 RCP<const Tpetra_Map> map =
918 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
920 RCP<Tpetra_CrsGraph> graph =
921 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
922 Array<GlobalOrdinal> columnIndices(2);
923 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
924 const size_t num_my_row = myGIDs.size();
925 for (
size_t i=0; i<num_my_row; ++i) {
927 columnIndices[0] = row;
930 columnIndices[1] = row+1;
933 graph->insertGlobalIndices(row, columnIndices(0,ncol));
935 graph->fillComplete();
936 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
939 Array<Scalar> vals(2);
941 for (
size_t i=0; i<num_my_row; ++i) {
943 columnIndices[0] = row;
945 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
951 columnIndices[1] = row+1;
953 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
958 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
960 matrix->fillComplete();
963 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
964 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
965 for (
size_t i=0; i<num_my_row; ++i) {
968 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
974 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
975 matrix->apply(*
x, *
y);
978 RCP<const Tpetra_Map> flat_x_map, flat_y_map;
979 RCP<const Tpetra_CrsGraph> flat_graph =
981 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
985 RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
986 RCP<Flat_Tpetra_Vector> flat_x =
988 RCP<Flat_Tpetra_Vector> flat_y =
990 flat_matrix->apply(*flat_x, *flat_y);
996 BaseScalar tol = 1.0e-14;
997 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
998 ArrayRCP<Scalar> y2_view = y2->get1dViewNonConst();
999 for (
size_t i=0; i<num_my_row; ++i) {
1000 TEST_EQUALITY( y_view[i].size(),
VectorSize );
1001 TEST_EQUALITY( y2_view[i].size(),
VectorSize );
1016 using Teuchos::ArrayView;
1017 using Teuchos::Array;
1018 using Teuchos::ArrayRCP;
1019 using Teuchos::ParameterList;
1024 typedef Teuchos::Comm<int> Tpetra_Comm;
1025 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1026 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1027 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1028 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1031 if ( !Kokkos::is_initialized() )
1032 Kokkos::initialize();
1036 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1037 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1038 RCP<const Tpetra_Map> map =
1039 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1041 RCP<Tpetra_CrsGraph> graph =
1042 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1043 Array<GlobalOrdinal> columnIndices(3);
1044 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1045 const size_t num_my_row = myGIDs.size();
1046 for (
size_t i=0; i<num_my_row; ++i) {
1048 if (row == 0 || row == nrow-1) {
1049 columnIndices[0] = row;
1050 graph->insertGlobalIndices(row, columnIndices(0,1));
1053 columnIndices[0] = row-1;
1054 columnIndices[1] = row;
1055 columnIndices[2] = row+1;
1056 graph->insertGlobalIndices(row, columnIndices(0,3));
1059 graph->fillComplete();
1060 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1063 Array<Scalar> vals(3);
1066 a_val.fastAccessCoeff(
j) =
1067 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1069 for (
size_t i=0; i<num_my_row; ++i) {
1071 if (row == 0 || row == nrow-1) {
1072 columnIndices[0] = row;
1074 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1077 columnIndices[0] = row-1;
1078 columnIndices[1] = row;
1079 columnIndices[2] = row+1;
1080 vals[0] =
Scalar(-1.0) * a_val;
1081 vals[1] =
Scalar(2.0) * a_val;
1082 vals[2] =
Scalar(-1.0) * a_val;
1083 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1086 matrix->fillComplete();
1088 matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
1089 Teuchos::VERB_EXTREME);
1092 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1093 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1096 b_val.fastAccessCoeff(
j) =
1097 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1099 for (
size_t i=0; i<num_my_row; ++i) {
1101 if (row == 0 || row == nrow-1)
1104 b_view[i] = -
Scalar(b_val * h * h);
1108 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1109 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1110 typedef typename BST::mag_type base_mag_type;
1111 typedef typename Tpetra_Vector::mag_type mag_type;
1112 base_mag_type btol = 1e-9;
1113 mag_type tol = btol;
1116 out.getOStream().get());
1117 TEST_EQUALITY_CONST( solved,
true );
1124 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1126 for (
size_t i=0; i<num_my_row; ++i) {
1128 BaseScalar xx = row * h;
1130 val.fastAccessCoeff(
j) =
1131 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1133 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1139 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1141 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1145 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), btol);
1150 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2) 1160 using Teuchos::ArrayView;
1161 using Teuchos::Array;
1162 using Teuchos::ArrayRCP;
1163 using Teuchos::ParameterList;
1164 using Teuchos::getParametersFromXmlFile;
1169 typedef Teuchos::Comm<int> Tpetra_Comm;
1170 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1171 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1172 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1173 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1176 if ( !Kokkos::is_initialized() )
1177 Kokkos::initialize();
1181 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1182 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1183 RCP<const Tpetra_Map> map =
1184 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1186 RCP<Tpetra_CrsGraph> graph =
1187 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1188 Array<GlobalOrdinal> columnIndices(3);
1189 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1190 const size_t num_my_row = myGIDs.size();
1191 for (
size_t i=0; i<num_my_row; ++i) {
1193 if (row == 0 || row == nrow-1) {
1194 columnIndices[0] = row;
1195 graph->insertGlobalIndices(row, columnIndices(0,1));
1198 columnIndices[0] = row-1;
1199 columnIndices[1] = row;
1200 columnIndices[2] = row+1;
1201 graph->insertGlobalIndices(row, columnIndices(0,3));
1204 graph->fillComplete();
1205 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1208 Array<Scalar> vals(3);
1211 a_val.fastAccessCoeff(
j) =
1212 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1214 for (
size_t i=0; i<num_my_row; ++i) {
1216 if (row == 0 || row == nrow-1) {
1217 columnIndices[0] = row;
1219 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1222 columnIndices[0] = row-1;
1223 columnIndices[1] = row;
1224 columnIndices[2] = row+1;
1225 vals[0] =
Scalar(-1.0) * a_val;
1226 vals[1] =
Scalar(2.0) * a_val;
1227 vals[2] =
Scalar(-1.0) * a_val;
1228 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1231 matrix->fillComplete();
1234 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1235 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1238 b_val.fastAccessCoeff(
j) =
1239 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1241 for (
size_t i=0; i<num_my_row; ++i) {
1243 if (row == 0 || row == nrow-1)
1246 b_view[i] = -
Scalar(b_val * h * h);
1250 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1251 RCP<OP> matrix_op = matrix;
1252 RCP<ParameterList> muelu_params =
1253 getParametersFromXmlFile(
"muelu_cheby.xml");
1255 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1258 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1259 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1260 typedef typename BST::mag_type base_mag_type;
1261 typedef typename Tpetra_Vector::mag_type mag_type;
1262 base_mag_type btol = 1e-9;
1263 mag_type tol = btol;
1266 out.getOStream().get());
1267 TEST_EQUALITY_CONST( solved,
true );
1274 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1276 for (
size_t i=0; i<num_my_row; ++i) {
1278 BaseScalar xx = row * h;
1280 val.fastAccessCoeff(
j) =
1281 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1283 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1288 if (BST::magnitude(v.coeff(
j)) < btol)
1289 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1290 if (BST::magnitude(
val.coeff(
j)) < btol)
1291 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1295 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), btol);
1307 #if defined(HAVE_STOKHOS_BELOS) 1317 using Teuchos::ArrayView;
1318 using Teuchos::Array;
1319 using Teuchos::ArrayRCP;
1320 using Teuchos::ParameterList;
1325 typedef Teuchos::Comm<int> Tpetra_Comm;
1326 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1327 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1328 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1329 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1332 if ( !Kokkos::is_initialized() )
1333 Kokkos::initialize();
1337 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1338 RCP<const Tpetra_Map> map =
1339 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1341 RCP<Tpetra_CrsGraph> graph =
1342 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1343 Array<GlobalOrdinal> columnIndices(2);
1344 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1345 const size_t num_my_row = myGIDs.size();
1346 for (
size_t i=0; i<num_my_row; ++i) {
1348 columnIndices[0] = row;
1350 if (row != nrow-1) {
1351 columnIndices[1] = row+1;
1354 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1356 graph->fillComplete();
1357 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1360 Array<Scalar> vals(2);
1362 for (
size_t i=0; i<num_my_row; ++i) {
1364 columnIndices[0] = row;
1366 val.fastAccessCoeff(
j) =
j+1;
1370 if (row != nrow-1) {
1371 columnIndices[1] = row+1;
1373 val.fastAccessCoeff(
j) =
j+1;
1377 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1379 matrix->fillComplete();
1382 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1383 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1384 for (
size_t i=0; i<num_my_row; ++i) {
1389 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1390 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1391 typedef BaseScalar BelosScalar;
1393 typedef Scalar BelosScalar;
1395 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1396 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1397 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1398 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1399 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1400 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1401 typename ST::magnitudeType tol = 1e-9;
1402 belosParams->set(
"Flexible Gmres",
false);
1403 belosParams->set(
"Num Blocks", 100);
1404 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1405 belosParams->set(
"Maximum Iterations", 100);
1406 belosParams->set(
"Verbosity", 33);
1407 belosParams->set(
"Output Style", 1);
1408 belosParams->set(
"Output Frequency", 1);
1409 belosParams->set(
"Output Stream", out.getOStream());
1410 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1411 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1412 problem->setProblem();
1413 Belos::ReturnType ret = solver->solve();
1414 TEST_EQUALITY_CONST( ret, Belos::Converged );
1426 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1427 for (
size_t i=0; i<num_my_row; ++i) {
1431 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
1436 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1441 if (ST::magnitude(v.coeff(
j)) < tol)
1442 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1446 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1458 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) 1469 using Teuchos::ArrayView;
1470 using Teuchos::Array;
1471 using Teuchos::ArrayRCP;
1472 using Teuchos::ParameterList;
1477 typedef Teuchos::Comm<int> Tpetra_Comm;
1478 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1479 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1480 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1481 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1484 if ( !Kokkos::is_initialized() )
1485 Kokkos::initialize();
1489 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1490 RCP<const Tpetra_Map> map =
1491 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1493 RCP<Tpetra_CrsGraph> graph =
1494 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1495 Array<GlobalOrdinal> columnIndices(2);
1496 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1497 const size_t num_my_row = myGIDs.size();
1498 for (
size_t i=0; i<num_my_row; ++i) {
1500 columnIndices[0] = row;
1502 if (row != nrow-1) {
1503 columnIndices[1] = row+1;
1506 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1508 graph->fillComplete();
1509 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1512 Array<Scalar> vals(2);
1514 for (
size_t i=0; i<num_my_row; ++i) {
1516 columnIndices[0] = row;
1518 val.fastAccessCoeff(
j) =
j+1;
1522 if (row != nrow-1) {
1523 columnIndices[1] = row+1;
1525 val.fastAccessCoeff(
j) =
j+1;
1529 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
1531 matrix->fillComplete();
1534 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1535 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1536 for (
size_t i=0; i<num_my_row; ++i) {
1541 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
1542 Ifpack2::Factory factory;
1543 RCP<Prec> M = factory.create<Tpetra_CrsMatrix>(
"RILUK", matrix);
1548 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1549 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1550 typedef BaseScalar BelosScalar;
1552 typedef Scalar BelosScalar;
1554 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1555 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1556 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1557 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1558 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1559 problem->setRightPrec(M);
1560 problem->setProblem();
1561 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1562 typename ST::magnitudeType tol = 1e-9;
1563 belosParams->set(
"Flexible Gmres",
false);
1564 belosParams->set(
"Num Blocks", 100);
1565 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1566 belosParams->set(
"Maximum Iterations", 100);
1567 belosParams->set(
"Verbosity", 33);
1568 belosParams->set(
"Output Style", 1);
1569 belosParams->set(
"Output Frequency", 1);
1570 belosParams->set(
"Output Stream", out.getOStream());
1572 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1573 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1574 Belos::ReturnType ret = solver->solve();
1575 TEST_EQUALITY_CONST( ret, Belos::Converged );
1587 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1588 for (
size_t i=0; i<num_my_row; ++i) {
1592 val.fastAccessCoeff(
j) = BaseScalar(1.0) / BaseScalar(
j+1);
1597 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1602 if (ST::magnitude(v.coeff(
j)) < tol)
1603 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1607 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1619 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) 1629 using Teuchos::ArrayView;
1630 using Teuchos::Array;
1631 using Teuchos::ArrayRCP;
1632 using Teuchos::ParameterList;
1633 using Teuchos::getParametersFromXmlFile;
1638 typedef Teuchos::Comm<int> Tpetra_Comm;
1639 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1640 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1641 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1642 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1645 if ( !Kokkos::is_initialized() )
1646 Kokkos::initialize();
1650 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1651 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1652 RCP<const Tpetra_Map> map =
1653 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1655 RCP<Tpetra_CrsGraph> graph =
1656 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1657 Array<GlobalOrdinal> columnIndices(3);
1658 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1659 const size_t num_my_row = myGIDs.size();
1660 for (
size_t i=0; i<num_my_row; ++i) {
1662 if (row == 0 || row == nrow-1) {
1663 columnIndices[0] = row;
1664 graph->insertGlobalIndices(row, columnIndices(0,1));
1667 columnIndices[0] = row-1;
1668 columnIndices[1] = row;
1669 columnIndices[2] = row+1;
1670 graph->insertGlobalIndices(row, columnIndices(0,3));
1673 graph->fillComplete();
1674 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1677 Array<Scalar> vals(3);
1680 a_val.fastAccessCoeff(
j) =
1681 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1683 for (
size_t i=0; i<num_my_row; ++i) {
1685 if (row == 0 || row == nrow-1) {
1686 columnIndices[0] = row;
1688 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1691 columnIndices[0] = row-1;
1692 columnIndices[1] = row;
1693 columnIndices[2] = row+1;
1694 vals[0] =
Scalar(-1.0) * a_val;
1695 vals[1] =
Scalar(2.0) * a_val;
1696 vals[2] =
Scalar(-1.0) * a_val;
1697 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1700 matrix->fillComplete();
1703 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1704 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1707 b_val.fastAccessCoeff(
j) =
1708 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1710 for (
size_t i=0; i<num_my_row; ++i) {
1712 if (row == 0 || row == nrow-1)
1715 b_view[i] = -
Scalar(b_val * h * h);
1719 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1720 RCP<ParameterList> muelu_params =
1721 getParametersFromXmlFile(
"muelu_cheby.xml");
1722 RCP<OP> matrix_op = matrix;
1724 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(matrix_op, *muelu_params);
1727 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1728 #ifdef HAVE_STOKHOS_ENSEMBLE_REDUCT 1729 typedef BaseScalar BelosScalar;
1731 typedef Scalar BelosScalar;
1733 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1734 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1735 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1736 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1737 problem->setRightPrec(M);
1738 problem->setProblem();
1739 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1740 typename ST::magnitudeType tol = 1e-9;
1741 belosParams->set(
"Num Blocks", 100);
1742 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1743 belosParams->set(
"Maximum Iterations", 100);
1744 belosParams->set(
"Verbosity", 33);
1745 belosParams->set(
"Output Style", 1);
1746 belosParams->set(
"Output Frequency", 1);
1747 belosParams->set(
"Output Stream", out.getOStream());
1750 belosParams->set(
"Implicit Residual Scaling",
"None");
1752 RCP<Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true> > solver =
1753 rcp(
new Belos::PseudoBlockCGSolMgr<BelosScalar,MV,OP,true>(problem, belosParams));
1754 Belos::ReturnType ret = solver->solve();
1755 TEST_EQUALITY_CONST( ret, Belos::Converged );
1757 #ifndef HAVE_STOKHOS_ENSEMBLE_REDUCT 1759 std::vector<int> ensemble_iterations =
1760 solver->getResidualStatusTest()->getEnsembleIterations();
1761 out <<
"Ensemble iterations = ";
1763 out << ensemble_iterations[i] <<
" ";
1772 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1774 for (
size_t i=0; i<num_my_row; ++i) {
1776 BaseScalar xx = row * h;
1778 val.fastAccessCoeff(
j) =
1779 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1781 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1786 if (ST::magnitude(v.coeff(
j)) < tol)
1787 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1788 if (ST::magnitude(
val.coeff(
j)) < tol)
1789 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1793 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1806 #if defined(HAVE_STOKHOS_AMESOS2) 1816 using Teuchos::ArrayView;
1817 using Teuchos::Array;
1818 using Teuchos::ArrayRCP;
1819 using Teuchos::ParameterList;
1824 typedef Teuchos::Comm<int> Tpetra_Comm;
1825 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1826 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1827 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
1828 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1829 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1832 if ( !Kokkos::is_initialized() )
1833 Kokkos::initialize();
1837 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1838 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1839 RCP<const Tpetra_Map> map =
1840 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1842 RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map,
size_t(3));
1843 Array<GlobalOrdinal> columnIndices(3);
1844 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1845 const size_t num_my_row = myGIDs.size();
1846 for (
size_t i=0; i<num_my_row; ++i) {
1848 if (row == 0 || row == nrow-1) {
1849 columnIndices[0] = row;
1850 graph->insertGlobalIndices(row, columnIndices(0,1));
1853 columnIndices[0] = row-1;
1854 columnIndices[1] = row;
1855 columnIndices[2] = row+1;
1856 graph->insertGlobalIndices(row, columnIndices(0,3));
1859 graph->fillComplete();
1860 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1863 Array<Scalar> vals(3);
1866 a_val.fastAccessCoeff(
j) =
1867 BaseScalar(1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1869 for (
size_t i=0; i<num_my_row; ++i) {
1871 if (row == 0 || row == nrow-1) {
1872 columnIndices[0] = row;
1874 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1877 columnIndices[0] = row-1;
1878 columnIndices[1] = row;
1879 columnIndices[2] = row+1;
1880 vals[0] =
Scalar(-1.0) * a_val;
1881 vals[1] =
Scalar(2.0) * a_val;
1882 vals[2] =
Scalar(-1.0) * a_val;
1883 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1886 matrix->fillComplete();
1889 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1890 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1893 b_val.fastAccessCoeff(
j) =
1894 BaseScalar(-1.0) + BaseScalar(
j) / BaseScalar(
VectorSize);
1896 for (
size_t i=0; i<num_my_row; ++i) {
1898 if (row == 0 || row == nrow-1)
1901 b_view[i] = -
Scalar(b_val * h * h);
1905 typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
1906 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1907 std::string solver_name;
1908 #if defined(HAVE_AMESOS2_BASKER) 1909 solver_name =
"basker";
1910 #elif defined(HAVE_AMESOS2_KLU2) 1911 solver_name =
"klu2";
1912 #elif defined(HAVE_AMESOS2_SUPERLUDIST) 1913 solver_name =
"superlu_dist";
1914 #elif defined(HAVE_AMESOS2_SUPERLUMT) 1915 solver_name =
"superlu_mt";
1916 #elif defined(HAVE_AMESOS2_SUPERLU) 1917 solver_name =
"superlu";
1918 #elif defined(HAVE_AMESOS2_PARDISO_MKL) 1919 solver_name =
"pardisomkl";
1920 #elif defined(HAVE_AMESOS2_LAPACK) 1921 solver_name =
"lapack";
1922 #elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL) 1923 solver_name =
"lapack";
1929 out <<
"Solving linear system with " << solver_name << std::endl;
1930 RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
1931 solver_name, matrix,
x, b);
1938 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1939 typename ST::magnitudeType tol = 1e-9;
1940 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1942 for (
size_t i=0; i<num_my_row; ++i) {
1944 BaseScalar xx = row * h;
1946 val.fastAccessCoeff(
j) =
1947 BaseScalar(0.5) * (b_val.coeff(
j)/a_val.coeff(
j)) * xx * (xx - BaseScalar(1.0));
1949 TEST_EQUALITY( x_view[i].size(),
VectorSize );
1954 if (ST::magnitude(v.coeff(
j)) < tol)
1955 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1956 if (ST::magnitude(
val.coeff(
j)) < tol)
1957 val.fastAccessCoeff(
j) = BaseScalar(0.0);
1961 TEST_FLOATING_EQUALITY(v.coeff(
j),
val.coeff(
j), tol);
1973 #define CRSMATRIX_MP_VECTOR_TESTS_SLGN(S, LO, GO, N) \ 1974 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorAdd, S, LO, GO, N ) \ 1975 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, VectorDot, S, LO, GO, N ) \ 1976 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorAdd, S, LO, GO, N ) \ 1977 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDot, S, LO, GO, N ) \ 1978 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MultiVectorDotSub, S, LO, GO, N ) \ 1979 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixVectorMultiply, S, LO, GO, N ) \ 1980 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, MatrixMultiVectorMultiply, S, LO, GO, N ) \ 1981 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Flatten, S, LO, GO, N ) \ 1982 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimpleCG, S, LO, GO, N ) \ 1983 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, SimplePCG_Muelu, S, LO, GO, N ) \ 1984 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES, S, LO, GO, N ) \ 1985 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosGMRES_RILUK, S, LO, GO, N ) \ 1986 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, BelosCG_Muelu, S, LO, GO, N ) \ 1987 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_MP, Amesos2, S, LO, GO, N ) 1989 #define CRSMATRIX_MP_VECTOR_TESTS_N_SFS(N) \ 1990 typedef Stokhos::DeviceForNode<N>::type Device; \ 1991 typedef Stokhos::StaticFixedStorage<int,double,VectorSize,Device::execution_space> SFS; \ 1992 CRSMATRIX_MP_VECTOR_TESTS_SLGN(SFS, int, int, N) 1994 #define CRSMATRIX_MP_VECTOR_TESTS_N(N) \ 1995 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)
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)