42#include "Teuchos_UnitTestHelpers.hpp"
46#include "Teuchos_XMLParameterListCoreHelpers.hpp"
52#include "Tpetra_Core.hpp"
53#include "Tpetra_Map.hpp"
54#include "Tpetra_MultiVector.hpp"
55#include "Tpetra_Vector.hpp"
56#include "Tpetra_CrsGraph.hpp"
57#include "Tpetra_CrsMatrix.hpp"
61#ifdef HAVE_STOKHOS_BELOS
63#include "BelosLinearProblem.hpp"
64#include "BelosPseudoBlockGmresSolMgr.hpp"
65#include "BelosPseudoBlockCGSolMgr.hpp"
69#ifdef HAVE_STOKHOS_IFPACK2
71#include "Ifpack2_Factory.hpp"
75#ifdef HAVE_STOKHOS_MUELU
77#include "MueLu_CreateTpetraPreconditioner.hpp"
81#ifdef HAVE_STOKHOS_AMESOS2
83#include "Amesos2_Factory.hpp"
96#define USE_SCALAR_MEAN_BASED_PREC 1
98template <
typename scalar,
typename ordinal>
101 const ordinal nStoch,
102 const ordinal iColFEM,
103 const ordinal iStoch )
105 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
106 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
107 return X_fem + X_stoch;
111template <
typename scalar,
typename ordinal>
115 const ordinal nStoch,
116 const ordinal iColFEM,
118 const ordinal iStoch)
120 const scalar X_fem = 100.0 + scalar(iColFEM) / scalar(nFEM);
121 const scalar X_stoch = 1.0 + scalar(iStoch) / scalar(nStoch);
122 const scalar X_col = 0.01 + scalar(iVec) / scalar(nVec);
123 return X_fem + X_stoch + X_col;
127template <
typename scalar,
typename ordinal>
130 const ordinal nStoch,
131 const ordinal iRowFEM,
132 const ordinal iColFEM,
133 const ordinal iStoch )
135 const scalar A_fem = ( 10.0 + scalar(iRowFEM) / scalar(nFEM) ) +
136 ( 5.0 + scalar(iColFEM) / scalar(nFEM) );
138 const scalar A_stoch = ( 1.0 + scalar(iStoch) / scalar(nStoch) );
140 return A_fem + A_stoch;
144template <
typename kokkos_cijk_type,
typename ordinal_type>
150 using Teuchos::Array;
152 typedef typename kokkos_cijk_type::value_type value_type;
160 Array< RCP<const one_d_basis> > bases(
stoch_dim);
162 bases[i] = rcp(
new legendre_basis(
poly_ord,
true));
163 RCP<const product_basis> basis = rcp(
new product_basis(bases));
166 RCP<Cijk> cijk = basis->computeTripleProductTensor();
169 kokkos_cijk_type kokkos_cijk =
170 Stokhos::create_product_tensor<execution_space>(*basis, *cijk);
191 using Teuchos::ArrayView;
192 using Teuchos::Array;
193 using Teuchos::ArrayRCP;
195 typedef typename Storage::value_type BaseScalar;
197 typedef typename Scalar::cijk_type Cijk;
199 typedef Teuchos::Comm<int> Tpetra_Comm;
200 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
201 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
204 if ( !Kokkos::is_initialized() )
205 Kokkos::initialize();
213 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
217 RCP<const Tpetra_Map> map =
218 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
220 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
221 const size_t num_my_row = myGIDs.size();
224 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
225 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
226 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
227 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
228 Scalar val1(cijk), val2(cijk);
229 for (
size_t i=0; i<num_my_row; ++i) {
232 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
233 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
238 x1_view = Teuchos::null;
239 x2_view = Teuchos::null;
244 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
245 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
251 ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
253 BaseScalar tol = 1.0e-14;
254 for (
size_t i=0; i<num_my_row; ++i) {
257 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
258 nrow, pce_size, row,
j);
259 val.fastAccessCoeff(
j) = alpha.coeff(0)*v + 0.12345*beta.coeff(0)*v;
261 TEST_EQUALITY( y_view[i].size(), pce_size );
278 using Teuchos::ArrayView;
279 using Teuchos::Array;
280 using Teuchos::ArrayRCP;
282 typedef typename Storage::value_type BaseScalar;
284 typedef typename Scalar::cijk_type Cijk;
286 typedef Teuchos::Comm<int> Tpetra_Comm;
287 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
288 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
289 typedef typename Tpetra_Vector::dot_type dot_type;
292 if ( !Kokkos::is_initialized() )
293 Kokkos::initialize();
297 setGlobalCijkTensor(cijk);
301 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
305 RCP<const Tpetra_Map> map =
306 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
308 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
309 const size_t num_my_row = myGIDs.size();
312 RCP<Tpetra_Vector> x1 = Tpetra::createVector<Scalar>(map);
313 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
314 ArrayRCP<Scalar> x1_view = x1->get1dViewNonConst();
315 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
316 Scalar val1(cijk), val2(cijk);
317 for (
size_t i=0; i<num_my_row; ++i) {
320 val1.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
321 val2.fastAccessCoeff(
j) = 0.12345 * generate_vector_coefficient<BaseScalar,size_t>(nrow, pce_size, row,
j);
326 x1_view = Teuchos::null;
327 x2_view = Teuchos::null;
330 dot_type dot = x1->dot(*x2);
335 dot_type local_val(0);
336 for (
size_t i=0; i<num_my_row; ++i) {
339 BaseScalar v = generate_vector_coefficient<BaseScalar,size_t>(
340 nrow, pce_size, row,
j);
341 local_val += 0.12345 * v * v;
347 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, local_val,
348 Teuchos::outArg(
val));
350 out <<
"dot = " << dot <<
" expected = " <<
val << std::endl;
352 BaseScalar tol = 1.0e-14;
353 TEST_FLOATING_EQUALITY( dot,
val, tol );
367 using Teuchos::ArrayView;
368 using Teuchos::Array;
369 using Teuchos::ArrayRCP;
371 typedef typename Storage::value_type BaseScalar;
373 typedef typename Scalar::cijk_type Cijk;
375 typedef Teuchos::Comm<int> Tpetra_Comm;
376 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
377 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
380 if ( !Kokkos::is_initialized() )
381 Kokkos::initialize();
385 setGlobalCijkTensor(cijk);
389 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
393 RCP<const Tpetra_Map> map =
394 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
396 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
397 const size_t num_my_row = myGIDs.size();
401 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
402 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
403 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
404 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
405 Scalar val1(cijk), val2(cijk);
406 for (
size_t i=0; i<num_my_row; ++i) {
408 for (
size_t j=0;
j<ncol; ++
j) {
411 generate_multi_vector_coefficient<BaseScalar,size_t>(
412 nrow, ncol, pce_size, row,
j, k);
413 val1.fastAccessCoeff(k) = v;
414 val2.fastAccessCoeff(k) = 0.12345 * v;
416 x1_view[
j][i] = val1;
417 x2_view[
j][i] = val2;
420 x1_view = Teuchos::null;
421 x2_view = Teuchos::null;
426 RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
427 y->update(alpha, *x1, beta, *x2,
Scalar(0.0));
433 ArrayRCP< ArrayRCP<Scalar> > y_view = y->get2dViewNonConst();
435 BaseScalar tol = 1.0e-14;
436 for (
size_t i=0; i<num_my_row; ++i) {
438 for (
size_t j=0;
j<ncol; ++
j) {
440 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
441 nrow, ncol, pce_size, row,
j, k);
442 val.fastAccessCoeff(k) = alpha.coeff(0)*v + 0.12345*beta.coeff(0)*v;
444 TEST_EQUALITY( y_view[
j][i].size(), pce_size );
447 val.fastAccessCoeff(k), tol );
463 using Teuchos::ArrayView;
464 using Teuchos::Array;
465 using Teuchos::ArrayRCP;
467 typedef typename Storage::value_type BaseScalar;
469 typedef typename Scalar::cijk_type Cijk;
471 typedef Teuchos::Comm<int> Tpetra_Comm;
472 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
473 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
474 typedef typename Tpetra_MultiVector::dot_type dot_type;
477 if ( !Kokkos::is_initialized() )
478 Kokkos::initialize();
482 setGlobalCijkTensor(cijk);
486 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
490 RCP<const Tpetra_Map> map =
491 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
493 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
494 const size_t num_my_row = myGIDs.size();
498 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
499 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
500 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
501 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
502 Scalar val1(cijk), val2(cijk);
503 for (
size_t i=0; i<num_my_row; ++i) {
505 for (
size_t j=0;
j<ncol; ++
j) {
508 generate_multi_vector_coefficient<BaseScalar,size_t>(
509 nrow, ncol, pce_size, row,
j, k);
510 val1.fastAccessCoeff(k) = v;
511 val2.fastAccessCoeff(k) = 0.12345 * v;
513 x1_view[
j][i] = val1;
514 x2_view[
j][i] = val2;
517 x1_view = Teuchos::null;
518 x2_view = Teuchos::null;
521 Array<dot_type> dots(ncol);
522 x1->dot(*x2, dots());
527 Array<dot_type> local_vals(ncol, dot_type(0));
528 for (
size_t i=0; i<num_my_row; ++i) {
530 for (
size_t j=0;
j<ncol; ++
j) {
532 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
533 nrow, ncol, pce_size, row,
j, k);
534 local_vals[
j] += 0.12345 * v * v;
540 Array<dot_type> vals(ncol, dot_type(0));
541 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM, Teuchos::as<int>(ncol),
542 local_vals.getRawPtr(), vals.getRawPtr());
544 BaseScalar tol = 1.0e-14;
545 for (
size_t j=0;
j<ncol; ++
j) {
546 out <<
"dots(" <<
j <<
") = " << dots[
j]
547 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
548 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
563 using Teuchos::ArrayView;
564 using Teuchos::Array;
565 using Teuchos::ArrayRCP;
567 typedef typename Storage::value_type BaseScalar;
569 typedef typename Scalar::cijk_type Cijk;
571 typedef Teuchos::Comm<int> Tpetra_Comm;
572 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
573 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
574 typedef typename Tpetra_MultiVector::dot_type dot_type;
577 if ( !Kokkos::is_initialized() )
578 Kokkos::initialize();
582 setGlobalCijkTensor(cijk);
586 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
590 RCP<const Tpetra_Map> map =
591 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
593 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
594 const size_t num_my_row = myGIDs.size();
598 RCP<Tpetra_MultiVector> x1 = Tpetra::createMultiVector<Scalar>(map, ncol);
599 RCP<Tpetra_MultiVector> x2 = Tpetra::createMultiVector<Scalar>(map, ncol);
600 ArrayRCP< ArrayRCP<Scalar> > x1_view = x1->get2dViewNonConst();
601 ArrayRCP< ArrayRCP<Scalar> > x2_view = x2->get2dViewNonConst();
602 Scalar val1(cijk), val2(cijk);
603 for (
size_t i=0; i<num_my_row; ++i) {
605 for (
size_t j=0;
j<ncol; ++
j) {
608 generate_multi_vector_coefficient<BaseScalar,size_t>(
609 nrow, ncol, pce_size, row,
j, k);
610 val1.fastAccessCoeff(k) = v;
611 val2.fastAccessCoeff(k) = 0.12345 * v;
613 x1_view[
j][i] = val1;
614 x2_view[
j][i] = val2;
617 x1_view = Teuchos::null;
618 x2_view = Teuchos::null;
622 Teuchos::Array<size_t> cols(ncol_sub);
623 cols[0] = 4; cols[1] = 2;
624 RCP<const Tpetra_MultiVector> x1_sub = x1->subView(cols());
625 RCP<const Tpetra_MultiVector> x2_sub = x2->subView(cols());
628 Array<dot_type> dots(ncol_sub);
629 x1_sub->dot(*x2_sub, dots());
634 Array<dot_type> local_vals(ncol_sub, dot_type(0));
635 for (
size_t i=0; i<num_my_row; ++i) {
637 for (
size_t j=0;
j<ncol_sub; ++
j) {
639 BaseScalar v = generate_multi_vector_coefficient<BaseScalar,size_t>(
640 nrow, ncol, pce_size, row, cols[
j], k);
641 local_vals[
j] += 0.12345 * v * v;
647 Array<dot_type> vals(ncol_sub, dot_type(0));
648 Teuchos::reduceAll(*comm, Teuchos::REDUCE_SUM,
649 Teuchos::as<int>(ncol_sub), local_vals.getRawPtr(),
652 BaseScalar tol = 1.0e-14;
653 for (
size_t j=0;
j<ncol_sub; ++
j) {
654 out <<
"dots(" <<
j <<
") = " << dots[
j]
655 <<
" expected(" <<
j <<
") = " << vals[
j] << std::endl;
656 TEST_FLOATING_EQUALITY( dots[
j], vals[
j], tol );
671 using Teuchos::ArrayView;
672 using Teuchos::Array;
673 using Teuchos::ArrayRCP;
675 typedef typename Storage::value_type BaseScalar;
677 typedef typename Scalar::cijk_type Cijk;
679 typedef Teuchos::Comm<int> Tpetra_Comm;
680 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
681 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
682 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
683 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
686 if ( !Kokkos::is_initialized() )
687 Kokkos::initialize();
691 setGlobalCijkTensor(cijk);
696 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
697 RCP<const Tpetra_Map> map =
698 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
700 RCP<Tpetra_CrsGraph> graph =
701 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
702 Array<GlobalOrdinal> columnIndices(2);
703 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
704 const size_t num_my_row = myGIDs.size();
705 for (
size_t i=0; i<num_my_row; ++i) {
707 columnIndices[0] = row;
710 columnIndices[1] = row+1;
713 graph->insertGlobalIndices(row, columnIndices(0,ncol));
715 graph->fillComplete();
716 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
719 Array<Scalar> vals(2);
721 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
723 const size_t num_col = row == nrow - 1 ? 1 : 2;
724 for (
size_t local_col=0; local_col<num_col; ++local_col) {
726 columnIndices[local_col] = col;
728 val.fastAccessCoeff(k) =
729 generate_matrix_coefficient<BaseScalar,size_t>(
730 nrow, pce_size, row, col, k);
731 vals[local_col] =
val;
733 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
735 matrix->fillComplete();
738 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
739 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
740 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
743 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
744 nrow, pce_size, row,
j);
745 x_view[local_row] =
val;
747 x_view = Teuchos::null;
756 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
757 matrix->apply(*x, *y);
763 ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
764 BaseScalar tol = 1.0e-14;
765 typename Cijk::HostMirror host_cijk =
766 Kokkos::create_mirror_view(cijk);
767 Kokkos::deep_copy(host_cijk, cijk);
768 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
770 const size_t num_col = row == nrow - 1 ? 1 : 2;
772 for (
size_t local_col=0; local_col<num_col; ++local_col) {
776 const LocalOrdinal entry_beg = host_cijk.entry_begin(i);
779 for (
LocalOrdinal entry = entry_beg; entry < entry_end; ++entry) {
782 const BaseScalar a_j =
783 generate_matrix_coefficient<BaseScalar,size_t>(
784 nrow, pce_size, row, col,
j);
785 const BaseScalar a_k =
786 generate_matrix_coefficient<BaseScalar,size_t>(
787 nrow, pce_size, row, col, k);
788 const BaseScalar x_j =
789 generate_vector_coefficient<BaseScalar,size_t>(
790 nrow, pce_size, col,
j);
791 const BaseScalar x_k =
792 generate_vector_coefficient<BaseScalar,size_t>(
793 nrow, pce_size, col, k);
794 tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
796 val.fastAccessCoeff(i) += tmp;
799 TEST_EQUALITY( y_view[local_row].size(), pce_size );
802 val.fastAccessCoeff(i), tol );
817 using Teuchos::ArrayView;
818 using Teuchos::Array;
819 using Teuchos::ArrayRCP;
821 typedef typename Storage::value_type BaseScalar;
823 typedef typename Scalar::cijk_type Cijk;
825 typedef Teuchos::Comm<int> Tpetra_Comm;
826 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
827 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
828 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
829 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
832 if ( !Kokkos::is_initialized() )
833 Kokkos::initialize();
837 setGlobalCijkTensor(cijk);
842 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
843 RCP<const Tpetra_Map> map =
844 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
846 RCP<Tpetra_CrsGraph> graph =
847 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
848 Array<GlobalOrdinal> columnIndices(2);
849 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
850 const size_t num_my_row = myGIDs.size();
851 for (
size_t i=0; i<num_my_row; ++i) {
853 columnIndices[0] = row;
856 columnIndices[1] = row+1;
859 graph->insertGlobalIndices(row, columnIndices(0,ncol));
861 graph->fillComplete();
862 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
865 Array<Scalar> vals(2);
867 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
869 const size_t num_col = row == nrow - 1 ? 1 : 2;
870 for (
size_t local_col=0; local_col<num_col; ++local_col) {
872 columnIndices[local_col] = col;
874 val.fastAccessCoeff(k) =
875 generate_matrix_coefficient<BaseScalar,size_t>(
876 nrow, pce_size, row, col, k);
877 vals[local_col] =
val;
879 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
881 matrix->fillComplete();
885 RCP<Tpetra_MultiVector> x = Tpetra::createMultiVector<Scalar>(map, ncol);
886 ArrayRCP< ArrayRCP<Scalar> > x_view = x->get2dViewNonConst();
887 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
889 for (
size_t col=0; col<ncol; ++col) {
892 generate_multi_vector_coefficient<BaseScalar,size_t>(
893 nrow, ncol, pce_size, row, col, k);
894 val.fastAccessCoeff(k) = v;
896 x_view[col][local_row] =
val;
899 x_view = Teuchos::null;
908 RCP<Tpetra_MultiVector> y = Tpetra::createMultiVector<Scalar>(map, ncol);
909 matrix->apply(*x, *y);
915 ArrayRCP< ArrayRCP<Scalar> > y_view = y->get2dViewNonConst();
916 BaseScalar tol = 1.0e-14;
917 typename Cijk::HostMirror host_cijk =
918 Kokkos::create_mirror_view(cijk);
919 Kokkos::deep_copy(host_cijk, cijk);
920 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
922 for (
size_t xcol=0; xcol<ncol; ++xcol) {
923 const size_t num_col = row == nrow - 1 ? 1 : 2;
925 for (
size_t local_col=0; local_col<num_col; ++local_col) {
929 const LocalOrdinal entry_beg = host_cijk.entry_begin(i);
932 for (
LocalOrdinal entry = entry_beg; entry < entry_end; ++entry) {
935 const BaseScalar a_j =
936 generate_matrix_coefficient<BaseScalar,size_t>(
937 nrow, pce_size, row, col,
j);
938 const BaseScalar a_k =
939 generate_matrix_coefficient<BaseScalar,size_t>(
940 nrow, pce_size, row, col, k);
941 const BaseScalar x_j =
942 generate_multi_vector_coefficient<BaseScalar,size_t>(
943 nrow, ncol, pce_size, col, xcol,
j);
944 const BaseScalar x_k =
945 generate_multi_vector_coefficient<BaseScalar,size_t>(
946 nrow, ncol, pce_size, col, xcol, k);
947 tmp += host_cijk.value(entry) * ( a_j * x_k + a_k * x_j );
949 val.fastAccessCoeff(i) += tmp;
952 TEST_EQUALITY( y_view[xcol][local_row].size(), pce_size );
955 val.fastAccessCoeff(i), tol );
971 using Teuchos::ArrayView;
972 using Teuchos::Array;
973 using Teuchos::ArrayRCP;
975 typedef typename Storage::value_type BaseScalar;
977 typedef typename Scalar::cijk_type Cijk;
979 typedef Teuchos::Comm<int> Tpetra_Comm;
980 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
981 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
982 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
983 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
985 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
986 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
989 if ( !Kokkos::is_initialized() )
990 Kokkos::initialize();
994 setGlobalCijkTensor(cijk);
999 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1000 RCP<const Tpetra_Map> map =
1001 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1003 RCP<Tpetra_CrsGraph> graph =
1004 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
1005 Array<GlobalOrdinal> columnIndices(2);
1006 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1007 const size_t num_my_row = myGIDs.size();
1008 for (
size_t i=0; i<num_my_row; ++i) {
1010 columnIndices[0] = row;
1012 if (row != nrow-1) {
1013 columnIndices[1] = row+1;
1016 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1018 graph->fillComplete();
1019 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1022 Array<Scalar> vals(2);
1024 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1026 const size_t num_col = row == nrow - 1 ? 1 : 2;
1027 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1029 columnIndices[local_col] = col;
1031 val.fastAccessCoeff(k) =
1032 generate_matrix_coefficient<BaseScalar,size_t>(
1033 nrow, pce_size, row, col, k);
1034 vals[local_col] =
val;
1036 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1038 matrix->fillComplete();
1041 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1043 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1044 for (
size_t i=0; i<num_my_row; ++i) {
1047 val.fastAccessCoeff(
j) = generate_vector_coefficient<BaseScalar,size_t>(
1048 nrow, pce_size, row,
j);
1054 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
1055 matrix->apply(*x, *y);
1072 RCP<const Tpetra_Map> flat_x_map, flat_y_map;
1073 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1076 cijk_graph, pce_size);
1077 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1081 RCP<Tpetra_Vector> y2 = Tpetra::createVector<Scalar>(map);
1083 RCP<Flat_Tpetra_Vector> flat_x =
1085 RCP<Flat_Tpetra_Vector> flat_y =
1087 flat_matrix->apply(*flat_x, *flat_y);
1108 BaseScalar tol = 1.0e-14;
1109 ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
1110 ArrayRCP<Scalar> y2_view = y2->get1dViewNonConst();
1111 for (
size_t i=0; i<num_my_row; ++i) {
1112 TEST_EQUALITY( y_view[i].size(), pce_size );
1113 TEST_EQUALITY( y2_view[i].size(), pce_size );
1131 using Teuchos::ArrayView;
1132 using Teuchos::Array;
1133 using Teuchos::ArrayRCP;
1134 using Teuchos::ParameterList;
1136 typedef typename Storage::value_type BaseScalar;
1138 typedef typename Scalar::cijk_type Cijk;
1140 typedef Teuchos::Comm<int> Tpetra_Comm;
1141 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1142 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1143 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1144 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1147 if ( !Kokkos::is_initialized() )
1148 Kokkos::initialize();
1152 setGlobalCijkTensor(cijk);
1157 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1158 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1159 RCP<const Tpetra_Map> map =
1160 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1162 RCP<Tpetra_CrsGraph> graph =
1163 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
1164 Array<GlobalOrdinal> columnIndices(3);
1165 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1166 const size_t num_my_row = myGIDs.size();
1167 for (
size_t i=0; i<num_my_row; ++i) {
1169 if (row == 0 || row == nrow-1) {
1170 columnIndices[0] = row;
1171 graph->insertGlobalIndices(row, columnIndices(0,1));
1174 columnIndices[0] = row-1;
1175 columnIndices[1] = row;
1176 columnIndices[2] = row+1;
1177 graph->insertGlobalIndices(row, columnIndices(0,3));
1180 graph->fillComplete();
1181 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1184 Array<Scalar> vals(3);
1187 a_val.fastAccessCoeff(
j) =
1188 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1190 for (
size_t i=0; i<num_my_row; ++i) {
1192 if (row == 0 || row == nrow-1) {
1193 columnIndices[0] = row;
1195 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1198 columnIndices[0] = row-1;
1199 columnIndices[1] = row;
1200 columnIndices[2] = row+1;
1201 vals[0] =
Scalar(1.0) * a_val;
1202 vals[1] =
Scalar(-2.0) * a_val;
1203 vals[2] =
Scalar(1.0) * a_val;
1204 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1207 matrix->fillComplete();
1213 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1215 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1218 b_val.fastAccessCoeff(
j) =
1219 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1221 for (
size_t i=0; i<num_my_row; ++i) {
1223 if (row == 0 || row == nrow-1)
1226 b_view[i] = b_val * (h*h);
1234 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
1235 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1236 typedef typename BST::mag_type base_mag_type;
1237 typedef typename Tpetra_Vector::mag_type mag_type;
1238 base_mag_type btol = 1e-9;
1239 mag_type tol = btol;
1242 out.getOStream().get());
1243 TEST_EQUALITY_CONST( solved,
true );
1249 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1250 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1251 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1252 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1255 cijk_graph, pce_size);
1256 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1258 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1260 RCP<Flat_Tpetra_Vector> flat_x =
1262 RCP<Flat_Tpetra_Vector> flat_b =
1265 tol, max_its, out.getOStream().get());
1266 TEST_EQUALITY_CONST( solved_flat,
true );
1270 ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
1271 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1272 for (
size_t i=0; i<num_my_row; ++i) {
1273 TEST_EQUALITY( x_view[i].size(), pce_size );
1274 TEST_EQUALITY( x2_view[i].size(), pce_size );
1280 if (
j < v.size() && BST::abs(v.coeff(
j)) < btol)
1281 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1282 if (
j < v2.size() && BST::abs(v2.coeff(
j)) < btol)
1283 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1287 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
1295#if defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION)
1308 using Teuchos::ArrayView;
1309 using Teuchos::Array;
1310 using Teuchos::ArrayRCP;
1311 using Teuchos::ParameterList;
1312 using Teuchos::getParametersFromXmlFile;
1314 typedef typename Storage::value_type BaseScalar;
1316 typedef typename Scalar::cijk_type Cijk;
1318 typedef Teuchos::Comm<int> Tpetra_Comm;
1319 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1320 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1321 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1322 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1325 if ( !Kokkos::is_initialized() )
1326 Kokkos::initialize();
1335 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1336 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1337 RCP<const Tpetra_Map> map =
1338 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1340 RCP<Tpetra_CrsGraph> graph =
1341 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
1342 Array<GlobalOrdinal> columnIndices(3);
1343 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1344 const size_t num_my_row = myGIDs.size();
1345 for (
size_t i=0; i<num_my_row; ++i) {
1347 if (row == 0 || row == nrow-1) {
1348 columnIndices[0] = row;
1349 graph->insertGlobalIndices(row, columnIndices(0,1));
1352 columnIndices[0] = row-1;
1353 columnIndices[1] = row;
1354 columnIndices[2] = row+1;
1355 graph->insertGlobalIndices(row, columnIndices(0,3));
1358 graph->fillComplete();
1359 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1362 Array<Scalar> vals(3);
1365 a_val.fastAccessCoeff(
j) =
1366 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1368 for (
size_t i=0; i<num_my_row; ++i) {
1370 if (row == 0 || row == nrow-1) {
1371 columnIndices[0] = row;
1373 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1376 columnIndices[0] = row-1;
1377 columnIndices[1] = row;
1378 columnIndices[2] = row+1;
1379 vals[0] =
Scalar(1.0) * a_val;
1380 vals[1] =
Scalar(-2.0) * a_val;
1381 vals[2] =
Scalar(1.0) * a_val;
1382 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1385 matrix->fillComplete();
1388 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1390 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1393 b_val.fastAccessCoeff(
j) =
1394 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1396 for (
size_t i=0; i<num_my_row; ++i) {
1398 if (row == 0 || row == nrow-1)
1401 b_view[i] = b_val * (h*h);
1405 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1406 typedef Kokkos::Details::ArithTraits<BaseScalar> BST;
1407 typedef typename BST::mag_type base_mag_type;
1408 typedef typename Tpetra_Vector::mag_type mag_type;
1409 base_mag_type btol = 1e-9;
1410 mag_type tol = btol;
1414 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1415 RCP<ParameterList> muelu_params =
1416 getParametersFromXmlFile(
"muelu_cheby.xml");
1417#if USE_SCALAR_MEAN_BASED_PREC
1418 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_OP;
1419 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1420 RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1422 RCP<Scalar_OP> mean_matrix_op = mean_matrix;
1423 RCP<Scalar_OP> M_s =
1424 MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1428 Stokhos::create_mean_based_product_tensor<typename Storage::execution_space,typename Storage::ordinal_type,BaseScalar>();
1431 RCP<OP> mean_matrix_op = mean_matrix;
1433 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
1439 out.getOStream().get());
1440 TEST_EQUALITY_CONST( solved,
true );
1447 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1449 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1450 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1451 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1452 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1455 cijk_graph, pce_size);
1456 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1458 RCP<Flat_Tpetra_Vector> flat_x =
1460 RCP<Flat_Tpetra_Vector> flat_b =
1468 tol, max_its, out.getOStream().get());
1469 TEST_EQUALITY_CONST( solved_flat,
true );
1473 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1474 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1475 for (
size_t i=0; i<num_my_row; ++i) {
1476 TEST_EQUALITY( x_view[i].size(), pce_size );
1477 TEST_EQUALITY( x2_view[i].size(), pce_size );
1483 if (
j < v.size() && BST::abs(v.coeff(
j)) < btol)
1484 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1485 if (
j < v2.size() && BST::abs(v2.coeff(
j)) < btol)
1486 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1490 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
1505#if defined(HAVE_STOKHOS_BELOS)
1515 using Teuchos::ArrayView;
1516 using Teuchos::Array;
1517 using Teuchos::ArrayRCP;
1518 using Teuchos::ParameterList;
1520 typedef typename Storage::value_type BaseScalar;
1522 typedef typename Scalar::cijk_type Cijk;
1524 typedef Teuchos::Comm<int> Tpetra_Comm;
1525 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1526 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1527 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1528 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1531 if ( !Kokkos::is_initialized() )
1532 Kokkos::initialize();
1541 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1542 RCP<const Tpetra_Map> map =
1543 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1545 RCP<Tpetra_CrsGraph> graph =
1546 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
1547 Array<GlobalOrdinal> columnIndices(2);
1548 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1549 const size_t num_my_row = myGIDs.size();
1550 for (
size_t i=0; i<num_my_row; ++i) {
1552 columnIndices[0] = row;
1554 if (row != nrow-1) {
1555 columnIndices[1] = row+1;
1558 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1560 graph->fillComplete();
1561 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1564 Array<Scalar> vals(2);
1566 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1568 const size_t num_col = row == nrow - 1 ? 1 : 2;
1569 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1571 columnIndices[local_col] = col;
1573 val.fastAccessCoeff(k) =
1574 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1575 vals[local_col] =
val;
1577 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1579 matrix->fillComplete();
1582 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1584 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1585 for (
size_t i=0; i<num_my_row; ++i) {
1591 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1592 typedef BaseScalar BelosScalar;
1593 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1594 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1595 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1596 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1597 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
1598 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1599 typename ST::magnitudeType tol = 1e-9;
1600 belosParams->set(
"Flexible Gmres",
false);
1601 belosParams->set(
"Num Blocks", 100);
1602 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1603 belosParams->set(
"Maximum Iterations", 100);
1604 belosParams->set(
"Verbosity", 33);
1605 belosParams->set(
"Output Style", 1);
1606 belosParams->set(
"Output Frequency", 1);
1607 belosParams->set(
"Output Stream", out.getOStream());
1608 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1609 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1610 problem->setProblem();
1611 Belos::ReturnType ret = solver->solve();
1612 TEST_EQUALITY_CONST( ret, Belos::Converged );
1618 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1619 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1620 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1621 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1624 cijk_graph, pce_size);
1625 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1627 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1629 RCP<Flat_Tpetra_Vector> flat_x =
1631 RCP<Flat_Tpetra_Vector> flat_b =
1633 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
1634 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1635 typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1636 RCP< FBLinProb > flat_problem =
1637 rcp(
new FBLinProb(flat_matrix, flat_x, flat_b));
1638 RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1639 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1641 flat_problem->setProblem();
1642 Belos::ReturnType flat_ret = flat_solver->solve();
1643 TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
1646 typename ST::magnitudeType btol = 100*tol;
1647 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1648 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1649 for (
size_t i=0; i<num_my_row; ++i) {
1650 TEST_EQUALITY( x_view[i].size(), pce_size );
1651 TEST_EQUALITY( x2_view[i].size(), pce_size );
1657 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
1658 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1659 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
1660 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1664 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
1681#if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2)
1692 using Teuchos::ArrayView;
1693 using Teuchos::Array;
1694 using Teuchos::ArrayRCP;
1695 using Teuchos::ParameterList;
1697 typedef typename Storage::value_type BaseScalar;
1699 typedef typename Scalar::cijk_type Cijk;
1701 typedef Teuchos::Comm<int> Tpetra_Comm;
1702 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1703 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1704 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1705 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1708 if ( !Kokkos::is_initialized() )
1709 Kokkos::initialize();
1718 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1719 RCP<const Tpetra_Map> map =
1720 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1722 RCP<Tpetra_CrsGraph> graph =
1723 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
1724 Array<GlobalOrdinal> columnIndices(2);
1725 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1726 const size_t num_my_row = myGIDs.size();
1727 for (
size_t i=0; i<num_my_row; ++i) {
1729 columnIndices[0] = row;
1731 if (row != nrow-1) {
1732 columnIndices[1] = row+1;
1735 graph->insertGlobalIndices(row, columnIndices(0,ncol));
1737 graph->fillComplete();
1738 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1741 Array<Scalar> vals(2);
1743 for (
size_t local_row=0; local_row<num_my_row; ++local_row) {
1745 const size_t num_col = row == nrow - 1 ? 1 : 2;
1746 for (
size_t local_col=0; local_col<num_col; ++local_col) {
1748 columnIndices[local_col] = col;
1750 val.fastAccessCoeff(k) =
1751 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(k+1);
1752 vals[local_col] =
val;
1754 matrix->replaceGlobalValues(row, columnIndices(0,num_col), vals(0,num_col));
1756 matrix->fillComplete();
1759 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1761 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1762 for (
size_t i=0; i<num_my_row; ++i) {
1767 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
1768 RCP<ParameterList> belosParams = rcp(
new ParameterList);
1769 typedef BaseScalar BelosScalar;
1770 typedef Teuchos::ScalarTraits<BaseScalar> ST;
1771 typename ST::magnitudeType tol = 1e-9;
1773 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
1774 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1775 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
1777#if USE_SCALAR_MEAN_BASED_PREC
1778 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1779 RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1781 typedef Ifpack2::Preconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Prec;
1782 Ifpack2::Factory factory;
1783 RCP<Scalar_Prec> M_s = factory.create<Scalar_Tpetra_CrsMatrix>(
"RILUK", mean_matrix);
1789 typedef Ifpack2::Preconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node> Prec;
1790 Ifpack2::Factory factory;
1791 RCP<Prec> M = factory.create<Tpetra_CrsMatrix>(
"RILUK", mean_matrix);
1797 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
1798 problem->setRightPrec(M);
1799 problem->setProblem();
1800 belosParams->set(
"Flexible Gmres",
false);
1801 belosParams->set(
"Num Blocks", 100);
1802 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
1803 belosParams->set(
"Maximum Iterations", 100);
1804 belosParams->set(
"Verbosity", 33);
1805 belosParams->set(
"Output Style", 1);
1806 belosParams->set(
"Output Frequency", 1);
1807 belosParams->set(
"Output Stream", out.getOStream());
1809 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
1810 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
1811 Belos::ReturnType ret = solver->solve();
1812 TEST_EQUALITY_CONST( ret, Belos::Converged );
1816 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
1818 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
1819 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
1820 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
1821 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
1824 cijk_graph, pce_size);
1825 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
1827 RCP<Flat_Tpetra_Vector> flat_x =
1829 RCP<Flat_Tpetra_Vector> flat_b =
1831 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
1832 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FOP;
1833 typedef Belos::LinearProblem<BelosScalar,FMV,FOP> FBLinProb;
1834 RCP< FBLinProb > flat_problem =
1835 rcp(
new FBLinProb(flat_matrix, flat_x, flat_b));
1836 RCP<Belos::SolverManager<BelosScalar,FMV,FOP> > flat_solver =
1837 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,FMV,FOP>(flat_problem,
1839 flat_problem->setProblem();
1840 Belos::ReturnType flat_ret = flat_solver->solve();
1841 TEST_EQUALITY_CONST( flat_ret, Belos::Converged );
1844 typename ST::magnitudeType btol = 100*tol;
1845 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
1846 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
1847 for (
size_t i=0; i<num_my_row; ++i) {
1848 TEST_EQUALITY( x_view[i].size(), pce_size );
1849 TEST_EQUALITY( x2_view[i].size(), pce_size );
1855 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
1856 v.fastAccessCoeff(
j) = BaseScalar(0.0);
1857 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
1858 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
1862 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
1877#if defined(HAVE_STOKHOS_BELOS) && defined(HAVE_STOKHOS_IFPACK2) && defined(HAVE_STOKHOS_MUELU) && defined(HAVE_STOKHOS_AMESOS2) && defined(HAVE_TPETRA_EXPLICIT_INSTANTIATION)
1890 using Teuchos::ArrayView;
1891 using Teuchos::Array;
1892 using Teuchos::ArrayRCP;
1893 using Teuchos::ParameterList;
1894 using Teuchos::getParametersFromXmlFile;
1896 typedef typename Storage::value_type BaseScalar;
1898 typedef typename Scalar::cijk_type Cijk;
1900 typedef Teuchos::Comm<int> Tpetra_Comm;
1901 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
1902 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
1903 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
1904 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
1907 if ( !Kokkos::is_initialized() )
1908 Kokkos::initialize();
1917 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
1918 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
1919 RCP<const Tpetra_Map> map =
1920 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
1922 RCP<Tpetra_CrsGraph> graph =
1923 rcp(
new Tpetra_CrsGraph(map,
size_t(3)));
1924 Array<GlobalOrdinal> columnIndices(3);
1925 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
1926 const size_t num_my_row = myGIDs.size();
1927 for (
size_t i=0; i<num_my_row; ++i) {
1929 if (row == 0 || row == nrow-1) {
1930 columnIndices[0] = row;
1931 graph->insertGlobalIndices(row, columnIndices(0,1));
1934 columnIndices[0] = row-1;
1935 columnIndices[1] = row;
1936 columnIndices[2] = row+1;
1937 graph->insertGlobalIndices(row, columnIndices(0,3));
1940 graph->fillComplete();
1941 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
1944 Array<Scalar> vals(3);
1947 a_val.fastAccessCoeff(
j) =
1948 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
1950 for (
size_t i=0; i<num_my_row; ++i) {
1952 if (row == 0 || row == nrow-1) {
1953 columnIndices[0] = row;
1955 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
1958 columnIndices[0] = row-1;
1959 columnIndices[1] = row;
1960 columnIndices[2] = row+1;
1961 vals[0] =
Scalar(1.0) * a_val;
1962 vals[1] =
Scalar(-2.0) * a_val;
1963 vals[2] =
Scalar(1.0) * a_val;
1964 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
1967 matrix->fillComplete();
1970 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
1972 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
1975 b_val.fastAccessCoeff(
j) =
1976 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
1978 for (
size_t i=0; i<num_my_row; ++i) {
1980 if (row == 0 || row == nrow-1)
1983 b_view[i] = b_val * (h*h);
1988 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
1989 RCP<ParameterList> muelu_params =
1990 getParametersFromXmlFile(
"muelu_cheby.xml");
1991#if USE_SCALAR_MEAN_BASED_PREC
1992 typedef Tpetra::Operator<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_OP;
1993 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Scalar_Tpetra_CrsMatrix;
1994 RCP<Scalar_Tpetra_CrsMatrix> mean_matrix =
1996 RCP<Scalar_OP> mean_matrix_op = mean_matrix;
1997 RCP<Scalar_OP> M_s =
1998 MueLu::CreateTpetraPreconditioner<BaseScalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
2002 Stokhos::create_mean_based_product_tensor<typename Storage::execution_space,typename Storage::ordinal_type,BaseScalar>();
2005 RCP<OP> mean_matrix_op = mean_matrix;
2007 MueLu::CreateTpetraPreconditioner<Scalar,LocalOrdinal,GlobalOrdinal,Node>(mean_matrix_op, *muelu_params);
2012 typedef Teuchos::ScalarTraits<BaseScalar> ST;
2013 typedef BaseScalar BelosScalar;
2014 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
2015 typedef Belos::LinearProblem<BelosScalar,MV,OP> BLinProb;
2016 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2017 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
2018 problem->setRightPrec(M);
2019 problem->setProblem();
2020 RCP<ParameterList> belosParams = rcp(
new ParameterList);
2021 typename ST::magnitudeType tol = 1e-9;
2022 belosParams->set(
"Num Blocks", 100);
2023 belosParams->set(
"Convergence Tolerance", BelosScalar(tol));
2024 belosParams->set(
"Maximum Iterations", 100);
2025 belosParams->set(
"Verbosity", 33);
2026 belosParams->set(
"Output Style", 1);
2027 belosParams->set(
"Output Frequency", 1);
2028 belosParams->set(
"Output Stream", out.getOStream());
2030 RCP<Belos::SolverManager<BelosScalar,MV,OP> > solver =
2031 rcp(
new Belos::PseudoBlockGmresSolMgr<BelosScalar,MV,OP>(problem, belosParams));
2034 Belos::ReturnType ret = solver->solve();
2035 TEST_EQUALITY_CONST( ret, Belos::Converged );
2041 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
2042 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2043 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2044 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2047 cijk_graph, pce_size);
2048 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2050 RCP<Tpetra_Vector> x2 = Tpetra::createVector<Scalar>(map);
2052 RCP<Flat_Tpetra_Vector> flat_x =
2054 RCP<Flat_Tpetra_Vector> flat_b =
2056 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> FMV;
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 );
2072 typename ST::magnitudeType btol = 100*tol;
2073 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
2074 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2075 for (
size_t i=0; i<num_my_row; ++i) {
2076 TEST_EQUALITY( x_view[i].size(), pce_size );
2077 TEST_EQUALITY( x2_view[i].size(), pce_size );
2083 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
2084 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2085 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
2086 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
2090 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
2106#if defined(HAVE_STOKHOS_AMESOS2)
2116 using Teuchos::ArrayView;
2117 using Teuchos::Array;
2118 using Teuchos::ArrayRCP;
2119 using Teuchos::ParameterList;
2121 typedef typename Storage::value_type BaseScalar;
2123 typedef typename Scalar::cijk_type Cijk;
2125 typedef Teuchos::Comm<int> Tpetra_Comm;
2126 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
2127 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
2128 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_MultiVector;
2129 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
2130 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
2133 if ( !Kokkos::is_initialized() )
2134 Kokkos::initialize();
2143 BaseScalar h = 1.0 /
static_cast<BaseScalar
>(nrow-1);
2144 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
2145 RCP<const Tpetra_Map> map =
2146 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
2148 RCP<Tpetra_CrsGraph> graph = Tpetra::createCrsGraph(map,
size_t(3));
2149 Array<GlobalOrdinal> columnIndices(3);
2150 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
2151 const size_t num_my_row = myGIDs.size();
2152 for (
size_t i=0; i<num_my_row; ++i) {
2154 if (row == 0 || row == nrow-1) {
2155 columnIndices[0] = row;
2156 graph->insertGlobalIndices(row, columnIndices(0,1));
2159 columnIndices[0] = row-1;
2160 columnIndices[1] = row;
2161 columnIndices[2] = row+1;
2162 graph->insertGlobalIndices(row, columnIndices(0,3));
2165 graph->fillComplete();
2166 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
2169 Array<Scalar> vals(3);
2172 a_val.fastAccessCoeff(
j) =
2173 BaseScalar(1.0) + BaseScalar(1.0) / BaseScalar(
j+1);
2175 for (
size_t i=0; i<num_my_row; ++i) {
2177 if (row == 0 || row == nrow-1) {
2178 columnIndices[0] = row;
2180 matrix->replaceGlobalValues(row, columnIndices(0,1), vals(0,1));
2183 columnIndices[0] = row-1;
2184 columnIndices[1] = row;
2185 columnIndices[2] = row+1;
2186 vals[0] =
Scalar(1.0) * a_val;
2187 vals[1] =
Scalar(-2.0) * a_val;
2188 vals[2] =
Scalar(1.0) * a_val;
2189 matrix->replaceGlobalValues(row, columnIndices(0,3), vals(0,3));
2192 matrix->fillComplete();
2195 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
2197 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
2200 b_val.fastAccessCoeff(
j) =
2201 BaseScalar(2.0) - BaseScalar(1.0) / BaseScalar(
j+1);
2203 for (
size_t i=0; i<num_my_row; ++i) {
2205 if (row == 0 || row == nrow-1)
2208 b_view[i] = b_val * (h*h);
2213 typedef Amesos2::Solver<Tpetra_CrsMatrix,Tpetra_MultiVector> Solver;
2214 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
2215 std::string solver_name;
2216#if defined(HAVE_AMESOS2_BASKER)
2217 solver_name =
"basker";
2218#elif defined(HAVE_AMESOS2_KLU2)
2219 solver_name =
"klu2";
2220#elif defined(HAVE_AMESOS2_SUPERLUDIST)
2221 solver_name =
"superlu_dist";
2222#elif defined(HAVE_AMESOS2_SUPERLUMT)
2223 solver_name =
"superlu_mt";
2224#elif defined(HAVE_AMESOS2_SUPERLU)
2225 solver_name =
"superlu";
2226#elif defined(HAVE_AMESOS2_PARDISO_MKL)
2227 solver_name =
"pardisomkl";
2228#elif defined(HAVE_AMESOS2_LAPACK)
2229 solver_name =
"lapack";
2230#elif defined(HAVE_AMESOS2_CHOLMOD) && defined (HAVE_AMESOS2_EXPERIMENTAL)
2231 solver_name =
"lapack";
2237 out <<
"Solving linear system with " << solver_name << std::endl;
2239 RCP<Solver> solver = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(
2240 solver_name, matrix, x, b);
2247 typedef Tpetra::Vector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_Vector;
2248 typedef Tpetra::MultiVector<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_MultiVector;
2249 typedef Tpetra::CrsMatrix<BaseScalar,LocalOrdinal,GlobalOrdinal,Node> Flat_Tpetra_CrsMatrix;
2250 RCP<const Tpetra_Map> flat_x_map, flat_b_map;
2251 RCP<const Tpetra_CrsGraph> flat_graph, cijk_graph;
2254 cijk_graph, pce_size);
2255 RCP<Flat_Tpetra_CrsMatrix> flat_matrix =
2257 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();
2270 typedef Kokkos::Details::ArithTraits<BaseScalar> ST;
2271 typename ST::mag_type btol = 1e-12;
2272 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
2273 ArrayRCP<Scalar> x2_view = x2->get1dViewNonConst();
2274 for (
size_t i=0; i<num_my_row; ++i) {
2275 TEST_EQUALITY( x_view[i].size(), pce_size );
2276 TEST_EQUALITY( x2_view[i].size(), pce_size );
2282 if (
j < v.size() && ST::magnitude(v.coeff(
j)) < btol)
2283 v.fastAccessCoeff(
j) = BaseScalar(0.0);
2284 if (
j < v2.size() && ST::magnitude(v2.coeff(
j)) < btol)
2285 v2.fastAccessCoeff(
j) = BaseScalar(0.0);
2289 TEST_FLOATING_EQUALITY(v.coeff(
j), v2.coeff(
j), btol);
2304#define CRSMATRIX_UQ_PCE_TESTS_SLGN(S, LO, GO, N) \
2305 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorAdd, S, LO, GO, N ) \
2306 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, VectorDot, S, LO, GO, N ) \
2307 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorAdd, S, LO, GO, N ) \
2308 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDot, S, LO, GO, N ) \
2309 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MultiVectorDotSub, S, LO, GO, N ) \
2310 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixVectorMultiply, S, LO, GO, N ) \
2311 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, MatrixMultiVectorMultiply, S, LO, GO, N ) \
2312 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Flatten, S, LO, GO, N ) \
2313 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimpleCG, S, LO, GO, N ) \
2314 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, SimplePCG_Muelu, S, LO, GO, N ) \
2315 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES, S, LO, GO, N ) \
2316 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosGMRES_RILUK, S, LO, GO, N ) \
2317 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, BelosCG_Muelu, S, LO, GO, N ) \
2318 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix_PCE, Amesos2, S, LO, GO, N )
2320#define CRSMATRIX_UQ_PCE_TESTS_N(N) \
2321 typedef Stokhos::DeviceForNode2<N>::type Device; \
2322 typedef Stokhos::DynamicStorage<int,double,Device::execution_space> DS; \
2323 using default_local_ordinal_type = Tpetra::Map<>::local_ordinal_type; \
2324 using default_global_ordinal_type = Tpetra::Map<>::global_ordinal_type; \
2325 CRSMATRIX_UQ_PCE_TESTS_SLGN(DS, default_local_ordinal_type, default_global_ordinal_type, N)
expr1 expr1 expr1 expr2 expr1 expr1 c expr2 expr1 c fastAccessCoeff(j) - expr2.val(j)
Kokkos::DefaultExecutionSpace execution_space
kokkos_cijk_type build_cijk(ordinal_type stoch_dim, ordinal_type poly_ord)
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_PCE, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node)
scalar generate_multi_vector_coefficient(const ordinal nFEM, const ordinal nVec, const ordinal nStoch, const ordinal iColFEM, const ordinal iVec, const ordinal iStoch)
scalar generate_vector_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iColFEM, const ordinal iStoch)
scalar generate_matrix_coefficient(const ordinal nFEM, const ordinal nStoch, const ordinal iRowFEM, const ordinal iColFEM, const ordinal iStoch)
Multivariate orthogonal polynomial basis generated from a total-order complete-polynomial tensor prod...
Legendre polynomial basis.
Abstract base class for 1-D orthogonal polynomials.
Data structure storing a sparse 3-tensor C(i,j,k) in a a compressed format.
KokkosClassic::DefaultNode::DefaultNodeType Node
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< is_view_uq_pce< view_type >::value, typenameCijkType< view_type >::type >::type cijk(const view_type &view)
void setGlobalCijkTensor(const cijk_type &cijk)
Teuchos::RCP< Tpetra::CrsMatrix< typename Scalar::value_type, LO, GO, N > > build_mean_scalar_matrix(const Tpetra::CrsMatrix< Scalar, LO, GO, N > &A)
Teuchos::RCP< Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_pce_graph(const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &graph, const CijkType &cijk, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_domain_map, Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_range_map, Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &cijk_graph, const size_t matrix_pce_size)
bool CG_Solve(const Matrix &A, Vector &x, const Vector &b, typename Vector::mag_type tol, Ordinal max_its, std::ostream *out=0)
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::CrsMatrix< Scalar, LO, GO, N > > build_mean_matrix(const Tpetra::CrsMatrix< Scalar, LO, GO, N > &A)
Teuchos::RCP< const Tpetra::MultiVector< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_vector_view(const Tpetra::MultiVector< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &vec, const Teuchos::RCP< const Tpetra::Map< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_map)
Teuchos::RCP< Tpetra::CrsMatrix< typename Storage::value_type, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > create_flat_matrix(const Tpetra::CrsMatrix< Sacado::UQ::PCE< Storage >, LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > &mat, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &flat_graph, const Teuchos::RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Kokkos::Compat::KokkosDeviceWrapperNode< Device > > > &cijk_graph, const CijkType &cijk_dev)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x