2 #include "Teuchos_UnitTestHarness.hpp" 3 #include "Teuchos_UnitTestRepository.hpp" 4 #include "Teuchos_UnitTestHelpers.hpp" 5 #include "Teuchos_GlobalMPISession.hpp" 8 #include "Tpetra_DefaultPlatform.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;
31 typedef Teuchos::Comm<int> Tpetra_Comm;
32 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
34 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
35 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
39 RCP<const Tpetra_Comm> comm =
40 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
41 RCP<const Tpetra_Map> map =
42 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
44 RCP<Tpetra_CrsGraph> graph =
45 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
46 Array<GlobalOrdinal> columnIndices(2);
47 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
48 const size_t num_my_row = myGIDs.size();
49 for (
size_t i=0; i<num_my_row; ++i) {
51 columnIndices[0] = row;
58 graph->insertGlobalIndices(row, columnIndices(0,ncol));
60 graph->fillComplete();
61 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
64 Array<Scalar> vals(2);
65 for (
size_t i=0; i<num_my_row; ++i) {
67 columnIndices[0] = row;
76 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
78 matrix->fillComplete();
81 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
84 matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
85 Teuchos::VERB_EXTREME);
87 x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
88 Teuchos::VERB_EXTREME);
91 RCP<Tpetra_Vector>
y = Tpetra::createVector<Scalar>(map);
92 matrix->apply(*
x, *
y);
94 y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,
false))),
95 Teuchos::VERB_EXTREME);
98 ArrayRCP<Scalar> y_view =
y->get1dViewNonConst();
99 ArrayRCP<Scalar> x_view =
y->get1dViewNonConst();
100 for (
size_t i=0; i<num_my_row; ++i) {
107 TEST_EQUALITY( y_view[i],
val );
111 #if defined(HAVE_STOKHOS_BELOS) 121 using Teuchos::ArrayView;
122 using Teuchos::Array;
123 using Teuchos::ArrayRCP;
124 using Teuchos::ParameterList;
126 typedef Teuchos::Comm<int> Tpetra_Comm;
127 typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
129 typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
130 typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
134 RCP<const Tpetra_Comm> comm =
135 Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
136 RCP<const Tpetra_Map> map =
137 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
139 RCP<Tpetra_CrsGraph> graph =
140 rcp(
new Tpetra_CrsGraph(map,
size_t(2), Tpetra::StaticProfile));
141 Array<GlobalOrdinal> columnIndices(2);
142 ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
143 const size_t num_my_row = myGIDs.size();
144 for (
size_t i=0; i<num_my_row; ++i) {
146 columnIndices[0] = row;
152 graph->insertGlobalIndices(row, columnIndices(0,ncol));
154 graph->fillComplete();
155 RCP<Tpetra_CrsMatrix> matrix = rcp(
new Tpetra_CrsMatrix(graph));
158 Array<Scalar> vals(2);
159 for (
size_t i=0; i<num_my_row; ++i) {
161 columnIndices[0] = row;
170 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
172 matrix->fillComplete();
175 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
176 ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
177 for (
size_t i=0; i<num_my_row; ++i) {
182 typedef Teuchos::ScalarTraits<Scalar> ST;
184 typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
185 typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
186 RCP<Tpetra_Vector>
x = Tpetra::createVector<Scalar>(map);
187 RCP< BLinProb > problem = rcp(
new BLinProb(matrix,
x, b));
188 RCP<ParameterList> belosParams = rcp(
new ParameterList);
189 typename ST::magnitudeType tol = 1e-9;
190 belosParams->set(
"Flexible Gmres",
false);
191 belosParams->set(
"Num Blocks", 100);
192 belosParams->set(
"Convergence Tolerance",
Scalar(tol));
193 belosParams->set(
"Maximum Iterations", 100);
194 belosParams->set(
"Verbosity", 33);
195 belosParams->set(
"Output Style", 1);
196 belosParams->set(
"Output Frequency", 1);
197 belosParams->set(
"Output Stream", out.getOStream());
198 RCP<Belos::SolverManager<Scalar,MV,OP> > solver =
199 rcp(
new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem, belosParams));
200 problem->setProblem();
201 Belos::ReturnType ret = solver->solve();
202 TEST_EQUALITY_CONST( ret, Belos::Converged );
213 ArrayRCP<Scalar> x_view =
x->get1dViewNonConst();
215 for (
size_t i=0; i<num_my_row; ++i) {
225 if (ST::magnitude(v) < tol)
228 TEST_FLOATING_EQUALITY(v,
val, tol);
240 typedef KokkosClassic::DefaultNode::DefaultNodeType
Node;
243 Tpetra_CrsMatrix, MatVec,
double,
int,
int,
Node )
248 Teuchos::GlobalMPISession mpiSession(&argc, &
argv);
251 int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc,
argv);
KokkosClassic::DefaultNode::DefaultNodeType Node
Node int main(int argc, char *argv[])
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix, MatVec, double, int, int, Node) TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(Tpetra_CrsMatrix
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix, MatVec, Scalar, LocalOrdinal, GlobalOrdinal, Node)
KokkosClassic::DefaultNode::DefaultNodeType Node
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y