Zoltan2
Loading...
Searching...
No Matches
Bug9500.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
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
20 ColorerTest(const Teuchos::RCP<const Teuchos::Comm<int> > &comm, int multiple)
21 {
22 int me = comm->getRank();
23 int np = comm->getSize();
24
25 // Create non-symmetrix matrix with non-contiguous row map -- only even GIDs
26 size_t myNrows = 4;
27 Teuchos::Array<gno_t> myRows(myNrows);
28 for (size_t i = 0; i < myNrows; i++) {
29 myRows[i] = multiple * (me * myNrows + i);
30 }
31
32 Tpetra::global_size_t dummy =
33 Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
34 Teuchos::RCP<const map_t> map = rcp(new map_t(dummy, myRows, 0, comm));
35
36 size_t nnz = 2;
37 JBlock = rcp(new matrix_t(map, nnz));
38 Teuchos::Array<gno_t> myCols(nnz);
39 Teuchos::Array<zscalar_t> myVals(nnz);
40
41 for (size_t i = 0; i < myNrows; i++) {
42 auto gid = map->getGlobalElement(i);
43 size_t cnt = 0;
44 myCols[cnt++] = gid;
45 if (gid+multiple <= map->getMaxAllGlobalIndex())
46 myCols[cnt++] = gid+multiple;
47 JBlock->insertGlobalValues(gid, myCols(0,cnt), myVals(0, cnt));
48 }
49 JBlock->fillComplete();
50
51 // Fill JBlock with random numbers for a better test.
52 using IST = typename Kokkos::Details::ArithTraits<zscalar_t>::val_type;
53 using pool_type =
54 Kokkos::Random_XorShift64_Pool<execution_space_t>;
55 pool_type rand_pool(static_cast<uint64_t>(me));
56
57 Kokkos::fill_random(JBlock->getLocalMatrixDevice().values, rand_pool,
58 static_cast<IST>(1.), static_cast<IST>(9999.));
59
60 Teuchos::FancyOStream foo(Teuchos::rcp(&std::cout,false));
61 JBlock->describe(foo, Teuchos::VERB_EXTREME);
62
63 // Build same matrix with cyclic domain and range maps
64 // To make range and domain maps differ for square matrices,
65 // keep some processors empty in the cyclic maps
66
67 size_t nIndices = std::max(JBlock->getGlobalNumCols(),
68 JBlock->getGlobalNumRows());
69 Teuchos::Array<gno_t> indices(nIndices);
70
71 Teuchos::RCP<const map_t> vMapCyclic =
72 getCyclicMap(JBlock->getGlobalNumCols(), indices, np-1,
73 multiple, comm);
74 Teuchos::RCP<const map_t> wMapCyclic =
75 getCyclicMap(JBlock->getGlobalNumRows(), indices, np-2,
76 multiple, comm);
77
78 // Make JCyclic: same matrix with different Domain and Range maps
79 RCP<const graph_t> block_graph = JBlock->getCrsGraph();
80 RCP<graph_t> cyclic_graph = rcp(new graph_t(*block_graph));
81 cyclic_graph->resumeFill();
82 cyclic_graph->fillComplete(vMapCyclic, wMapCyclic);
83 JCyclic = rcp(new matrix_t(cyclic_graph));
84 JCyclic->resumeFill();
85 TEUCHOS_ASSERT(block_graph->getLocalNumRows() ==
86 cyclic_graph->getLocalNumRows());
87 {
88 auto val_s = JBlock->getLocalMatrixHost().values;
89 auto val_d = JCyclic->getLocalMatrixHost().values;
90 TEUCHOS_ASSERT(val_s.extent(0) == val_d.extent(0));
91 Kokkos::deep_copy(val_d, val_s);
92 }
93 JCyclic->fillComplete();
94 JCyclic->describe(foo, Teuchos::VERB_EXTREME);
95 }
96
98 bool run(const char* testname, Teuchos::ParameterList &params) {
99
100 bool ok = true;
101
102 params.set("symmetric", false);
103
104 // test with default maps
105 ok = buildAndCheckSeedMatrix(testname, params, true);
106
107 // test with cyclic maps
108 ok &= buildAndCheckSeedMatrix(testname, params, false);
109
110 return ok;
111 }
112
115 const char *testname,
116 Teuchos::ParameterList &params,
117 const bool useBlock
118 )
119 {
120 int ierr = 0;
121
122 // Pick matrix depending on useBlock flag
123 Teuchos::RCP<matrix_t> J = (useBlock ? JBlock : JCyclic);
124 int me = J->getRowMap()->getComm()->getRank();
125
126 std::cout << "Running " << testname << " with "
127 << (useBlock ? "Block maps" : "Cyclic maps")
128 << std::endl;
129
130 // Create a colorer
132 colorer.computeColoring(params);
133
134 // Check coloring
135 if (!colorer.checkColoring()) {
136 std::cout << testname << " with "
137 << (useBlock ? "Block maps" : "Cyclic maps")
138 << " FAILED: invalid coloring returned"
139 << std::endl;
140 return false;
141 }
142
143 // Compute seed matrix V -- the application wants this matrix
144 // Dense matrix of 0/1 indicating the compression via coloring
145
146 const int numColors = colorer.getNumColors();
147
148 // Compute the seed matrix; applications want this seed matrix
149
150 multivector_t V(J->getDomainMap(), numColors);
151 colorer.computeSeedMatrix(V);
152
153 // To test the result...
154 // Compute the compressed matrix W
155 multivector_t W(J->getRangeMap(), numColors);
156
157 J->apply(V, W);
158
159 // Reconstruct matrix from compression vector
160 Teuchos::RCP<matrix_t> Jp = rcp(new matrix_t(*J, Teuchos::Copy));
161 Jp->setAllToScalar(static_cast<zscalar_t>(-1.));
162
163 colorer.reconstructMatrix(W, *Jp);
164
165 // Check that values of J = values of Jp
166 auto J_local_matrix = J->getLocalMatrixDevice();
167 auto Jp_local_matrix = Jp->getLocalMatrixDevice();
168 const size_t num_local_nz = J->getLocalNumEntries();
169
170 Kokkos::parallel_reduce(
171 "TpetraCrsColorer::testReconstructedMatrix()",
172 Kokkos::RangePolicy<execution_space_t>(0, num_local_nz),
173 KOKKOS_LAMBDA(const size_t nz, int &errorcnt) {
174 if (J_local_matrix.values(nz) != Jp_local_matrix.values(nz)) {
175 printf("Error in nonzero comparison %zu: %g != %g",
176 nz, J_local_matrix.values(nz), Jp_local_matrix.values(nz));
177 errorcnt++;
178 }
179 },
180 ierr);
181
182
183 if (ierr > 0) {
184 std::cout << testname << " FAILED on rank " << me << " with "
185 << (useBlock ? "Block maps" : "Cyclic maps")
186 << std::endl;
187 params.print();
188 }
189
190 return (ierr == 0);
191 }
192
193private:
194
196 // Return a map that is cyclic (like dealing rows to processors)
197 Teuchos::RCP<const map_t> getCyclicMap(
198 size_t nIndices,
199 Teuchos::Array<gno_t> &indices,
200 int mapNumProc,
201 int multiple,
202 const Teuchos::RCP<const Teuchos::Comm<int> > &comm)
203 {
204 size_t cnt = 0;
205 int me = comm->getRank();
206 int np = comm->getSize();
207 if (mapNumProc > np) mapNumProc = np; // corner case: bad input
208 if (mapNumProc <= 0) mapNumProc = 1; // corner case: np is too small
209
210 for (size_t i = 0; i < nIndices; i++)
211 if (me == int(i % np)) indices[cnt++] = multiple*i;
212
213 Tpetra::global_size_t dummy =
214 Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
215
216 return rcp(new map_t(dummy, indices(0,cnt), 0, comm));
217 }
218
220 // Input matrix -- built in Constructor
221 bool symmetric; // User can specify whether matrix is symmetric
222 Teuchos::RCP<matrix_t> JBlock; // has Trilinos default domain and range maps
223 Teuchos::RCP<matrix_t> JCyclic; // has cyclic domain and range maps
224};
225
226
228int doTheTest(const Teuchos::RCP<const Teuchos::Comm<int> > &comm, int multiple)
229{
230 bool ok = true;
231 int ierr = 0;
232
233 ColorerTest testColorer(comm, multiple);
234
235 // Set parameters and compute coloring
236 {
237 Teuchos::ParameterList coloring_params;
238 std::string matrixType = "Jacobian";
239 bool symmetrize = true; // Zoltan's TRANSPOSE symmetrization, if needed
240
241 coloring_params.set("MatrixType", matrixType);
242 coloring_params.set("symmetrize", symmetrize);
243
244 ok = testColorer.run("Test One", coloring_params);
245 if (!ok) ierr++;
246 }
247
248 {
249 Teuchos::ParameterList coloring_params;
250 std::string matrixType = "Jacobian";
251 bool symmetrize = false; // colorer's BIPARTITE symmetrization, if needed
252
253 coloring_params.set("MatrixType", matrixType);
254 coloring_params.set("symmetrize", symmetrize);
255
256 ok = testColorer.run("Test Two", coloring_params);
257 if (!ok) ierr++;
258 }
259
260 {
261 Teuchos::ParameterList coloring_params;
262 std::string matrixType = "Jacobian";
263
264 coloring_params.set("MatrixType", matrixType);
265
266 ok = testColorer.run("Test Three", coloring_params);
267 if (!ok) ierr++;
268 }
269 return ierr;
270}
271
273int main(int narg, char **arg)
274{
275 Tpetra::ScopeGuard scope(&narg, &arg);
276 Teuchos::RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
277 int ierr = 0;
278
279 ierr += doTheTest(comm, 1); // Contiguous row map
280 ierr += doTheTest(comm, 2); // Non-contiguous row map --
281 // only even-numbered indices
282 ierr += doTheTest(comm, 5); // Indices spaced wider than rows/proc
283
284 int gerr;
285 Teuchos::reduceAll<int, int>(*comm, Teuchos::REDUCE_SUM, 1, &ierr, &gerr);
286 if (comm->getRank() == 0) {
287 if (gerr == 0)
288 std::cout << "TEST PASSED" << std::endl;
289 else
290 std::cout << "TEST FAILED" << std::endl;
291 }
292}
int doTheTest(const Teuchos::RCP< const Teuchos::Comm< int > > &comm, int multiple)
Definition: Bug9500.cpp:228
common code used by tests
float zscalar_t
int main()
ColorerTest(const Teuchos::RCP< const Teuchos::Comm< int > > &comm, int multiple)
Definition: Bug9500.cpp:20
bool run(const char *testname, Teuchos::ParameterList &params)
Definition: Bug9500.cpp:98
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