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_2_DECL(
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), Tpetra::StaticProfile));
47  Array<GlobalOrdinal> columnIndices(2);
48  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
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  ArrayRCP<Scalar> y_view = y->get1dViewNonConst();
100  ArrayRCP<Scalar> x_view = y->get1dViewNonConst();
101  for (size_t i=0; i<num_my_row; ++i) {
102  //const GlobalOrdinal row = myGIDs[i];
103  Scalar val = 2.0;
104 
105  // if (row != nrow-1) {
106  // val += 2.0;
107  // }
108  TEST_EQUALITY( y_view[i], val );
109  }
110 }
111 
112 #if defined(HAVE_STOKHOS_BELOS)
113 
114 //
115 // Test Belos GMRES solve for a simple banded upper-triangular matrix
116 //
117 TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(
118  Tpetra_CrsMatrix, BelosGMRES, Scalar, Node )
119 {
120  using Teuchos::RCP;
121  using Teuchos::rcp;
122  using Teuchos::ArrayView;
123  using Teuchos::Array;
124  using Teuchos::ArrayRCP;
125  using Teuchos::ParameterList;
126  using LocalOrdinal = Tpetra::Map<>::local_ordinal_type;
127  using GlobalOrdinal = Tpetra::Map<>::global_ordinal_type;
128 
129  typedef Teuchos::Comm<int> Tpetra_Comm;
130  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
131  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
132  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
133  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
134 
135  // Build banded matrix
136  GlobalOrdinal nrow = 10;
137  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
138  RCP<const Tpetra_Map> map =
139  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
140  nrow, comm);
141  RCP<Tpetra_CrsGraph> graph =
142  rcp(new Tpetra_CrsGraph(map, size_t(2), Tpetra::StaticProfile));
143  Array<GlobalOrdinal> columnIndices(2);
144  ArrayView<const GlobalOrdinal> myGIDs = map->getNodeElementList();
145  const size_t num_my_row = myGIDs.size();
146  for (size_t i=0; i<num_my_row; ++i) {
147  const GlobalOrdinal row = myGIDs[i];
148  columnIndices[0] = row;
149  size_t ncol = 1;
150  // if (row != nrow-1) {
151  // columnIndices[1] = row+1;
152  // ncol = 2;
153  // }
154  graph->insertGlobalIndices(row, columnIndices(0,ncol));
155  }
156  graph->fillComplete();
157  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
158 
159  // Set values in matrix
160  Array<Scalar> vals(2);
161  for (size_t i=0; i<num_my_row; ++i) {
162  const GlobalOrdinal row = myGIDs[i];
163  columnIndices[0] = row;
164  vals[0] = Scalar(2.0);
165  size_t ncol = 1;
166 
167  // if (row != nrow-1) {
168  // columnIndices[1] = row+1;
169  // vals[1] = Scalar(2.0);
170  // ncol = 2;
171  // }
172  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
173  }
174  matrix->fillComplete();
175 
176  // Fill RHS vector
177  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
178  ArrayRCP<Scalar> b_view = b->get1dViewNonConst();
179  for (size_t i=0; i<num_my_row; ++i) {
180  b_view[i] = Scalar(1.0);
181  }
182 
183  // Solve
184  typedef Teuchos::ScalarTraits<Scalar> ST;
185  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
186  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
187  typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
188  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
189  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
190  RCP<ParameterList> belosParams = rcp(new ParameterList);
191  typename ST::magnitudeType tol = 1e-9;
192  belosParams->set("Flexible Gmres", false);
193  belosParams->set("Num Blocks", 100);
194  belosParams->set("Convergence Tolerance", Scalar(tol));
195  belosParams->set("Maximum Iterations", 100);
196  belosParams->set("Verbosity", 33);
197  belosParams->set("Output Style", 1);
198  belosParams->set("Output Frequency", 1);
199  belosParams->set("Output Stream", out.getOStream());
200  RCP<Belos::SolverManager<Scalar,MV,OP> > solver =
201  rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem, belosParams));
202  problem->setProblem();
203  Belos::ReturnType ret = solver->solve();
204  TEST_EQUALITY_CONST( ret, Belos::Converged );
205 
206  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
207  // Teuchos::VERB_EXTREME);
208 
209  // Check -- Correct answer is:
210  // [ 0, 0, ..., 0 ]
211  // [ 1, 1/2, ..., 1/VectorSize ]
212  // [ 0, 0, ..., 0 ]
213  // [ 1, 1/2, ..., 1/VectorSize ]
214  // ....
215  ArrayRCP<Scalar> x_view = x->get1dViewNonConst();
216  Scalar val = Scalar(0.5);
217  for (size_t i=0; i<num_my_row; ++i) {
218  // const GlobalOrdinal row = myGIDs[i];
219  // if (row % 2) {
220  // val = Scalar(0.5);
221  // }
222  // else
223  // val = Scalar(0.0);
224 
225  // Set small values to zero
226  Scalar v = x_view[i];
227  if (ST::magnitude(v) < tol)
228  v = Scalar(0.0);
229 
230  TEST_FLOATING_EQUALITY(v, val, tol);
231  }
232 }
233 
234 #else
235 
236 TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(
237  Tpetra_CrsMatrix, BelosGMRES, Scalar, Node )
238 {}
239 
240 #endif
241 
242 using Node = Tpetra::Map<>::node_type;
243 
244 TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(
245  Tpetra_CrsMatrix, MatVec, double, Node )
246 TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(
247  Tpetra_CrsMatrix, BelosGMRES, double, Node )
248 
249 int main( int argc, char* argv[] ) {
250  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
251 
252  // Run tests
253  int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
254 
255  return ret;
256 }
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:286
Teuchos_UnitTestHelpers.hpp
TEST_EQUALITY_CONST
#define TEST_EQUALITY_CONST(v1, v2)
node_type
Kokkos::Serial node_type
Definition: Stokhos_SacadoTraitsUnitTest.cpp:59
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_UNIT_TEST_TEMPLATE_2_DECL
TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(Kokkos_SG_SpMv, CrsProductTensorCijk, Scalar, Device)
Definition: Stokhos_KokkosArrayKernelsUnitTest_Host.hpp:43
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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
TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT
TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(Tpetra_CrsMatrix, MatVec, double, Node) TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(Tpetra_CrsMatrix
val
expr val()
Teuchos_GlobalMPISession.hpp
BelosGMRES
BelosGMRES
Definition: tpetra_mat_vec.cpp:247
TEST_EQUALITY
#define TEST_EQUALITY(v1, v2)
Stokhos_ConfigDefs.h
Generated by   doxygen 1.8.5