Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
tpetra_mat_vec.cpp
Go to the documentation of this file.
1// Testing utilties
2#include "Teuchos_UnitTestHarness.hpp"
3#include "Teuchos_UnitTestRepository.hpp"
4#include "Teuchos_UnitTestHelpers.hpp"
5#include "Teuchos_GlobalMPISession.hpp"
6
7// Tpetra
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"
13
14// Belos solver
15#include "Stokhos_ConfigDefs.h"
16#ifdef HAVE_STOKHOS_BELOS
17#include "BelosTpetraAdapter.hpp"
18#include "BelosLinearProblem.hpp"
19#include "BelosPseudoBlockGmresSolMgr.hpp"
20#endif
21
23 Tpetra_CrsMatrix, MatVec, Scalar, Node )
24{
25 using Teuchos::RCP;
26 using Teuchos::rcp;
27 using Teuchos::ArrayView;
28 using Teuchos::Array;
29 using Teuchos::ArrayRCP;
30 using LocalOrdinal = Tpetra::Map<>::local_ordinal_type;
31 using GlobalOrdinal = Tpetra::Map<>::global_ordinal_type;
32
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;
38
39 // Build banded matrix
40 GlobalOrdinal nrow = 10;
41 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
42 RCP<const Tpetra_Map> map =
43 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
44 nrow, comm);
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) {
51 const GlobalOrdinal row = myGIDs[i];
52 columnIndices[0] = row;
53 size_t ncol = 1;
54
55 // if (row != nrow-1) {
56 // columnIndices[1] = row+1;
57 // ncol = 2;
58 // }
59 graph->insertGlobalIndices(row, columnIndices(0,ncol));
60 }
61 graph->fillComplete();
62 RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
63
64 // Set values in matrix
65 Array<Scalar> vals(2);
66 for (size_t i=0; i<num_my_row; ++i) {
67 const GlobalOrdinal row = myGIDs[i];
68 columnIndices[0] = row;
69 vals[0] = Scalar(2.0);
70 size_t ncol = 1;
71
72 // if (row != nrow-1) {
73 // columnIndices[1] = row+1;
74 // vals[1] = Scalar(2.0);
75 // ncol = 2;
76 // }
77 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
78 }
79 matrix->fillComplete();
80
81 // Fill vector
82 RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
83 x->putScalar(1.0);
84
85 matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
86 Teuchos::VERB_EXTREME);
87
88 x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
89 Teuchos::VERB_EXTREME);
90
91 // Multiply
92 RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
93 matrix->apply(*x, *y);
94
95 y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
96 Teuchos::VERB_EXTREME);
97
98 // Check
99 auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
100 for (size_t i=0; i<num_my_row; ++i) {
101 //const GlobalOrdinal row = myGIDs[i];
102 Scalar val = 2.0;
103
104 // if (row != nrow-1) {
105 // val += 2.0;
106 // }
107 TEST_EQUALITY( y_view(i,0), val );
108 }
109}
110
111#if defined(HAVE_STOKHOS_BELOS)
112
113//
114// Test Belos GMRES solve for a simple banded upper-triangular matrix
115//
117 Tpetra_CrsMatrix, BelosGMRES, Scalar, Node )
118{
119 using Teuchos::RCP;
120 using Teuchos::rcp;
121 using Teuchos::ArrayView;
122 using Teuchos::Array;
123 using Teuchos::ArrayRCP;
124 using Teuchos::ParameterList;
125 using LocalOrdinal = Tpetra::Map<>::local_ordinal_type;
126 using GlobalOrdinal = Tpetra::Map<>::global_ordinal_type;
127
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;
133
134 // Build banded matrix
135 GlobalOrdinal nrow = 10;
136 RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
137 RCP<const Tpetra_Map> map =
138 Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
139 nrow, comm);
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) {
146 const GlobalOrdinal row = myGIDs[i];
147 columnIndices[0] = row;
148 size_t ncol = 1;
149 // if (row != nrow-1) {
150 // columnIndices[1] = row+1;
151 // ncol = 2;
152 // }
153 graph->insertGlobalIndices(row, columnIndices(0,ncol));
154 }
155 graph->fillComplete();
156 RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
157
158 // Set values in matrix
159 Array<Scalar> vals(2);
160 for (size_t i=0; i<num_my_row; ++i) {
161 const GlobalOrdinal row = myGIDs[i];
162 columnIndices[0] = row;
163 vals[0] = Scalar(2.0);
164 size_t ncol = 1;
165
166 // if (row != nrow-1) {
167 // columnIndices[1] = row+1;
168 // vals[1] = Scalar(2.0);
169 // ncol = 2;
170 // }
171 matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
172 }
173 matrix->fillComplete();
174
175 // Fill RHS vector
176 RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
177 {
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);
181 }
182 }
183
184 // Solve
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 );
206
207 // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
208 // Teuchos::VERB_EXTREME);
209
210 // Check -- Correct answer is:
211 // [ 0, 0, ..., 0 ]
212 // [ 1, 1/2, ..., 1/VectorSize ]
213 // [ 0, 0, ..., 0 ]
214 // [ 1, 1/2, ..., 1/VectorSize ]
215 // ....
216 auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
217 Scalar val = Scalar(0.5);
218 for (size_t i=0; i<num_my_row; ++i) {
219 // const GlobalOrdinal row = myGIDs[i];
220 // if (row % 2) {
221 // val = Scalar(0.5);
222 // }
223 // else
224 // val = Scalar(0.0);
225
226 // Set small values to zero
227 Scalar v = x_view(i,0);
228 if (ST::magnitude(v) < tol)
229 v = Scalar(0.0);
230
231 TEST_FLOATING_EQUALITY(v, val, tol);
232 }
233}
234
235#else
236
238 Tpetra_CrsMatrix, BelosGMRES, Scalar, Node )
239{}
240
241#endif
242
243using Node = Tpetra::Map<>::node_type;
244
246 Tpetra_CrsMatrix, MatVec, double, Node )
248 Tpetra_CrsMatrix, BelosGMRES, double, Node )
249
250int main( int argc, char* argv[] ) {
251 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
252
253 // Run tests
254 int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
255
256 return ret;
257}
expr val()
KokkosClassic::DefaultNode::DefaultNodeType Node
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:265
TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(Tpetra_CrsMatrix, MatVec, double, Node) TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(Tpetra_CrsMatrix
BelosGMRES
Node int main(int argc, char *argv[])
TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(Tpetra_CrsMatrix, MatVec, Scalar, Node)