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)));
47 Array<GlobalOrdinal> columnIndices(2);
48 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
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 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
100 for (
size_t i=0; i<num_my_row; ++i) {
107 TEST_EQUALITY( y_view(i,0),
val );
111#if defined(HAVE_STOKHOS_BELOS)
121 using Teuchos::ArrayView;
122 using Teuchos::Array;
123 using Teuchos::ArrayRCP;
124 using Teuchos::ParameterList;
128 typedef Teuchos::Comm<int> Tpetra_Comm;
129 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
130 typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
131 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
132 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
136 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
137 RCP<const Tpetra_Map> map =
138 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
140 RCP<Tpetra_CrsGraph> graph =
141 rcp(
new Tpetra_CrsGraph(map,
size_t(2)));
142 Array<GlobalOrdinal> columnIndices(2);
143 ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
144 const size_t num_my_row = myGIDs.size();
145 for (
size_t i=0; i<num_my_row; ++i) {
147 columnIndices[0] = row;
153 graph->insertGlobalIndices(row, columnIndices(0,ncol));
155 graph->fillComplete();
156 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
159 Array<Scalar> vals(2);
160 for (
size_t i=0; i<num_my_row; ++i) {
162 columnIndices[0] = row;
171 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
173 matrix->fillComplete();
176 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
178 auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
179 for (
size_t i=0; i<num_my_row; ++i) {
180 b_view(i, 0) =
Scalar(1.0);
185 typedef Teuchos::ScalarTraits<Scalar> ST;
186 typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
187 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
188 typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
189 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
190 RCP< BLinProb > problem = rcp(
new BLinProb(matrix, x, b));
191 RCP<ParameterList> belosParams = rcp(
new ParameterList);
192 typename ST::magnitudeType tol = 1e-9;
193 belosParams->set(
"Flexible Gmres",
false);
194 belosParams->set(
"Num Blocks", 100);
195 belosParams->set(
"Convergence Tolerance",
Scalar(tol));
196 belosParams->set(
"Maximum Iterations", 100);
197 belosParams->set(
"Verbosity", 33);
198 belosParams->set(
"Output Style", 1);
199 belosParams->set(
"Output Frequency", 1);
200 belosParams->set(
"Output Stream", out.getOStream());
201 RCP<Belos::SolverManager<Scalar,MV,OP> > solver =
202 rcp(
new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem, belosParams));
203 problem->setProblem();
204 Belos::ReturnType ret = solver->solve();
205 TEST_EQUALITY_CONST( ret, Belos::Converged );
216 auto x_view =
x->getLocalViewHost(Tpetra::Access::ReadOnly);
218 for (
size_t i=0; i<num_my_row; ++i) {
228 if (ST::magnitude(v) < tol)
231 TEST_FLOATING_EQUALITY(v,
val, tol);
243using Node = Tpetra::Map<>::node_type;
246 Tpetra_CrsMatrix, MatVec,
double,
Node )
251 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
254 int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
KokkosClassic::DefaultNode::DefaultNodeType Node
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(Tpetra_CrsMatrix, MatVec, double, Node) TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(Tpetra_CrsMatrix
Node int main(int argc, char *argv[])
TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(Tpetra_CrsMatrix, MatVec, Scalar, Node)