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