42 #include "Teuchos_UnitTestHelpers.hpp" 46 #include "Teuchos_XMLParameterListCoreHelpers.hpp" 52 #include "Tpetra_ConfigDefs.hpp" 53 #include "Tpetra_DefaultPlatform.hpp" 54 #include "Tpetra_Map.hpp" 55 #include "Tpetra_MultiVector.hpp" 56 #include "Tpetra_Vector.hpp" 57 #include "Tpetra_CrsGraph.hpp" 58 #include "Tpetra_CrsMatrix.hpp" 62 #ifdef HAVE_STOKHOS_BELOS 64 #include "BelosLinearProblem.hpp" 65 #include "BelosPseudoBlockGmresSolMgr.hpp" 66 #include "BelosPseudoBlockCGSolMgr.hpp" 70 #ifdef HAVE_STOKHOS_IFPACK2 72 #include "Ifpack2_Factory.hpp" 76 #ifdef HAVE_STOKHOS_MUELU 78 #include "MueLu_CreateTpetraPreconditioner.hpp" 82 #ifdef HAVE_STOKHOS_AMESOS2 84 #include "Amesos2_Factory.hpp" 92 template <
typename scalar,
typename ordinal>
96 const ordinal iColFEM,
97 const ordinal iStoch )
99 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
100 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
101 return X_fem + X_stoch;
105 template <
typename scalar,
typename ordinal>
109 const ordinal nStoch,
110 const ordinal iColFEM,
112 const ordinal iStoch)
114 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
115 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
116 const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
117 return X_fem + X_stoch + X_col;
121 template <
typename scalar,
typename ordinal>
124 const ordinal nStoch,
125 const ordinal iRowFEM,
126 const ordinal iColFEM,
127 const ordinal iStoch )
129 const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
130 ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
132 const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
134 return A_fem + A_stoch;
138 template <
typename kokkos_cijk_type,
typename ordinal_type>
144 using Teuchos::Array;
154 Array< RCP<const one_d_basis> > bases(
stoch_dim);
156 bases[i] = rcp(
new legendre_basis(
poly_ord,
true));
157 RCP<const product_basis> basis = rcp(
new product_basis(bases));
160 RCP<Cijk>
cijk = basis->computeTripleProductTensor();
163 kokkos_cijk_type kokkos_cijk =
164 Stokhos::create_product_tensor<execution_space>(*basis, *
cijk);
185 using Teuchos::ArrayView;
186 using Teuchos::Array;
187 using Teuchos::ArrayRCP;
192 typedef typename Scalar::cijk_type Cijk;
194 typedef Teuchos::Comm<int> Tpetra_Comm;
195 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
199 if (!Kokkos::HostSpace::execution_space::is_initialized())
200 Kokkos::HostSpace::execution_space::initialize();
201 if (!execution_space::is_initialized())
202 execution_space::initialize();
210 RCP<const Tpetra_Comm> comm =
211 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
215 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
216 RCP<const Tpetra_Map> map =
217 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
219 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
220 const size_t num_my_row = myGIDs.size();
223 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
224 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
225 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
226 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
228 for (
size_t i=0; i<num_my_row; ++i) {
231 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
232 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
237 x1_view = Teuchos::null;
238 x2_view = Teuchos::null;
243 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
244 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
250 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
252 BaseScalar tol = 1.0e-14;
253 for (
size_t i=0; i<num_my_row; ++i) {
256 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
257 nrow, pce_size, row,
j);
258 val.fastAccessCoeff(
j) = alpha.coeff(0)*v + 0.12345*beta.coeff(0)*v;
260 TEST_EQUALITY( y_view[i].size(), pce_size );
262 TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(
j),
val.fastAccessCoeff(
j), tol );
277 using Teuchos::ArrayView;
278 using Teuchos::Array;
279 using Teuchos::ArrayRCP;
284 typedef typename Scalar::cijk_type Cijk;
286 typedef Teuchos::Comm<int> Tpetra_Comm;
287 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
289 typedef typename Tpetra_Vector::dot_type dot_type;
292 if (!Kokkos::HostSpace::execution_space::is_initialized())
293 Kokkos::HostSpace::execution_space::initialize();
294 if (!Device::is_initialized())
295 Device::initialize();
303 RCP<const Tpetra_Comm> comm =
304 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
308 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
309 RCP<const Tpetra_Map> map =
310 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
312 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
313 const size_t num_my_row = myGIDs.size();
316 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
317 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
318 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
319 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
321 for (
size_t i=0; i<num_my_row; ++i) {
324 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
325 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
330 x1_view = Teuchos::null;
331 x2_view = Teuchos::null;
334 dot_type
dot = x1->dot(*x2);
339 dot_type local_val(0);
340 for (
size_t i=0; i<num_my_row; ++i) {
343 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
344 nrow, pce_size, row,
j);
345 local_val += 0.12345 * v * v;
351 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
352 Teuchos::outArg(
val));
354 out <<
"dot = " <<
dot <<
" expected = " <<
val << std::endl;
356 BaseScalar tol = 1.0e-14;
357 TEST_FLOATING_EQUALITY(
dot,
val, tol );
371 using Teuchos::ArrayView;
372 using Teuchos::Array;
373 using Teuchos::ArrayRCP;
378 typedef typename Scalar::cijk_type Cijk;
380 typedef Teuchos::Comm<int> Tpetra_Comm;
381 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
385 if (!Kokkos::HostSpace::execution_space::is_initialized())
386 Kokkos::HostSpace::execution_space::initialize();
387 if (!Device::is_initialized())
388 Device::initialize();
396 RCP<const Tpetra_Comm> comm =
397 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
401 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
402 RCP<const Tpetra_Map> map =
403 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
405 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
406 const size_t num_my_row = myGIDs.size();
410 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
411 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
412 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
413 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
415 for (
size_t i=0; i<num_my_row; ++i) {
417 for (
size_t j=0;
j<ncol; ++
j) {
420 generate_multi_vector_coefficient<BaseScalar,size_t>(
421 nrow, ncol, pce_size, row,
j, k);
422 val1.fastAccessCoeff(k) = v;
423 val2.fastAccessCoeff(k) = 0.12345 * v;
425 x1_view[
j][i] = val1;
426 x2_view[
j][i] = val2;
429 x1_view = Teuchos::null;
430 x2_view = Teuchos::null;
435 RCP<Tpetra_MultiVector>
y = Tpetra::createMultiVector<Scalar>(map, ncol);
436 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
442 ArrayRCP< ArrayRCP<Scalar> > y_view =
y->get2dViewNonConst();
444 BaseScalar tol = 1.0e-14;
445 for (
size_t i=0; i<num_my_row; ++i) {
447 for (
size_t j=0;
j<ncol; ++
j) {
449 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
450 nrow, ncol, pce_size, row,
j, k);
451 val.fastAccessCoeff(k) = alpha.coeff(0)*v + 0.12345*beta.coeff(0)*v;
453 TEST_EQUALITY( y_view[
j][i].size(), pce_size );
455 TEST_FLOATING_EQUALITY( y_view[
j][i].fastAccessCoeff(k),
456 val.fastAccessCoeff(k), tol );
472 using Teuchos::ArrayView;
473 using Teuchos::Array;
474 using Teuchos::ArrayRCP;
479 typedef typename Scalar::cijk_type Cijk;
481 typedef Teuchos::Comm<int> Tpetra_Comm;
482 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
484 typedef typename Tpetra_MultiVector::dot_type dot_type;
487 if (!Kokkos::HostSpace::execution_space::is_initialized())
488 Kokkos::HostSpace::execution_space::initialize();
489 if (!Device::is_initialized())
490 Device::initialize();
498 RCP<const Tpetra_Comm> comm =
499 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
503 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
504 RCP<const Tpetra_Map> map =
505 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
507 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
508 const size_t num_my_row = myGIDs.size();
512 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
513 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
514 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
515 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
517 for (
size_t i=0; i<num_my_row; ++i) {
519 for (
size_t j=0;
j<ncol; ++
j) {
522 generate_multi_vector_coefficient<BaseScalar,size_t>(
523 nrow, ncol, pce_size, row,
j, k);
524 val1.fastAccessCoeff(k) = v;
525 val2.fastAccessCoeff(k) = 0.12345 * v;
527 x1_view[
j][i] = val1;
528 x2_view[
j][i] = val2;
531 x1_view = Teuchos::null;
532 x2_view = Teuchos::null;
535 Array<dot_type> dots(ncol);
536 x1->dot(*x2, dots());
541 Array<dot_type> local_vals(ncol, dot_type(0));
542 for (
size_t i=0; i<num_my_row; ++i) {
544 for (
size_t j=0;
j<ncol; ++
j) {
546 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
547 nrow, ncol, pce_size, row,
j, k);
548 local_vals[
j] += 0.12345 * v * v;
554 Array<dot_type> vals(ncol, dot_type(0));
555 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
556 local_vals.getRawPtr(), vals.getRawPtr());
558 BaseScalar tol = 1.0e-14;
559 for (
size_t j=0;
j<ncol; ++
j) {
560 out <<
"dots(" <<
j <<
") = " << dots[
j]
561 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
562 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
577 using Teuchos::ArrayView;
578 using Teuchos::Array;
579 using Teuchos::ArrayRCP;
584 typedef typename Scalar::cijk_type Cijk;
586 typedef Teuchos::Comm<int> Tpetra_Comm;
587 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
589 typedef typename Tpetra_MultiVector::dot_type dot_type;
592 if (!Kokkos::HostSpace::execution_space::is_initialized())
593 Kokkos::HostSpace::execution_space::initialize();
594 if (!Device::is_initialized())
595 Device::initialize();
603 RCP<const Tpetra_Comm> comm =
604 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
608 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
609 RCP<const Tpetra_Map> map =
610 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
612 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
613 const size_t num_my_row = myGIDs.size();
617 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
618 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
619 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
620 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
622 for (
size_t i=0; i<num_my_row; ++i) {
624 for (
size_t j=0;
j<ncol; ++
j) {
627 generate_multi_vector_coefficient<BaseScalar,size_t>(
628 nrow, ncol, pce_size, row,
j, k);
629 val1.fastAccessCoeff(k) = v;
630 val2.fastAccessCoeff(k) = 0.12345 * v;
632 x1_view[
j][i] = val1;
633 x2_view[
j][i] = val2;
636 x1_view = Teuchos::null;
637 x2_view = Teuchos::null;
641 Teuchos::Array<size_t> cols(ncol_sub);
642 cols[0] = 4; cols[1] = 2;
643 RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
644 RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
647 Array<dot_type> dots(ncol_sub);
648 x1_sub->dot(*x2_sub, dots());
653 Array<dot_type> local_vals(ncol_sub, dot_type(0));
654 for (
size_t i=0; i<num_my_row; ++i) {
656 for (
size_t j=0;
j<ncol_sub; ++
j) {
658 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
659 nrow, ncol, pce_size, row, cols[
j], k);
660 local_vals[
j] += 0.12345 * v * v;
666 Array<dot_type> vals(ncol_sub, dot_type(0));
667 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
668 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
671 BaseScalar tol = 1.0e-14;
672 for (
size_t j=0;
j<ncol_sub; ++
j) {
673 out <<
"dots(" <<
j <<
") = " << dots[
j]
674 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
675 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
690 using Teuchos::ArrayView;
691 using Teuchos::Array;
692 using Teuchos::ArrayRCP;
697 typedef typename Scalar::cijk_type Cijk;
699 typedef Teuchos::Comm<int> Tpetra_Comm;
700 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
702 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
703 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
706 if (!Kokkos::HostSpace::execution_space::is_initialized())
707 Kokkos::HostSpace::execution_space::initialize();
708 if (!Device::is_initialized())
709 Device::initialize();
718 RCP<const Tpetra_Comm> comm =
719 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
720 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
721 RCP<const Tpetra_Map> map =
722 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
724 RCP<Tpetra_CrsGraph> graph =
725 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
726 Array<GlobalOrdinal> columnIndices(2);
727 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
728 const size_t num_my_row = myGIDs.size();
729 for (
size_t i=0; i<num_my_row; ++i) {
731 columnIndices[0] = row;
734 columnIndices[1] = row+1;
737 graph->insertGlobalIndices(row, columnIndices(0,ncol));
739 graph->fillComplete();
740 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
743 Array<Scalar> vals(2);
745 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
747 const size_t num_col = row == nrow - 1 ? 1 : 2;
748 for (
size_t local_col=0; local_col<num_col; ++local_col) {
750 columnIndices[local_col] = col;
752 val.fastAccessCoeff(k) =
753 generate_matrix_coefficient<BaseScalar,size_t>(
754 nrow, pce_size, row, col, k);
755 vals[local_col] =
val;
757 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
759 matrix->fillComplete();
762 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
763 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
764 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
767 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
768 nrow, pce_size, row,
j);
769 x_view[local_row] =
val;
771 x_view = Teuchos::null;
780 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
781 matrix->apply(*
x, *
y);
787 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
788 BaseScalar tol = 1.0e-14;
789 typename Cijk::HostMirror host_cijk =
792 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
794 const size_t num_col = row == nrow - 1 ? 1 : 2;
796 for (
size_t local_col=0; local_col<num_col; ++local_col) {
800 const LocalOrdinal entry_beg = host_cijk.entry_begin(i);
803 for (
LocalOrdinal entry = entry_beg; entry < entry_end; ++entry) {
806 const BaseScalar a_j =
807 generate_matrix_coefficient<BaseScalar,size_t>(
808 nrow, pce_size, row, col,
j);
809 const BaseScalar a_k =
810 generate_matrix_coefficient<BaseScalar,size_t>(
811 nrow, pce_size, row, col, k);
812 const BaseScalar x_j =
813 generate_vector_coefficient<BaseScalar,size_t>(
814 nrow, pce_size, col,
j);
815 const BaseScalar x_k =
816 generate_vector_coefficient<BaseScalar,size_t>(
817 nrow, pce_size, col, k);
818 tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
820 val.fastAccessCoeff(i) += tmp;
823 TEST_EQUALITY( y_view[local_row].size(), pce_size );
825 TEST_FLOATING_EQUALITY( y_view[local_row].fastAccessCoeff(i),
826 val.fastAccessCoeff(i), tol );
841 using Teuchos::ArrayView;
842 using Teuchos::Array;
843 using Teuchos::ArrayRCP;
848 typedef typename Scalar::cijk_type Cijk;
850 typedef Teuchos::Comm<int> Tpetra_Comm;
851 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
853 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
854 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
857 if (!Kokkos::HostSpace::execution_space::is_initialized())
858 Kokkos::HostSpace::execution_space::initialize();
859 if (!Device::is_initialized())
860 Device::initialize();
869 RCP<const Tpetra_Comm> comm =
870 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
871 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
872 RCP<const Tpetra_Map> map =
873 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
875 RCP<Tpetra_CrsGraph> graph =
876 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
877 Array<GlobalOrdinal> columnIndices(2);
878 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
879 const size_t num_my_row = myGIDs.size();
880 for (
size_t i=0; i<num_my_row; ++i) {
882 columnIndices[0] = row;
885 columnIndices[1] = row+1;
888 graph->insertGlobalIndices(row, columnIndices(0,ncol));
890 graph->fillComplete();
891 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
894 Array<Scalar> vals(2);
896 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
898 const size_t num_col = row == nrow - 1 ? 1 : 2;
899 for (
size_t local_col=0; local_col<num_col; ++local_col) {
901 columnIndices[local_col] = col;
903 val.fastAccessCoeff(k) =
904 generate_matrix_coefficient<BaseScalar,size_t>(
905 nrow, pce_size, row, col, k);
906 vals[local_col] =
val;
908 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
910 matrix->fillComplete();
914 RCP<Tpetra_MultiVector>
x = Tpetra::createMultiVector<Scalar>(map, ncol);
915 ArrayRCP< ArrayRCP<Scalar> > x_view =
x->get2dViewNonConst();
916 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
918 for (
size_t col=0; col<ncol; ++col) {
921 generate_multi_vector_coefficient<BaseScalar,size_t>(
922 nrow, ncol, pce_size, row, col, k);
923 val.fastAccessCoeff(k) = v;
925 x_view[col][local_row] =
val;
928 x_view = Teuchos::null;
937 RCP<Tpetra_MultiVector>
y = Tpetra::createMultiVector<Scalar>(map, ncol);
938 matrix->apply(*
x, *
y);
944 ArrayRCP< ArrayRCP<Scalar> > y_view =
y->get2dViewNonConst();
945 BaseScalar tol = 1.0e-14;
946 typename Cijk::HostMirror host_cijk =
949 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
951 for (
size_t xcol=0; xcol<ncol; ++xcol) {
952 const size_t num_col = row == nrow - 1 ? 1 : 2;
954 for (
size_t local_col=0; local_col<num_col; ++local_col) {
958 const LocalOrdinal entry_beg = host_cijk.entry_begin(i);
961 for (
LocalOrdinal entry = entry_beg; entry < entry_end; ++entry) {
964 const BaseScalar a_j =
965 generate_matrix_coefficient<BaseScalar,size_t>(
966 nrow, pce_size, row, col,
j);
967 const BaseScalar a_k =
968 generate_matrix_coefficient<BaseScalar,size_t>(
969 nrow, pce_size, row, col, k);
970 const BaseScalar x_j =
971 generate_multi_vector_coefficient<BaseScalar,size_t>(
972 nrow, ncol, pce_size, col, xcol,
j);
973 const BaseScalar x_k =
974 generate_multi_vector_coefficient<BaseScalar,size_t>(
975 nrow, ncol, pce_size, col, xcol, k);
976 tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
978 val.fastAccessCoeff(i) += tmp;
981 TEST_EQUALITY( y_view[xcol][local_row].size(), pce_size );
983 TEST_FLOATING_EQUALITY( y_view[xcol][local_row].fastAccessCoeff(i),
984 val.fastAccessCoeff(i), tol );
1000 using Teuchos::ArrayView;
1001 using Teuchos::Array;
1002 using Teuchos::ArrayRCP;
1007 typedef typename Scalar::cijk_type Cijk;
1009 typedef Teuchos::Comm<int> Tpetra_Comm;
1010 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1012 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1013 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1016 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1019 if (!Kokkos::HostSpace::execution_space::is_initialized())
1020 Kokkos::HostSpace::execution_space::initialize();
1021 if (!Device::is_initialized())
1022 Device::initialize();
1031 RCP<const Tpetra_Comm> comm =
1032 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1033 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1034 RCP<const Tpetra_Map> map =
1035 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1037 RCP<Tpetra_CrsGraph> graph =
1038 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1039 Array<GlobalOrdinal> columnIndices(2);
1040 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1041 const size_t num_my_row = myGIDs.size();
1042 for (
size_t i=0; i<num_my_row; ++i) {
1044 columnIndices[0] = row;
1046 if (row != nrow-1) {
1047 columnIndices[1] = row+1;
1050 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1052 graph->fillComplete();
1053 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1056 Array<Scalar> vals(2);
1058 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1060 const size_t num_col = row == nrow - 1 ? 1 : 2;
1061 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1063 columnIndices[local_col] = col;
1065 val.fastAccessCoeff(k) =
1066 generate_matrix_coefficient<BaseScalar,size_t>(
1067 nrow, pce_size, row, col, k);
1068 vals[local_col] =
val;
1070 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1072 matrix->fillComplete();
1075 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1076 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1077 for (
size_t i=0; i<num_my_row; ++i) {
1080 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
1081 nrow, pce_size, row,
j);
1086 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
1087 matrix->apply(*
x, *
y);
1104 RCP<const Tpetra_Map> flat_x_map, flat_y_map;
1105 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1107 Stokhos::create_flat_pce_graph(*graph,
cijk, flat_x_map, flat_y_map,
1108 cijk_graph, pce_size);
1109 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1113 RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
1114 RCP<Flat_Tpetra_Vector> flat_x =
1116 RCP<Flat_Tpetra_Vector> flat_y =
1118 flat_matrix->apply(*flat_x, *flat_y);
1138 BaseScalar tol = 1.0e-14;
1139 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
1140 ArrayRCP<Scalar> y2_view = y2->get1dViewNonConst();
1141 for (
size_t i=0; i<num_my_row; ++i) {
1142 TEST_EQUALITY( y_view[i].size(), pce_size );
1143 TEST_EQUALITY( y2_view[i].size(), pce_size );
1145 TEST_FLOATING_EQUALITY( y_view[i].fastAccessCoeff(
j),
1146 y2_view[i].fastAccessCoeff(
j), tol );
1161 using Teuchos::ArrayView;
1162 using Teuchos::Array;
1163 using Teuchos::ArrayRCP;
1164 using Teuchos::ParameterList;
1169 typedef typename Scalar::cijk_type Cijk;
1171 typedef Teuchos::Comm<int> Tpetra_Comm;
1172 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1174 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1175 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1178 if (!Kokkos::HostSpace::execution_space::is_initialized())
1179 Kokkos::HostSpace::execution_space::initialize();
1180 if (!Device::is_initialized())
1181 Device::initialize();
1190 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1191 RCP<const Tpetra_Comm> comm =
1192 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1193 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1194 RCP<const Tpetra_Map> map =
1195 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1197 RCP<Tpetra_CrsGraph> graph =
1198 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1199 Array<GlobalOrdinal> columnIndices(3);
1200 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1201 const size_t num_my_row = myGIDs.size();
1202 for (
size_t i=0; i<num_my_row; ++i) {
1204 if (row == 0 || row == nrow-1) {
1205 columnIndices[0] = row;
1206 graph->insertGlobalIndices(row, columnIndices(0,1));
1209 columnIndices[0] = row-1;
1210 columnIndices[1] = row;
1211 columnIndices[2] = row+1;
1212 graph->insertGlobalIndices(row, columnIndices(0,3));
1215 graph->fillComplete();
1216 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1219 Array<Scalar> vals(3);
1222 a_val.fastAccessCoeff(
j) =
1223 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1225 for (
size_t i=0; i<num_my_row; ++i) {
1227 if (row == 0 || row == nrow-1) {
1228 columnIndices[0] = row;
1230 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1233 columnIndices[0] = row-1;
1234 columnIndices[1] = row;
1235 columnIndices[2] = row+1;
1236 vals[0] =
Scalar(1.0) * a_val;
1237 vals[1] =
Scalar(-2.0) * a_val;
1238 vals[2] =
Scalar(1.0) * a_val;
1239 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1242 matrix->fillComplete();
1248 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1249 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1252 b_val.fastAccessCoeff(
j) =
1253 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1255 for (
size_t i=0; i<num_my_row; ++i) {
1257 if (row == 0 || row == nrow-1)
1260 b_view[i] = b_val * (h*h);
1267 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1268 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1269 typedef typename BST::mag_type base_mag_type;
1270 typedef typename Tpetra_Vector::mag_type mag_type;
1271 base_mag_type btol = 1e-9;
1272 mag_type tol = btol;
1275 out.getOStream().get());
1276 TEST_EQUALITY_CONST( solved,
true );
1283 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1284 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1285 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1287 Stokhos::create_flat_pce_graph(*graph,
cijk, flat_x_map, flat_b_map,
1288 cijk_graph, pce_size);
1289 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1291 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1292 RCP<Flat_Tpetra_Vector> flat_x =
1294 RCP<Flat_Tpetra_Vector> flat_b =
1297 tol, max_its, out.getOStream().get());
1298 TEST_EQUALITY_CONST( solved_flat,
true );
1301 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1302 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1303 for (
size_t i=0; i<num_my_row; ++i) {
1304 TEST_EQUALITY( x_view[i].size(), pce_size );
1305 TEST_EQUALITY( x2_view[i].size(), pce_size );
1311 if (
j < v.size() &&
BST::abs(v.coeff(
j)) < btol)
1312 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1313 if (
j < v2.size() &&
BST::abs(v2.coeff(
j)) < btol)
1314 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1318 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
1326 #if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) 1339 using Teuchos::ArrayView;
1340 using Teuchos::Array;
1341 using Teuchos::ArrayRCP;
1342 using Teuchos::ParameterList;
1343 using Teuchos::getParametersFromXmlFile;
1348 typedef typename Scalar::cijk_type Cijk;
1350 typedef Teuchos::Comm<int> Tpetra_Comm;
1351 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1353 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1354 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1357 if (!Kokkos::HostSpace::execution_space::is_initialized())
1358 Kokkos::HostSpace::execution_space::initialize();
1359 if (!Device::is_initialized())
1360 Device::initialize();
1369 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1370 RCP<const Tpetra_Comm> comm =
1371 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1372 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1373 RCP<const Tpetra_Map> map =
1374 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1376 RCP<Tpetra_CrsGraph> graph =
1377 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1378 Array<GlobalOrdinal> columnIndices(3);
1379 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1380 const size_t num_my_row = myGIDs.size();
1381 for (
size_t i=0; i<num_my_row; ++i) {
1383 if (row == 0 || row == nrow-1) {
1384 columnIndices[0] = row;
1385 graph->insertGlobalIndices(row, columnIndices(0,1));
1388 columnIndices[0] = row-1;
1389 columnIndices[1] = row;
1390 columnIndices[2] = row+1;
1391 graph->insertGlobalIndices(row, columnIndices(0,3));
1394 graph->fillComplete();
1395 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1398 Array<Scalar> vals(3);
1401 a_val.fastAccessCoeff(
j) =
1402 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1404 for (
size_t i=0; i<num_my_row; ++i) {
1406 if (row == 0 || row == nrow-1) {
1407 columnIndices[0] = row;
1409 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1412 columnIndices[0] = row-1;
1413 columnIndices[1] = row;
1414 columnIndices[2] = row+1;
1415 vals[0] =
Scalar(1.0) * a_val;
1416 vals[1] =
Scalar(-2.0) * a_val;
1417 vals[2] =
Scalar(1.0) * a_val;
1418 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1421 matrix->fillComplete();
1424 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1425 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1428 b_val.fastAccessCoeff(
j) =
1429 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1431 for (
size_t i=0; i<num_my_row; ++i) {
1433 if (row == 0 || row == nrow-1)
1436 b_view[i] = b_val * (h*h);
1440 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1442 Stokhos::create_mean_based_product_tensor<Device,typename Storage::ordinal_type,BaseScalar>();
1444 RCP<ParameterList> muelu_params =
1445 getParametersFromXmlFile(
"muelu_cheby.xml");
1447 RCP<OP> mean_matrix_op = mean_matrix;
1449 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1453 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1454 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1455 typedef typename BST::mag_type base_mag_type;
1456 typedef typename Tpetra_Vector::mag_type mag_type;
1457 base_mag_type btol = 1e-9;
1458 mag_type tol = btol;
1461 out.getOStream().get());
1462 TEST_EQUALITY_CONST( solved,
true );
1469 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1470 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1471 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1473 Stokhos::create_flat_pce_graph(*graph,
cijk, flat_x_map, flat_b_map,
1474 cijk_graph, pce_size);
1475 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1477 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1478 RCP<Flat_Tpetra_Vector> flat_x =
1480 RCP<Flat_Tpetra_Vector> flat_b =
1488 tol, max_its, out.getOStream().get());
1489 TEST_EQUALITY_CONST( solved_flat,
true );
1492 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1493 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1494 for (
size_t i=0; i<num_my_row; ++i) {
1495 TEST_EQUALITY( x_view[i].size(), pce_size );
1496 TEST_EQUALITY( x2_view[i].size(), pce_size );
1502 if (
j < v.size() &&
BST::abs(v.coeff(
j)) < btol)
1503 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1504 if (
j < v2.size() &&
BST::abs(v2.coeff(
j)) < btol)
1505 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1509 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
1524 #if defined(HAVE_STOKHOS_BELOS) 1534 using Teuchos::ArrayView;
1535 using Teuchos::Array;
1536 using Teuchos::ArrayRCP;
1537 using Teuchos::ParameterList;
1542 typedef typename Scalar::cijk_type Cijk;
1544 typedef Teuchos::Comm<int> Tpetra_Comm;
1545 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1547 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1548 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1551 if (!Kokkos::HostSpace::execution_space::is_initialized())
1552 Kokkos::HostSpace::execution_space::initialize();
1553 if (!Device::is_initialized())
1554 Device::initialize();
1563 RCP<const Tpetra_Comm> comm =
1564 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1565 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1566 RCP<const Tpetra_Map> map =
1567 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1569 RCP<Tpetra_CrsGraph> graph =
1570 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1571 Array<GlobalOrdinal> columnIndices(2);
1572 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1573 const size_t num_my_row = myGIDs.size();
1574 for (
size_t i=0; i<num_my_row; ++i) {
1576 columnIndices[0] = row;
1578 if (row != nrow-1) {
1579 columnIndices[1] = row+1;
1582 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1584 graph->fillComplete();
1585 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1588 Array<Scalar> vals(2);
1590 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1592 const size_t num_col = row == nrow - 1 ? 1 : 2;
1593 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1595 columnIndices[local_col] = col;
1597 val.fastAccessCoeff(k) =
1598 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1599 vals[local_col] =
val;
1601 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1603 matrix->fillComplete();
1606 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1607 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1608 for (
size_t i=0; i<num_my_row; ++i) {
1613 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1614 typedef BaseScalar BelosScalar;
1616 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1617 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1618 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1619 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1620 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1621 typename ST::magnitudeType tol = 1e-9;
1622 belosParams->set(
"Flexible Gmres",
false);
1623 belosParams->set(
"Num Blocks", 100);
1624 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1625 belosParams->set(
"Maximum Iterations", 100);
1626 belosParams->set(
"Verbosity", 33);
1627 belosParams->set(
"Output Style", 1);
1628 belosParams->set(
"Output Frequency", 1);
1629 belosParams->set(
"Output Stream", out.getOStream());
1630 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1631 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1632 problem->setProblem();
1633 Belos::ReturnType ret = solver->solve();
1634 TEST_EQUALITY_CONST( ret, Belos::Converged );
1641 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1642 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1643 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1645 Stokhos::create_flat_pce_graph(*graph,
cijk, flat_x_map, flat_b_map,
1646 cijk_graph, pce_size);
1647 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1649 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1650 RCP<Flat_Tpetra_Vector> flat_x =
1652 RCP<Flat_Tpetra_Vector> flat_b =
1655 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1656 typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1657 RCP< FBLinProb > flat_problem =
1658 rcp(
new FBLinProb(flat_matrix, flat_x, flat_b));
1659 RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1660 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1662 flat_problem->setProblem();
1663 Belos::ReturnType flat_ret = flat_solver->solve();
1664 TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
1666 typename ST::magnitudeType btol = 100*tol;
1667 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1668 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1669 for (
size_t i=0; i<num_my_row; ++i) {
1670 TEST_EQUALITY( x_view[i].size(), pce_size );
1671 TEST_EQUALITY( x2_view[i].size(), pce_size );
1677 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
1678 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1679 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
1680 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1684 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
1701 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) 1712 using Teuchos::ArrayView;
1713 using Teuchos::Array;
1714 using Teuchos::ArrayRCP;
1715 using Teuchos::ParameterList;
1720 typedef typename Scalar::cijk_type Cijk;
1722 typedef Teuchos::Comm<int> Tpetra_Comm;
1723 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1725 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1726 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1729 if (!Kokkos::HostSpace::execution_space::is_initialized())
1730 Kokkos::HostSpace::execution_space::initialize();
1731 if (!Device::is_initialized())
1732 Device::initialize();
1741 RCP<const Tpetra_Comm> comm =
1742 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1743 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1744 RCP<const Tpetra_Map> map =
1745 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1747 RCP<Tpetra_CrsGraph> graph =
1748 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
1749 Array<GlobalOrdinal> columnIndices(2);
1750 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1751 const size_t num_my_row = myGIDs.size();
1752 for (
size_t i=0; i<num_my_row; ++i) {
1754 columnIndices[0] = row;
1756 if (row != nrow-1) {
1757 columnIndices[1] = row+1;
1760 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1762 graph->fillComplete();
1763 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1766 Array<Scalar> vals(2);
1768 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1770 const size_t num_col = row == nrow - 1 ? 1 : 2;
1771 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1773 columnIndices[local_col] = col;
1775 val.fastAccessCoeff(k) =
1776 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1777 vals[local_col] =
val;
1779 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1781 matrix->fillComplete();
1787 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1788 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1789 for (
size_t i=0; i<num_my_row; ++i) {
1794 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
1795 Ifpack2::Factory factory;
1796 RCP<Prec> M = factory.create<Tpetra_CrsMatrix>(
"RILUK", mean_matrix);
1801 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1802 typedef BaseScalar BelosScalar;
1804 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1805 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1806 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1807 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
1808 problem->setRightPrec(M);
1809 problem->setProblem();
1810 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1811 typename ST::magnitudeType tol = 1e-9;
1812 belosParams->set(
"Flexible Gmres",
false);
1813 belosParams->set(
"Num Blocks", 100);
1814 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1815 belosParams->set(
"Maximum Iterations", 100);
1816 belosParams->set(
"Verbosity", 33);
1817 belosParams->set(
"Output Style", 1);
1818 belosParams->set(
"Output Frequency", 1);
1819 belosParams->set(
"Output Stream", out.getOStream());
1821 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1822 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1823 Belos::ReturnType ret = solver->solve();
1824 TEST_EQUALITY_CONST( ret, Belos::Converged );
1828 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1829 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1830 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1832 Stokhos::create_flat_pce_graph(*graph,
cijk, flat_x_map, flat_b_map,
1833 cijk_graph, pce_size);
1834 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1836 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1837 RCP<Flat_Tpetra_Vector> flat_x =
1839 RCP<Flat_Tpetra_Vector> flat_b =
1842 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1843 typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1844 RCP< FBLinProb > flat_problem =
1845 rcp(
new FBLinProb(flat_matrix, flat_x, flat_b));
1846 RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1847 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1849 flat_problem->setProblem();
1850 Belos::ReturnType flat_ret = flat_solver->solve();
1851 TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
1853 typename ST::magnitudeType btol = 100*tol;
1854 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1855 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1856 for (
size_t i=0; i<num_my_row; ++i) {
1857 TEST_EQUALITY( x_view[i].size(), pce_size );
1858 TEST_EQUALITY( x2_view[i].size(), pce_size );
1864 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
1865 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1866 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
1867 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1871 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
1886 #if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION) 1899 using Teuchos::ArrayView;
1900 using Teuchos::Array;
1901 using Teuchos::ArrayRCP;
1902 using Teuchos::ParameterList;
1903 using Teuchos::getParametersFromXmlFile;
1908 typedef typename Scalar::cijk_type Cijk;
1910 typedef Teuchos::Comm<int> Tpetra_Comm;
1911 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1913 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1914 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1917 if (!Kokkos::HostSpace::execution_space::is_initialized())
1918 Kokkos::HostSpace::execution_space::initialize();
1919 if (!Device::is_initialized())
1920 Device::initialize();
1929 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1930 RCP<const Tpetra_Comm> comm =
1931 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
1932 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
1933 RCP<const Tpetra_Map> map =
1934 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
1936 RCP<Tpetra_CrsGraph> graph =
1937 rcp(
new Tpetra_CrsGraph(map,
size_t(3), Tpetra::StaticProfile));
1938 Array<GlobalOrdinal> columnIndices(3);
1939 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
1940 const size_t num_my_row = myGIDs.size();
1941 for (
size_t i=0; i<num_my_row; ++i) {
1943 if (row == 0 || row == nrow-1) {
1944 columnIndices[0] = row;
1945 graph->insertGlobalIndices(row, columnIndices(0,1));
1948 columnIndices[0] = row-1;
1949 columnIndices[1] = row;
1950 columnIndices[2] = row+1;
1951 graph->insertGlobalIndices(row, columnIndices(0,3));
1954 graph->fillComplete();
1955 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1958 Array<Scalar> vals(3);
1961 a_val.fastAccessCoeff(
j) =
1962 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1964 for (
size_t i=0; i<num_my_row; ++i) {
1966 if (row == 0 || row == nrow-1) {
1967 columnIndices[0] = row;
1969 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1972 columnIndices[0] = row-1;
1973 columnIndices[1] = row;
1974 columnIndices[2] = row+1;
1975 vals[0] =
Scalar(1.0) * a_val;
1976 vals[1] =
Scalar(-2.0) * a_val;
1977 vals[2] =
Scalar(1.0) * a_val;
1978 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1981 matrix->fillComplete();
1984 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1985 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1988 b_val.fastAccessCoeff(
j) =
1989 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1991 for (
size_t i=0; i<num_my_row; ++i) {
1993 if (row == 0 || row == nrow-1)
1996 b_view[i] = b_val * (h*h);
2000 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
2002 Stokhos::create_mean_based_product_tensor<Device,typename Storage::ordinal_type,BaseScalar>();
2004 RCP<ParameterList> muelu_params =
2005 getParametersFromXmlFile(
"muelu_cheby.xml");
2007 RCP<OP> mean_matrix_op = mean_matrix;
2009 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
2013 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2014 typedef BaseScalar BelosScalar;
2016 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2017 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2018 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
2019 problem->setRightPrec(M);
2020 problem->setProblem();
2021 RCP<ParameterList> belosParams = rcp(
new ParameterList);
2022 typename ST::magnitudeType tol = 1e-9;
2023 belosParams->set(
"Num Blocks", 100);
2024 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
2025 belosParams->set(
"Maximum Iterations", 100);
2026 belosParams->set(
"Verbosity", 33);
2027 belosParams->set(
"Output Style", 1);
2028 belosParams->set(
"Output Frequency", 1);
2029 belosParams->set(
"Output Stream", out.getOStream());
2031 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2032 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2035 Belos::ReturnType ret = solver->solve();
2036 TEST_EQUALITY_CONST( ret, Belos::Converged );
2043 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2044 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2045 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2047 Stokhos::create_flat_pce_graph(*graph,
cijk, flat_x_map, flat_b_map,
2048 cijk_graph, pce_size);
2049 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2051 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2052 RCP<Flat_Tpetra_Vector> flat_x =
2054 RCP<Flat_Tpetra_Vector> flat_b =
2057 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
2058 typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
2059 RCP< FBLinProb > flat_problem =
2060 rcp(
new FBLinProb(flat_matrix, flat_x, flat_b));
2061 RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
2062 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
2067 flat_problem->setProblem();
2068 Belos::ReturnType flat_ret = flat_solver->solve();
2069 TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
2071 typename ST::magnitudeType btol = 100*tol;
2072 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
2073 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2074 for (
size_t i=0; i<num_my_row; ++i) {
2075 TEST_EQUALITY( x_view[i].size(), pce_size );
2076 TEST_EQUALITY( x2_view[i].size(), pce_size );
2082 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
2083 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2084 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
2085 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
2089 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
2105 #if defined(HAVE_STOKHOS_AMESOS2) 2115 using Teuchos::ArrayView;
2116 using Teuchos::Array;
2117 using Teuchos::ArrayRCP;
2118 using Teuchos::ParameterList;
2123 typedef typename Scalar::cijk_type Cijk;
2125 typedef Teuchos::Comm<int> Tpetra_Comm;
2126 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2129 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2130 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2133 if (!Kokkos::HostSpace::execution_space::is_initialized())
2134 Kokkos::HostSpace::execution_space::initialize();
2135 if (!Device::is_initialized())
2136 Device::initialize();
2145 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
2146 RCP<const Tpetra_Comm> comm =
2147 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
2148 RCP<Node> node = KokkosClassic::Details::getNode<Node>();
2149 RCP<const Tpetra_Map> map =
2150 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal>(
2152 RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map,
size_t(3));
2153 Array<GlobalOrdinal> columnIndices(3);
2154 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
2155 const size_t num_my_row = myGIDs.size();
2156 for (
size_t i=0; i<num_my_row; ++i) {
2158 if (row == 0 || row == nrow-1) {
2159 columnIndices[0] = row;
2160 graph->insertGlobalIndices(row, columnIndices(0,1));
2163 columnIndices[0] = row-1;
2164 columnIndices[1] = row;
2165 columnIndices[2] = row+1;
2166 graph->insertGlobalIndices(row, columnIndices(0,3));
2169 graph->fillComplete();
2170 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
2173 Array<Scalar> vals(3);
2176 a_val.fastAccessCoeff(
j) =
2177 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
2179 for (
size_t i=0; i<num_my_row; ++i) {
2181 if (row == 0 || row == nrow-1) {
2182 columnIndices[0] = row;
2184 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2187 columnIndices[0] = row-1;
2188 columnIndices[1] = row;
2189 columnIndices[2] = row+1;
2190 vals[0] =
Scalar(1.0) * a_val;
2191 vals[1] =
Scalar(-2.0) * a_val;
2192 vals[2] =
Scalar(1.0) * a_val;
2193 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2196 matrix->fillComplete();
2199 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2200 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2203 b_val.fastAccessCoeff(
j) =
2204 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
2206 for (
size_t i=0; i<num_my_row; ++i) {
2208 if (row == 0 || row == nrow-1)
2211 b_view[i] = b_val * (h*h);
2215 typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
2216 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2217 std::string solver_name;
2218 #if defined(HAVE_AMESOS2_BASKER) 2219 solver_name =
"basker";
2220 #elif defined(HAVE_AMESOS2_KLU2) 2221 solver_name =
"klu2";
2222 #elif defined(HAVE_AMESOS2_SUPERLUDIST) 2223 solver_name =
"superlu_dist";
2224 #elif defined(HAVE_AMESOS2_SUPERLUMT) 2225 solver_name =
"superlu_mt";
2226 #elif defined(HAVE_AMESOS2_SUPERLU) 2227 solver_name =
"superlu";
2228 #elif defined(HAVE_AMESOS2_PARDISO_MKL) 2229 solver_name =
"pardisomkl";
2230 #elif defined(HAVE_AMESOS2_LAPACK) 2231 solver_name =
"lapack";
2232 #elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL) 2233 solver_name =
"lapack";
2239 out <<
"Solving linear system with " << solver_name << std::endl;
2240 RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2241 solver_name, matrix,
x, b);
2250 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2251 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2252 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2254 Stokhos::create_flat_pce_graph(*graph,
cijk, flat_x_map, flat_b_map,
2255 cijk_graph, pce_size);
2256 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2258 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2259 RCP<Flat_Tpetra_Vector> flat_x =
2261 RCP<Flat_Tpetra_Vector> flat_b =
2263 typedef Amesos2::Solver<Flat_Tpetra_CrsMatrix,Flat_Tpetra_MultiVector> Flat_Solver;
2264 RCP<Flat_Solver> flat_solver =
2265 Amesos2::create<Flat_Tpetra_CrsMatrix,Flat_Tpetra_MultiVector>(
2266 solver_name, flat_matrix, flat_x, flat_b);
2267 flat_solver->solve();
2269 typedef Kokkos::Details::ArithTraits<BaseScalar> ST;
2270 typename ST::mag_type btol = 1e-12;
2271 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
2272 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2273 for (
size_t i=0; i<num_my_row; ++i) {
2274 TEST_EQUALITY( x_view[i].size(), pce_size );
2275 TEST_EQUALITY( x2_view[i].size(), pce_size );
2281 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
2282 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2283 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
2284 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
2288 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
2303 #define CRSMATRIX_UQ_PCE_TESTS_SLGN(S, LO, GO, N) \ 2304 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorAdd, S, LO, GO, N ) \ 2305 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorDot, S, LO, GO, N ) \ 2306 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorAdd, S, LO, GO, N ) \ 2307 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDot, S, LO, GO, N ) \ 2308 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDotSub, S, LO, GO, N ) \ 2309 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixVectorMultiply, S, LO, GO, N ) \ 2310 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixMultiVectorMultiply, S, LO, GO, N ) \ 2311 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Flatten, S, LO, GO, N ) \ 2312 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimpleCG, S, LO, GO, N ) \ 2313 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, S, LO, GO, N ) \ 2314 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES, S, LO, GO, N ) \ 2315 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, S, LO, GO, N ) \ 2316 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosCG_Muelu, S, LO, GO, N ) \ 2317 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Amesos2, S, LO, GO, N ) 2319 #define CRSMATRIX_UQ_PCE_TESTS_N(N) \ 2320 typedef Stokhos::DeviceForNode2<N>::type Device; \ 2321 typedef Stokhos::DynamicStorage<int,double,Device::execution_space> DS; \ 2322 CRSMATRIX_UQ_PCE_TESTS_SLGN(DS, int, int, N)
Teuchos::RCP< Tpetra::CrsMatrix< Scalar, LO, GO, N > > build_mean_matrix(const Tpetra::CrsMatrix< Scalar, LO, GO, N > &A)
bool CG_Solve(const Matrix &A, Vector &x, const Vector &b, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
Kokkos::DefaultExecutionSpace execution_space
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Node > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::MP::Vector< Storage >, LocalOrdinal, GlobalOrdinal, Node > &vec_const, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Node > > &flat_map)
bool PCG_Solve(const Matrix &A, Vector &x, const Vector &b, const Prec &M, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
KokkosClassic::DefaultNode::DefaultNodeType Node
scalar generate_multi_vector_coefficient(const ordinal nFEM, const ordinal nVec, const ordinal nStoch, const ordinal iColFEM, const ordinal iVec, const ordinal iStoch)
void deep_copy(const Stokhos::CrsMatrix< ValueType, DstDevice, Layout > &dst, const Stokhos::CrsMatrix< ValueType, SrcDevice, Layout > &src)
KOKKOS_INLINE_FUNCTION PCE< Storage > abs(const PCE< Storage > &a)
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)
void setGlobalCijkTensor(const cijk_type &cijk)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typename CijkType< view_type >::type >::type cijk(const view_type &view)
Legendre polynomial basis.
Abstract base class for 1-D orthogonal polynomials.
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)
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
kokkos_cijk_type build_cijk(ordinal_type stoch_dim, ordinal_type poly_ord)
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_PCE, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node)
Stokhos::CrsMatrix< ValueType, Device, Layout >::HostMirror create_mirror_view(const Stokhos::CrsMatrix< ValueType, Device, Layout > &A)