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