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 // @HEADER
2 // *****************************************************************************
3 // Stokhos Package
4 //
5 // Copyright 2009 NTESS and the Stokhos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 // Testing utilties
11 #include "Teuchos_UnitTestHarness.hpp"
12 #include "Teuchos_UnitTestRepository.hpp"
13 #include "Teuchos_UnitTestHelpers.hpp"
14 #include "Teuchos_GlobalMPISession.hpp"
15 
16 // Tpetra
17 #include "Tpetra_Core.hpp"
18 #include "Tpetra_Map.hpp"
19 #include "Tpetra_Vector.hpp"
20 #include "Tpetra_CrsGraph.hpp"
21 #include "Tpetra_CrsMatrix.hpp"
22 
23 // Belos solver
24 #include "Stokhos_ConfigDefs.h"
25 #ifdef HAVE_STOKHOS_BELOS
26 #include "BelosTpetraAdapter.hpp"
27 #include "BelosLinearProblem.hpp"
28 #include "BelosPseudoBlockGmresSolMgr.hpp"
29 #endif
30 
31 TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(
32  Tpetra_CrsMatrix, MatVec, Scalar, Node )
33 {
34  using Teuchos::RCP;
35  using Teuchos::rcp;
36  using Teuchos::ArrayView;
37  using Teuchos::Array;
38  using Teuchos::ArrayRCP;
39  using LocalOrdinal = Tpetra::Map<>::local_ordinal_type;
40  using GlobalOrdinal = Tpetra::Map<>::global_ordinal_type;
41 
42  typedef Teuchos::Comm<int> Tpetra_Comm;
43  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
44  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
45  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
46  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
47 
48  // Build banded matrix
49  GlobalOrdinal nrow = 10;
50  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
51  RCP<const Tpetra_Map> map =
52  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
53  nrow, comm);
54  RCP<Tpetra_CrsGraph> graph =
55  rcp(new Tpetra_CrsGraph(map, size_t(2)));
56  Array<GlobalOrdinal> columnIndices(2);
57  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
58  const size_t num_my_row = myGIDs.size();
59  for (size_t i=0; i<num_my_row; ++i) {
60  const GlobalOrdinal row = myGIDs[i];
61  columnIndices[0] = row;
62  size_t ncol = 1;
63 
64  // if (row != nrow-1) {
65  // columnIndices[1] = row+1;
66  // ncol = 2;
67  // }
68  graph->insertGlobalIndices(row, columnIndices(0,ncol));
69  }
70  graph->fillComplete();
71  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
72 
73  // Set values in matrix
74  Array<Scalar> vals(2);
75  for (size_t i=0; i<num_my_row; ++i) {
76  const GlobalOrdinal row = myGIDs[i];
77  columnIndices[0] = row;
78  vals[0] = Scalar(2.0);
79  size_t ncol = 1;
80 
81  // if (row != nrow-1) {
82  // columnIndices[1] = row+1;
83  // vals[1] = Scalar(2.0);
84  // ncol = 2;
85  // }
86  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
87  }
88  matrix->fillComplete();
89 
90  // Fill vector
91  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
92  x->putScalar(1.0);
93 
94  matrix->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
95  Teuchos::VERB_EXTREME);
96 
97  x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
98  Teuchos::VERB_EXTREME);
99 
100  // Multiply
101  RCP<Tpetra_Vector> y = Tpetra::createVector<Scalar>(map);
102  matrix->apply(*x, *y);
103 
104  y->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
105  Teuchos::VERB_EXTREME);
106 
107  // Check
108  auto y_view = y->getLocalViewHost(Tpetra::Access::ReadOnly);
109  for (size_t i=0; i<num_my_row; ++i) {
110  //const GlobalOrdinal row = myGIDs[i];
111  Scalar val = 2.0;
112 
113  // if (row != nrow-1) {
114  // val += 2.0;
115  // }
116  TEST_EQUALITY( y_view(i,0), val );
117  }
118 }
119 
120 #if defined(HAVE_STOKHOS_BELOS)
121 
122 //
123 // Test Belos GMRES solve for a simple banded upper-triangular matrix
124 //
125 TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(
126  Tpetra_CrsMatrix, BelosGMRES, Scalar, Node )
127 {
128  using Teuchos::RCP;
129  using Teuchos::rcp;
130  using Teuchos::ArrayView;
131  using Teuchos::Array;
132  using Teuchos::ArrayRCP;
133  using Teuchos::ParameterList;
134  using LocalOrdinal = Tpetra::Map<>::local_ordinal_type;
135  using GlobalOrdinal = Tpetra::Map<>::global_ordinal_type;
136 
137  typedef Teuchos::Comm<int> Tpetra_Comm;
138  typedef Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> Tpetra_Map;
139  typedef Tpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_Vector;
140  typedef Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsMatrix;
141  typedef Tpetra::CrsGraph<LocalOrdinal,GlobalOrdinal,Node> Tpetra_CrsGraph;
142 
143  // Build banded matrix
144  GlobalOrdinal nrow = 10;
145  RCP<const Tpetra_Comm> comm = Tpetra::getDefaultComm();
146  RCP<const Tpetra_Map> map =
147  Tpetra::createUniformContigMapWithNode<LocalOrdinal,GlobalOrdinal,Node>(
148  nrow, comm);
149  RCP<Tpetra_CrsGraph> graph =
150  rcp(new Tpetra_CrsGraph(map, size_t(2)));
151  Array<GlobalOrdinal> columnIndices(2);
152  ArrayView<const GlobalOrdinal> myGIDs = map->getLocalElementList();
153  const size_t num_my_row = myGIDs.size();
154  for (size_t i=0; i<num_my_row; ++i) {
155  const GlobalOrdinal row = myGIDs[i];
156  columnIndices[0] = row;
157  size_t ncol = 1;
158  // if (row != nrow-1) {
159  // columnIndices[1] = row+1;
160  // ncol = 2;
161  // }
162  graph->insertGlobalIndices(row, columnIndices(0,ncol));
163  }
164  graph->fillComplete();
165  RCP<Tpetra_CrsMatrix> matrix = rcp(new Tpetra_CrsMatrix(graph));
166 
167  // Set values in matrix
168  Array<Scalar> vals(2);
169  for (size_t i=0; i<num_my_row; ++i) {
170  const GlobalOrdinal row = myGIDs[i];
171  columnIndices[0] = row;
172  vals[0] = Scalar(2.0);
173  size_t ncol = 1;
174 
175  // if (row != nrow-1) {
176  // columnIndices[1] = row+1;
177  // vals[1] = Scalar(2.0);
178  // ncol = 2;
179  // }
180  matrix->replaceGlobalValues(row, columnIndices(0,ncol), vals(0,ncol));
181  }
182  matrix->fillComplete();
183 
184  // Fill RHS vector
185  RCP<Tpetra_Vector> b = Tpetra::createVector<Scalar>(map);
186  {
187  auto b_view = b->getLocalViewHost(Tpetra::Access::OverwriteAll);
188  for (size_t i=0; i<num_my_row; ++i) {
189  b_view(i, 0) = Scalar(1.0);
190  }
191  }
192 
193  // Solve
194  typedef Teuchos::ScalarTraits<Scalar> ST;
195  typedef Tpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> MV;
196  typedef Tpetra::Operator<Scalar,LocalOrdinal,GlobalOrdinal,Node> OP;
197  typedef Belos::LinearProblem<Scalar,MV,OP> BLinProb;
198  RCP<Tpetra_Vector> x = Tpetra::createVector<Scalar>(map);
199  RCP< BLinProb > problem = rcp(new BLinProb(matrix, x, b));
200  RCP<ParameterList> belosParams = rcp(new ParameterList);
201  typename ST::magnitudeType tol = 1e-9;
202  belosParams->set("Flexible Gmres", false);
203  belosParams->set("Num Blocks", 100);
204  belosParams->set("Convergence Tolerance", Scalar(tol));
205  belosParams->set("Maximum Iterations", 100);
206  belosParams->set("Verbosity", 33);
207  belosParams->set("Output Style", 1);
208  belosParams->set("Output Frequency", 1);
209  belosParams->set("Output Stream", out.getOStream());
210  RCP<Belos::SolverManager<Scalar,MV,OP> > solver =
211  rcp(new Belos::PseudoBlockGmresSolMgr<Scalar,MV,OP>(problem, belosParams));
212  problem->setProblem();
213  Belos::ReturnType ret = solver->solve();
214  TEST_EQUALITY_CONST( ret, Belos::Converged );
215 
216  // x->describe(*(Teuchos::fancyOStream(rcp(&std::cout,false))),
217  // Teuchos::VERB_EXTREME);
218 
219  // Check -- Correct answer is:
220  // [ 0, 0, ..., 0 ]
221  // [ 1, 1/2, ..., 1/VectorSize ]
222  // [ 0, 0, ..., 0 ]
223  // [ 1, 1/2, ..., 1/VectorSize ]
224  // ....
225  auto x_view = x->getLocalViewHost(Tpetra::Access::ReadOnly);
226  Scalar val = Scalar(0.5);
227  for (size_t i=0; i<num_my_row; ++i) {
228  // const GlobalOrdinal row = myGIDs[i];
229  // if (row % 2) {
230  // val = Scalar(0.5);
231  // }
232  // else
233  // val = Scalar(0.0);
234 
235  // Set small values to zero
236  Scalar v = x_view(i,0);
237  if (ST::magnitude(v) < tol)
238  v = Scalar(0.0);
239 
240  TEST_FLOATING_EQUALITY(v, val, tol);
241  }
242 }
243 
244 #else
245 
246 TEUCHOS_UNIT_TEST_TEMPLATE_2_DECL(
247  Tpetra_CrsMatrix, BelosGMRES, Scalar, Node )
248 {}
249 
250 #endif
251 
252 using Node = Tpetra::Map<>::node_type;
253 
254 TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(
255  Tpetra_CrsMatrix, MatVec, double, Node )
256 TEUCHOS_UNIT_TEST_TEMPLATE_2_INSTANT(
257  Tpetra_CrsMatrix, BelosGMRES, double, Node )
258 
259 int main( int argc, char* argv[] ) {
260  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
261 
262  // Run tests
263  int ret = Teuchos::UnitTestRepository::runUnitTestsFromMain(argc, argv);
264 
265  return ret;
266 }
argv
char * argv[]
Definition: Stokhos_HouseTriDiagUnitTest.cpp:252
Teuchos_UnitTestHelpers.hpp
TEST_EQUALITY_CONST
#define TEST_EQUALITY_CONST(v1, v2)
node_type
Kokkos::Serial node_type
Definition: Stokhos_SacadoTraitsUnitTest.cpp:29
Teuchos_UnitTestRepository.hpp
Scalar
pce_type Scalar
Definition: linear2d_diffusion_scalar_types.hpp:17
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:11
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
LocalOrdinal
int LocalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:18
Teuchos_UnitTestHarness.hpp
GlobalOrdinal
int GlobalOrdinal
Definition: linear2d_diffusion_scalar_types.hpp:19
TEST_FLOATING_EQUALITY
#define TEST_FLOATING_EQUALITY(v1, v2, tol)
main
int main(int argc, char **argv)
Definition: cijk_ltb_partition.cpp:25
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
Node
Tpetra::KokkosClassic::DefaultNode::DefaultNodeType Node
Definition: experimental/linear2d_diffusion_pce.cpp:47
BelosGMRES
BelosGMRES
Definition: tpetra_mat_vec.cpp:257
TEST_EQUALITY
#define TEST_EQUALITY(v1, v2)
Stokhos_ConfigDefs.h
Generated on Fri Nov 22 2024 09:18:16 for Stokhos Package Browser (Single Doxygen Collection) by   doxygen 1.8.5