2 #include "Teuchos_UnitTestHarness.hpp" 3 #include "Teuchos_UnitTestRepository.hpp" 4 #include "Teuchos_UnitTestHelpers.hpp" 5 #include "Teuchos_GlobalMPISession.hpp" 8 #include "Tpetra_Core.hpp" 9 #include "Tpetra_Map.hpp" 10 #include "Tpetra_Vector.hpp" 11 #include "Tpetra_CrsGraph.hpp" 12 #include "Tpetra_CrsMatrix.hpp" 16 #ifdef HAVE_STOKHOS_BELOS 17 #include "BelosTpetraAdapter.hpp" 18 #include "BelosLinearProblem.hpp" 19 #include "BelosPseudoBlockGmresSolMgr.hpp" 27 using Teuchos::ArrayView;
29 using Teuchos::ArrayRCP;
33 typedef Teuchos::Comm<int> Tpetra_Comm;
34 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
35 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
36 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
37 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
41 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
42 RCP<const Tpetra_Map> map =
43 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
45 RCP<Tpetra_CrsGraph> graph =
46 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
47 Array<GlobalOrdinal> columnIndices(2);
48 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
49 const size_t num_my_row = myGIDs.size();
50 for (
size_t i=0; i<num_my_row; ++i) {
52 columnIndices[0] = row;
59 graph->insertGlobalIndices(row, columnIndices(0,ncol));
61 graph->fillComplete();
62 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
65 Array<Scalar> vals(2);
66 for (
size_t i=0; i<num_my_row; ++i) {
68 columnIndices[0] = row;
77 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
79 matrix->fillComplete();
82 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
85 matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
86 Teuchos::VERB_EXTREME);
88 x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
89 Teuchos::VERB_EXTREME);
92 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
93 matrix->apply(*
x, *
y);
95 y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
96 Teuchos::VERB_EXTREME);
99 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
100 ArrayRCP<Scalar> x_view =
y->get1dViewNonConst();
101 for (
size_t i=0; i<num_my_row; ++i) {
108 TEST_EQUALITY( y_view[i],
val );
112 #if defined(HAVE_STOKHOS_BELOS) 122 using Teuchos::ArrayView;
123 using Teuchos::Array;
124 using Teuchos::ArrayRCP;
125 using Teuchos::ParameterList;
129 typedef Teuchos::Comm<int> Tpetra_Comm;
130 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
131 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
132 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
133 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
137 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
138 RCP<const Tpetra_Map> map =
139 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
141 RCP<Tpetra_CrsGraph> graph =
142 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
143 Array<GlobalOrdinal> columnIndices(2);
144 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
145 const size_t num_my_row = myGIDs.size();
146 for (
size_t i=0; i<num_my_row; ++i) {
148 columnIndices[0] = row;
154 graph->insertGlobalIndices(row, columnIndices(0,ncol));
156 graph->fillComplete();
157 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
160 Array<Scalar> vals(2);
161 for (
size_t i=0; i<num_my_row; ++i) {
163 columnIndices[0] = row;
172 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
174 matrix->fillComplete();
177 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
178 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
179 for (
size_t i=0; i<num_my_row; ++i) {
184 typedef Teuchos::ScalarTraits<Scalar> ST;
185 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
186 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
187 typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
188 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
189 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
190 RCP<ParameterList> belosParams = rcp(
new ParameterList);
191 typename ST::magnitudeType tol = 1e-9;
192 belosParams->set(
"Flexible Gmres",
false);
193 belosParams->set(
"Num Blocks", 100);
194 belosParams->set(
"Convergence Tolerance",
Scalar(tol));
195 belosParams->set(
"Maximum Iterations", 100);
196 belosParams->set(
"Verbosity", 33);
197 belosParams->set(
"Output Style", 1);
198 belosParams->set(
"Output Frequency", 1);
199 belosParams->set(
"Output Stream", out.getOStream());
200 RCP<Belos::SolverManager<Scalar,MV,OP> > solver =
201 rcp(
new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem, belosParams));
202 problem->setProblem();
203 Belos::ReturnType ret = solver->solve();
204 TEST_EQUALITY_CONST( ret, Belos::Converged );
215 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
217 for (
size_t i=0; i<num_my_row; ++i) {
227 if (ST::magnitude(v) < tol)
230 TEST_FLOATING_EQUALITY(v,
val, tol);
245 Tpetra_CrsMatrix, MatVec,
double,
Node )
250 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
253 int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
Node int main(int argc, char *argv[])
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
KokkosClassic::DefaultNode::DefaultNodeType Node
TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(Tpetra_CrsMatrix, MatVec, double, Node) TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(Tpetra_CrsMatrix
TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(Tpetra_CrsMatrix, MatVec, Scalar, Node)
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y