Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
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_DefaultPlatform.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 
22 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
23  Tpetra_CrsMatrix, MatVec, Scalar, LocalOrdinal, GlobalOrdinal, Node )
24 {
25  using Teuchos::RCP;
26  using Teuchos::rcp;
27  using Teuchos::ArrayView;
28  using Teuchos::Array;
29  using Teuchos::ArrayRCP;
30 
31  typedef Teuchos::Comm<int> Tpetra_Comm;
32  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
33  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
34  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
35  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
36 
37  // Build banded matrix
38  GlobalOrdinal nrow = 10;
39  RCP<const Tpetra_Comm> comm =
40  Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
41  RCP<const Tpetra_Map> map =
42  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
43  nrow, comm);
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) {
50  const GlobalOrdinal row = myGIDs[i];
51  columnIndices[0] = row;
52  size_t ncol = 1;
53 
54  // if (row != nrow-1) {
55  // columnIndices[1] = row+1;
56  // ncol = 2;
57  // }
58  graph->insertGlobalIndices(row, columnIndices(0,ncol));
59  }
60  graph->fillComplete();
61  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
62 
63  // Set values in matrix
64  Array<Scalar> vals(2);
65  for (size_t i=0; i<num_my_row; ++i) {
66  const GlobalOrdinal row = myGIDs[i];
67  columnIndices[0] = row;
68  vals[0] = Scalar(2.0);
69  size_t ncol = 1;
70 
71  // if (row != nrow-1) {
72  // columnIndices[1] = row+1;
73  // vals[1] = Scalar(2.0);
74  // ncol = 2;
75  // }
76  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
77  }
78  matrix->fillComplete();
79 
80  // Fill vector
81  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
82  x->putScalar(1.0);
83 
84  matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
85  Teuchos::VERB_EXTREME);
86 
87  x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
88  Teuchos::VERB_EXTREME);
89 
90  // Multiply
91  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
92  matrix->apply(*x, *y);
93 
94  y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
95  Teuchos::VERB_EXTREME);
96 
97  // Check
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) {
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], 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 //
116 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
117  Tpetra_CrsMatrix, BelosGMRES, Scalar, LocalOrdinal, GlobalOrdinal, 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 
126  typedef Teuchos::Comm<int> Tpetra_Comm;
127  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
128  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
129  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
130  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
131 
132  // Build banded matrix
133  GlobalOrdinal nrow = 10;
134  RCP<const Tpetra_Comm> comm =
135  Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
136  RCP<const Tpetra_Map> map =
137  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
138  nrow, comm);
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) {
145  const GlobalOrdinal row = myGIDs[i];
146  columnIndices[0] = row;
147  size_t ncol = 1;
148  // if (row != nrow-1) {
149  // columnIndices[1] = row+1;
150  // ncol = 2;
151  // }
152  graph->insertGlobalIndices(row, columnIndices(0,ncol));
153  }
154  graph->fillComplete();
155  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
156 
157  // Set values in matrix
158  Array<Scalar> vals(2);
159  for (size_t i=0; i<num_my_row; ++i) {
160  const GlobalOrdinal row = myGIDs[i];
161  columnIndices[0] = row;
162  vals[0] = Scalar(2.0);
163  size_t ncol = 1;
164 
165  // if (row != nrow-1) {
166  // columnIndices[1] = row+1;
167  // vals[1] = Scalar(2.0);
168  // ncol = 2;
169  // }
170  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
171  }
172  matrix->fillComplete();
173 
174  // Fill RHS vector
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) {
178  b_view[i] = Scalar(1.0);
179  }
180 
181  // Solve
182  typedef Teuchos::ScalarTraits<Scalar> ST;
183  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
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 );
203 
204  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
205  // Teuchos::VERB_EXTREME);
206 
207  // Check -- Correct answer is:
208  // [ 0, 0, ..., 0 ]
209  // [ 1, 1/2, ..., 1/VectorSize ]
210  // [ 0, 0, ..., 0 ]
211  // [ 1, 1/2, ..., 1/VectorSize ]
212  // ....
213  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
214  Scalar val = Scalar(0.5);
215  for (size_t i=0; i<num_my_row; ++i) {
216  // const GlobalOrdinal row = myGIDs[i];
217  // if (row % 2) {
218  // val = Scalar(0.5);
219  // }
220  // else
221  // val = Scalar(0.0);
222 
223  // Set small values to zero
224  Scalar v = x_view[i];
225  if (ST::magnitude(v) < tol)
226  v = Scalar(0.0);
227 
228  TEST_FLOATING_EQUALITY(v, val, tol);
229  }
230 }
231 
232 #else
233 
234 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
235  Tpetra_CrsMatrix, BelosGMRES, Scalar, LocalOrdinal, GlobalOrdinal, Node )
236 {}
237 
238 #endif
239 
240 typedef KokkosClassic::DefaultNode::DefaultNodeType Node;
241 
242 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(
243  Tpetra_CrsMatrix, MatVec, double, int, int, Node )
244 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(
245  Tpetra_CrsMatrix, BelosGMRES, double, int, int, Node )
246 
247 int main( int argc, char* argv[] ) {
248  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
249 
250  // Run tests
251  int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
252 
253  return ret;
254 }
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:286
Teuchos_UnitTestHelpers.hpp
val
expr val()
TEST_EQUALITY_CONST
#define TEST_EQUALITY_CONST(v1, v2)
Node
KokkosClassic::DefaultNode::DefaultNodeType Node
Definition: tpetra_mat_vec.cpp:240
main
Node int main(int argc, char *argv[])
Definition: tpetra_mat_vec.cpp:247
Teuchos_UnitTestRepository.hpp
Teuchos::UnitTestRepository::runUnitTestsFromMain
static int runUnitTestsFromMain(int argc, char *argv[])
cusp::detail::device::x
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType * x
Definition: csr_vector.h:260
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT
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
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix, MatVec, Scalar, LocalOrdinal, GlobalOrdinal, Node)
Definition: tpetra_mat_vec.cpp:22
Node
KokkosClassic::DefaultNode::DefaultNodeType Node
Definition: experimental/linear2d_diffusion_pce.cpp:79
Teuchos_UnitTestHarness.hpp
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
Teuchos_GlobalMPISession.hpp
int
int
Definition: tpetra_mat_vec.cpp:245
Scalar
pce_type Scalar
Definition: ifpack2_pce_instantiations.cpp:46
LocalOrdinal
int LocalOrdinal
Definition: ifpack2_pce_instantiations.cpp:47
TEST_EQUALITY
TEST_EQUALITY(rcp_dynamic_cast< const EnhancedNumberValidator< double > >(castedDep1->getValuesAndValidators().find("val1") ->second, true) ->getMax(), double1Vali->getMax())
BelosGMRES
BelosGMRES
Definition: tpetra_mat_vec.cpp:245
Stokhos_ConfigDefs.h
cusp::detail::device::y
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267
GlobalOrdinal
int GlobalOrdinal
Definition: ifpack2_pce_instantiations.cpp:48
Generated by   doxygen 1.8.14