Stokhos Package Browser (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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 = Tpetra::getDefaultComm();
40  RCP<const Tpetra_Map> map =
41  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
42  nrow, comm);
43  RCP<Tpetra_CrsGraph> graph =
44  rcp(new Tpetra_CrsGraph(map, size_t(2), Tpetra::StaticProfile));
45  Array<GlobalOrdinal> columnIndices(2);
46  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
47  const size_t num_my_row = myGIDs.size();
48  for (size_t i=0; i<num_my_row; ++i) {
49  const GlobalOrdinal row = myGIDs[i];
50  columnIndices[0] = row;
51  size_t ncol = 1;
52 
53  // if (row != nrow-1) {
54  // columnIndices[1] = row+1;
55  // ncol = 2;
56  // }
57  graph->insertGlobalIndices(row, columnIndices(0,ncol));
58  }
59  graph->fillComplete();
60  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
61 
62  // Set values in matrix
63  Array<Scalar> vals(2);
64  for (size_t i=0; i<num_my_row; ++i) {
65  const GlobalOrdinal row = myGIDs[i];
66  columnIndices[0] = row;
67  vals[0] = Scalar(2.0);
68  size_t ncol = 1;
69 
70  // if (row != nrow-1) {
71  // columnIndices[1] = row+1;
72  // vals[1] = Scalar(2.0);
73  // ncol = 2;
74  // }
75  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
76  }
77  matrix->fillComplete();
78 
79  // Fill vector
80  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
81  x->putScalar(1.0);
82 
83  matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
84  Teuchos::VERB_EXTREME);
85 
86  x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
87  Teuchos::VERB_EXTREME);
88 
89  // Multiply
90  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
91  matrix->apply(*x, *y);
92 
93  y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
94  Teuchos::VERB_EXTREME);
95 
96  // Check
97  ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
98  ArrayRCP<Scalar> x_view = y->get1dViewNonConst();
99  for (size_t i=0; i<num_my_row; ++i) {
100  //const GlobalOrdinal row = myGIDs[i];
101  Scalar val = 2.0;
102 
103  // if (row != nrow-1) {
104  // val += 2.0;
105  // }
106  TEST_EQUALITY( y_view[i], val );
107  }
108 }
109 
110 #if defined(HAVE_STOKHOS_BELOS)
111 
112 //
113 // Test Belos GMRES solve for a simple banded upper-triangular matrix
114 //
115 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
116  Tpetra_CrsMatrix, BelosGMRES, Scalar, LocalOrdinal, GlobalOrdinal, Node )
117 {
118  using Teuchos::RCP;
119  using Teuchos::rcp;
120  using Teuchos::ArrayView;
121  using Teuchos::Array;
122  using Teuchos::ArrayRCP;
123  using Teuchos::ParameterList;
124 
125  typedef Teuchos::Comm<int> Tpetra_Comm;
126  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
127  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
128  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
129  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
130 
131  // Build banded matrix
132  GlobalOrdinal nrow = 10;
133  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
134  RCP<const Tpetra_Map> map =
135  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
136  nrow, comm);
137  RCP<Tpetra_CrsGraph> graph =
138  rcp(new Tpetra_CrsGraph(map, size_t(2), Tpetra::StaticProfile));
139  Array<GlobalOrdinal> columnIndices(2);
140  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
141  const size_t num_my_row = myGIDs.size();
142  for (size_t i=0; i<num_my_row; ++i) {
143  const GlobalOrdinal row = myGIDs[i];
144  columnIndices[0] = row;
145  size_t ncol = 1;
146  // if (row != nrow-1) {
147  // columnIndices[1] = row+1;
148  // ncol = 2;
149  // }
150  graph->insertGlobalIndices(row, columnIndices(0,ncol));
151  }
152  graph->fillComplete();
153  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
154 
155  // Set values in matrix
156  Array<Scalar> vals(2);
157  for (size_t i=0; i<num_my_row; ++i) {
158  const GlobalOrdinal row = myGIDs[i];
159  columnIndices[0] = row;
160  vals[0] = Scalar(2.0);
161  size_t ncol = 1;
162 
163  // if (row != nrow-1) {
164  // columnIndices[1] = row+1;
165  // vals[1] = Scalar(2.0);
166  // ncol = 2;
167  // }
168  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
169  }
170  matrix->fillComplete();
171 
172  // Fill RHS vector
173  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
174  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
175  for (size_t i=0; i<num_my_row; ++i) {
176  b_view[i] = Scalar(1.0);
177  }
178 
179  // Solve
180  typedef Teuchos::ScalarTraits<Scalar> ST;
181  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
182  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
183  typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
184  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
185  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
186  RCP<ParameterList> belosParams = rcp(new ParameterList);
187  typename ST::magnitudeType tol = 1e-9;
188  belosParams->set("Flexible Gmres", false);
189  belosParams->set("Num Blocks", 100);
190  belosParams->set("Convergence Tolerance", Scalar(tol));
191  belosParams->set("Maximum Iterations", 100);
192  belosParams->set("Verbosity", 33);
193  belosParams->set("Output Style", 1);
194  belosParams->set("Output Frequency", 1);
195  belosParams->set("Output Stream", out.getOStream());
196  RCP<Belos::SolverManager<Scalar,MV,OP> > solver =
197  rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem, belosParams));
198  problem->setProblem();
199  Belos::ReturnType ret = solver->solve();
200  TEST_EQUALITY_CONST( ret, Belos::Converged );
201 
202  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
203  // Teuchos::VERB_EXTREME);
204 
205  // Check -- Correct answer is:
206  // [ 0, 0, ..., 0 ]
207  // [ 1, 1/2, ..., 1/VectorSize ]
208  // [ 0, 0, ..., 0 ]
209  // [ 1, 1/2, ..., 1/VectorSize ]
210  // ....
211  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
212  Scalar val = Scalar(0.5);
213  for (size_t i=0; i<num_my_row; ++i) {
214  // const GlobalOrdinal row = myGIDs[i];
215  // if (row % 2) {
216  // val = Scalar(0.5);
217  // }
218  // else
219  // val = Scalar(0.0);
220 
221  // Set small values to zero
222  Scalar v = x_view[i];
223  if (ST::magnitude(v) < tol)
224  v = Scalar(0.0);
225 
226  TEST_FLOATING_EQUALITY(v, val, tol);
227  }
228 }
229 
230 #else
231 
232 TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(
233  Tpetra_CrsMatrix, BelosGMRES, Scalar, LocalOrdinal, GlobalOrdinal, Node )
234 {}
235 
236 #endif
237 
238 typedef KokkosClassic::DefaultNode::DefaultNodeType Node;
239 
240 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(
241  Tpetra_CrsMatrix, MatVec, double, int, int, Node )
242 TEUCHOS_UNIT_TEST_TEMPLATE_4_INSTANT(
243  Tpetra_CrsMatrix, BelosGMRES, double, int, int, Node )
244 
245 int main( int argc, char* argv[] ) {
246  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
247 
248  // Run tests
249  int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
250 
251  return ret;
252 }
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:286
Teuchos_UnitTestHelpers.hpp
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL
TEUCHOS_UNIT_TEST_TEMPLATE_4_DECL(Tpetra_CrsMatrix_MP, VectorAdd, Storage, LocalOrdinal, GlobalOrdinal, Node)
Definition: Stokhos_TpetraCrsMatrixMPVectorUnitTest.hpp:124
TEST_EQUALITY_CONST
#define TEST_EQUALITY_CONST(v1, v2)
Teuchos_UnitTestRepository.hpp
Scalar
pce_type Scalar
Definition: linear2d_diffusion_scalar_types.hpp:49
Teuchos::UnitTestRepository::runUnitTestsFromMain
static int runUnitTestsFromMain(int argc, char *argv[])
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
LocalOrdinal
int LocalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:50
Node
KokkosClassic::DefaultNode::DefaultNodeType Node
Definition: experimental/linear2d_diffusion_pce.cpp:79
Teuchos_UnitTestHarness.hpp
GlobalOrdinal
int GlobalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:51
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
main
int main(int argc, char **argv)
Definition: cijk_ltb_partition.cpp:57
val
expr val()
Teuchos_GlobalMPISession.hpp
BelosGMRES
BelosGMRES
Definition: tpetra_mat_vec.cpp:243
TEST_EQUALITY
#define TEST_EQUALITY(v1, v2)
Stokhos_ConfigDefs.h
Generated by   doxygen 1.8.5