Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
TpetraCrsColorer.cpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Zoltan2: A package of combinatorial algorithms for scientific computing
4 //
5 // Copyright 2012 NTESS and the Zoltan2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #include "Tpetra_Core.hpp"
11 #include "Kokkos_Random.hpp"
12 #include "Zoltan2_TestHelpers.hpp"
14 
15 // Class to test the Colorer utility
16 class ColorerTest {
17 public:
18  using map_t = Tpetra::Map<>;
19  using gno_t = typename map_t::global_ordinal_type;
20  using graph_t = Tpetra::CrsGraph<>;
21  using matrix_t = Tpetra::CrsMatrix<zscalar_t>;
22  using multivector_t = Tpetra::MultiVector<zscalar_t>;
23  using execution_space_t = typename matrix_t::device_type::execution_space;
24 
26  // Construct the test:
27  // Read or generate a matrix (JBlock) with default range and domain maps
28  // Construct identical matrix (JCyclic) with cyclic range and domain maps
29 
30  ColorerTest(Teuchos::RCP<const Teuchos::Comm<int> > &comm,
31  int narg, char**arg)
32  : symmetric(false)
33  {
34  int me = comm->getRank();
35  int np = comm->getSize();
36 
37  // Process command line arguments
38  bool distributeInput = true;
39  size_t xdim = 10, ydim = 11, zdim = 12;
40 
41  Teuchos::CommandLineProcessor cmdp(false, false);
42  cmdp.setOption("file", &matrixFileName,
43  "Name of the Matrix Market file to use");
44  cmdp.setOption("xdim", &xdim,
45  "Number of nodes in x-direction for generated matrix");
46  cmdp.setOption("ydim", &ydim,
47  "Number of nodes in y-direction for generated matrix");
48  cmdp.setOption("zdim", &zdim,
49  "Number of nodes in z-direction for generated matrix");
50  cmdp.setOption("distribute", "no-distribute", &distributeInput,
51  "Should Zoltan2 distribute the matrix as it is read?");
52  cmdp.setOption("symmetric", "non-symmetric", &symmetric,
53  "Is the matrix symmetric?");
54  cmdp.parse(narg, arg);
55 
56  // Get and store a matrix
57  if (matrixFileName != "") {
58  // Read from a file
59  UserInputForTests uinput(".", matrixFileName, comm, true, distributeInput);
60  JBlock = uinput.getUITpetraCrsMatrix();
61  }
62  else {
63  // Generate a matrix
64  UserInputForTests uinput(xdim, ydim, zdim, string("Laplace3D"), comm,
65  true, distributeInput);
66  JBlock = uinput.getUITpetraCrsMatrix();
67  }
68 
69  // Build same matrix with cyclic domain and range maps
70  // To make range and domain maps differ for square matrices,
71  // keep some processors empty in the cyclic maps
72 
73  size_t nIndices = std::max(JBlock->getGlobalNumCols(),
74  JBlock->getGlobalNumRows());
75  Teuchos::Array<gno_t> indices(nIndices);
76 
77  Teuchos::RCP<const map_t> vMapCyclic =
78  getCyclicMap(JBlock->getGlobalNumCols(), indices, np-1, comm);
79  Teuchos::RCP<const map_t> wMapCyclic =
80  getCyclicMap(JBlock->getGlobalNumRows(), indices, np-2, comm);
81 
82  // Fill JBlock with random numbers for a better test.
83  JBlock->resumeFill();
84 
85 #if KOKKOS_VERSION >= 40799
86  using IST = typename KokkosKernels::ArithTraits<zscalar_t>::val_type;
87 #else
88  using IST = typename Kokkos::ArithTraits<zscalar_t>::val_type;
89 #endif
90  using pool_type =
91  Kokkos::Random_XorShift64_Pool<execution_space_t>;
92  pool_type rand_pool(static_cast<uint64_t>(me));
93 
94  Kokkos::fill_random(JBlock->getLocalMatrixDevice().values, rand_pool,
95  static_cast<IST>(1.), static_cast<IST>(9999.));
96  JBlock->fillComplete();
97 
98 
99  // Make JCyclic: same matrix with different Domain and Range maps
100  RCP<const graph_t> block_graph = JBlock->getCrsGraph();
101  RCP<graph_t> cyclic_graph = rcp(new graph_t(*block_graph));
102  cyclic_graph->resumeFill();
103  cyclic_graph->fillComplete(vMapCyclic, wMapCyclic);
104  JCyclic = rcp(new matrix_t(cyclic_graph));
105  JCyclic->resumeFill();
106  TEUCHOS_ASSERT(block_graph->getLocalNumRows() == cyclic_graph->getLocalNumRows());
107  {
108  auto val_s = JBlock->getLocalMatrixHost().values;
109  auto val_d = JCyclic->getLocalMatrixHost().values;
110  TEUCHOS_ASSERT(val_s.extent(0) == val_d.extent(0));
111  Kokkos::deep_copy(val_d, val_s);
112  }
113  JCyclic->fillComplete();
114  }
115 
117  bool run(const char *testname, Teuchos::ParameterList &params) {
118 
119  bool ok = true;
120 
121  params.set("symmetric", symmetric);
122  params.set("library", "zoltan");
123 
124  // test with default maps
125  ok = buildAndCheckSeedMatrix(testname, params, true);
126 
127  // test with cyclic maps
128  ok &= buildAndCheckSeedMatrix(testname, params, false);
129  // if (matrixFileName != "west0067") {
130 
131  params.set("library", "zoltan2");
132  // test with default maps
133  ok = buildAndCheckSeedMatrix(testname, params, true);
134 
135  // test with cyclic maps
136  ok &= buildAndCheckSeedMatrix(testname, params, false);
137  // }
138 
139  return ok;
140  }
141 
144  const char *testname,
145  Teuchos::ParameterList &params,
146  const bool useBlock
147  )
148  {
149  int ierr = 0;
150 
151  // Pick matrix depending on useBlock flag
152  Teuchos::RCP<matrix_t> J = (useBlock ? JBlock : JCyclic);
153  int me = J->getRowMap()->getComm()->getRank();
154 
155  std::cout << params.get("library", "zoltan2") << " Running " << testname << " with "
156  << (useBlock ? "Block maps" : "Cyclic maps")
157  << std::endl;
158 
159  // Create a colorer
161  colorer.computeColoring(params);
162 
163  // Check coloring
164  if (!colorer.checkColoring()) {
165  std::cout << testname << " with "
166  << (useBlock ? "Block maps" : "Cyclic maps")
167  << " FAILED: invalid coloring returned"
168  << std::endl;
169  return false;
170  }
171 
172  // Compute seed matrix V -- the application wants this matrix
173  // Dense matrix of 0/1 indicating the compression via coloring
174 
175  const int numColors = colorer.getNumColors();
176 
177  // Compute the seed matrix; applications want this seed matrix
178 
179  multivector_t V(J->getDomainMap(), numColors);
180  colorer.computeSeedMatrix(V);
181 
182  // To test the result...
183  // Compute the compressed matrix W
184  multivector_t W(J->getRangeMap(), numColors);
185 
186  J->apply(V, W);
187 
188  // Reconstruct matrix from compression vector
189  Teuchos::RCP<matrix_t> Jp = rcp(new matrix_t(*J, Teuchos::Copy));
190  Jp->setAllToScalar(static_cast<zscalar_t>(-1.));
191 
192  colorer.reconstructMatrix(W, *Jp);
193 
194  // Check that values of J = values of Jp
195  auto J_local_matrix = J->getLocalMatrixDevice();
196  auto Jp_local_matrix = Jp->getLocalMatrixDevice();
197  const size_t num_local_nz = J->getLocalNumEntries();
198 
199  Kokkos::parallel_reduce(
200  "TpetraCrsColorer::testReconstructedMatrix()",
201  Kokkos::RangePolicy<execution_space_t>(0, num_local_nz),
202  KOKKOS_LAMBDA(const size_t nz, int &errorcnt) {
203  if (J_local_matrix.values(nz) != Jp_local_matrix.values(nz)) {
204  Kokkos::printf("Error in nonzero comparison %zu: %g != %g",
205  nz, J_local_matrix.values(nz), Jp_local_matrix.values(nz));
206  errorcnt++;
207  }
208  },
209  ierr);
210 
211 
212  if (ierr > 0) {
213  std::cout << testname << " FAILED on rank " << me << " with "
214  << (useBlock ? "Block maps" : "Cyclic maps")
215  << std::endl;
216  params.print();
217  }
218 
219  return (ierr == 0);
220  }
221 
222 private:
223 
225  // Return a map that is cyclic (like dealing rows to processors)
226  Teuchos::RCP<const map_t> getCyclicMap(
227  size_t nIndices,
228  Teuchos::Array<gno_t> &indices,
229  int mapNumProc,
230  const Teuchos::RCP<const Teuchos::Comm<int> > &comm)
231  {
232  size_t cnt = 0;
233  int me = comm->getRank();
234  int np = comm->getSize();
235  if (mapNumProc > np) mapNumProc = np; // corner case: bad input
236  if (mapNumProc <= 0) mapNumProc = 1; // corner case: np is too small
237 
238  for (size_t i = 0; i < nIndices; i++)
239  if (me == int(i % np)) indices[cnt++] = i;
240 
241  Tpetra::global_size_t dummy =
242  Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
243 
244  return rcp(new map_t(dummy, indices(0,cnt), 0, comm));
245  }
246 
248  // Input matrix -- built in Constructor
249  bool symmetric; // User can specify whether matrix is symmetric
250  Teuchos::RCP<matrix_t> JBlock; // has Trilinos default domain and range maps
251  Teuchos::RCP<matrix_t> JCyclic; // has cyclic domain and range maps
252  std::string matrixFileName;
253 };
254 
255 
257 int main(int narg, char **arg)
258 {
259  Tpetra::ScopeGuard scope(&narg, &arg);
260  Teuchos::RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
261  bool ok = true;
262  int ierr = 0;
263 
264  ColorerTest testColorer(comm, narg, arg);
265 
266  // Set parameters and compute coloring
267  {
268  Teuchos::ParameterList coloring_params;
269  std::string matrixType = "Jacobian";
270  bool symmetrize = true; // Zoltan's TRANSPOSE symmetrization, if needed
271 
272  coloring_params.set("MatrixType", matrixType);
273  coloring_params.set("symmetrize", symmetrize);
274 
275  ok = testColorer.run("Test One", coloring_params);
276  if (!ok) ierr++;
277  }
278 
279  {
280  Teuchos::ParameterList coloring_params;
281  std::string matrixType = "Jacobian";
282  bool symmetrize = false; // colorer's BIPARTITE symmetrization, if needed
283 
284  coloring_params.set("MatrixType", matrixType);
285  coloring_params.set("symmetrize", symmetrize);
286 
287  ok = testColorer.run("Test Two", coloring_params);
288  if (!ok) ierr++;
289  }
290 
291  {
292  Teuchos::ParameterList coloring_params;
293  std::string matrixType = "Jacobian";
294 
295  coloring_params.set("MatrixType", matrixType);
296 
297  ok = testColorer.run("Test Three", coloring_params);
298  if (!ok) ierr++;
299  }
300 
301  int gerr;
302  Teuchos::reduceAll<int, int>(*comm, Teuchos::REDUCE_SUM, 1, &ierr, &gerr);
303  if (comm->getRank() == 0) {
304  if (gerr == 0)
305  std::cout << "TEST PASSED" << std::endl;
306  else
307  std::cout << "TEST FAILED" << std::endl;
308  }
309 
310 //Through cmake...
311 //Test cases -- UserInputForTests can generate Galeri or read files:
312 //- tri-diagonal matrix -- can check the number of colors
313 //- galeri matrix
314 //- read from file: symmetric matrix and non-symmetric matrix
315 
316 //Through code ...
317 //Test with fitted and non-fitted maps
318 //Call regular and fitted versions of functions
319 
320 //Through code ...
321 //Test both with and without Symmetrize --
322 //test both to exercise both sets of callbacks in Zoltan
323 // --matrixType = Jacobian/Hessian
324 // --symmetric, --no-symmetric
325 // --symmetrize, --no-symmetrize
326 
327 //Through cmake
328 //Test both with and without distributeInput
329 // --distribute, --no-distribute
330 
331 }
Tpetra::CrsMatrix< zscalar_t > matrix_t
Definition: Bug9500.cpp:21
typename matrix_t::device_type::execution_space execution_space_t
Definition: Bug9500.cpp:23
RCP< tcrsMatrix_t > getUITpetraCrsMatrix()
ColorerTest(Teuchos::RCP< const Teuchos::Comm< int > > &comm, int narg, char **arg)
int main(int narg, char **arg)
Definition: coloring1.cpp:164
common code used by tests
Tpetra::MultiVector< zscalar_t > multivector_t
Definition: Bug9500.cpp:22
bool run(const char *testname, Teuchos::ParameterList &params)
typename map_t::global_ordinal_type gno_t
Definition: Bug9500.cpp:19
bool buildAndCheckSeedMatrix(const char *testname, Teuchos::ParameterList &params, const bool useBlock)
Definition: Bug9500.cpp:126
Tpetra::global_size_t global_size_t
Tpetra::CrsGraph<> graph_t
Definition: Bug9500.cpp:20
Tpetra::Map<> map_t
Definition: Bug9500.cpp:18