2 #include "Tpetra_Core.hpp"
3 #include "Kokkos_Random.hpp"
11 using gno_t =
typename map_t::global_ordinal_type;
20 ColorerTest(
const Teuchos::RCP<
const Teuchos::Comm<int> > &comm,
int multiple)
22 int me = comm->getRank();
23 int np = comm->getSize();
27 Teuchos::Array<gno_t> myRows(myNrows);
28 for (
size_t i = 0; i < myNrows; i++) {
29 myRows[i] = multiple * (me * myNrows + i);
33 Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
34 Teuchos::RCP<const map_t> map = rcp(
new map_t(dummy, myRows, 0, comm));
37 JBlock = rcp(
new matrix_t(map, nnz));
38 Teuchos::Array<gno_t> myCols(nnz);
39 Teuchos::Array<zscalar_t> myVals(nnz);
41 for (
size_t i = 0; i < myNrows; i++) {
42 auto gid = map->getGlobalElement(i);
45 if (gid+multiple <= map->getMaxAllGlobalIndex())
46 myCols[cnt++] = gid+multiple;
47 JBlock->insertGlobalValues(gid, myCols(0,cnt), myVals(0, cnt));
49 JBlock->fillComplete();
52 using IST =
typename Kokkos::ArithTraits<zscalar_t>::val_type;
54 Kokkos::Random_XorShift64_Pool<execution_space_t>;
55 pool_type rand_pool(static_cast<uint64_t>(me));
57 Kokkos::fill_random(JBlock->getLocalMatrixDevice().values, rand_pool,
58 static_cast<IST
>(1.),
static_cast<IST
>(9999.));
60 Teuchos::FancyOStream foo(Teuchos::rcp(&std::cout,
false));
61 JBlock->describe(foo, Teuchos::VERB_EXTREME);
67 size_t nIndices = std::max(JBlock->getGlobalNumCols(),
68 JBlock->getGlobalNumRows());
69 Teuchos::Array<gno_t> indices(nIndices);
71 Teuchos::RCP<const map_t> vMapCyclic =
72 getCyclicMap(JBlock->getGlobalNumCols(), indices, np-1,
74 Teuchos::RCP<const map_t> wMapCyclic =
75 getCyclicMap(JBlock->getGlobalNumRows(), indices, np-2,
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());
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);
93 JCyclic->fillComplete();
94 JCyclic->describe(foo, Teuchos::VERB_EXTREME);
98 bool run(
const char* testname, Teuchos::ParameterList ¶ms) {
102 params.set(
"symmetric",
false);
115 const char *testname,
116 Teuchos::ParameterList ¶ms,
123 Teuchos::RCP<matrix_t> J = (useBlock ? JBlock : JCyclic);
124 int me = J->getRowMap()->getComm()->getRank();
126 std::cout <<
"Running " << testname <<
" with "
127 << (useBlock ?
"Block maps" :
"Cyclic maps")
132 colorer.computeColoring(params);
135 if (!colorer.checkColoring()) {
136 std::cout << testname <<
" with "
137 << (useBlock ?
"Block maps" :
"Cyclic maps")
138 <<
" FAILED: invalid coloring returned"
146 const int numColors = colorer.getNumColors();
151 colorer.computeSeedMatrix(V);
160 Teuchos::RCP<matrix_t> Jp = rcp(
new matrix_t(*J, Teuchos::Copy));
161 Jp->setAllToScalar(static_cast<zscalar_t>(-1.));
163 colorer.reconstructMatrix(W, *Jp);
166 auto J_local_matrix = J->getLocalMatrixDevice();
167 auto Jp_local_matrix = Jp->getLocalMatrixDevice();
168 const size_t num_local_nz = J->getLocalNumEntries();
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 Kokkos::printf(
"Error in nonzero comparison %zu: %g != %g",
176 nz, J_local_matrix.values(nz), Jp_local_matrix.values(nz));
184 std::cout << testname <<
" FAILED on rank " << me <<
" with "
185 << (useBlock ?
"Block maps" :
"Cyclic maps")
197 Teuchos::RCP<const map_t> getCyclicMap(
199 Teuchos::Array<gno_t> &indices,
202 const Teuchos::RCP<
const Teuchos::Comm<int> > &comm)
205 int me = comm->getRank();
206 int np = comm->getSize();
207 if (mapNumProc > np) mapNumProc = np;
208 if (mapNumProc <= 0) mapNumProc = 1;
210 for (
size_t i = 0; i < nIndices; i++)
211 if (me ==
int(i % np)) indices[cnt++] = multiple*i;
214 Teuchos::OrdinalTraits<Tpetra::global_size_t>::invalid();
216 return rcp(
new map_t(dummy, indices(0,cnt), 0, comm));
222 Teuchos::RCP<matrix_t> JBlock;
223 Teuchos::RCP<matrix_t> JCyclic;
228 int doTheTest(
const Teuchos::RCP<
const Teuchos::Comm<int> > &comm,
int multiple)
237 Teuchos::ParameterList coloring_params;
238 std::string matrixType =
"Jacobian";
239 bool symmetrize =
true;
241 coloring_params.set(
"MatrixType", matrixType);
242 coloring_params.set(
"symmetrize", symmetrize);
244 ok = testColorer.
run(
"Test One", coloring_params);
249 Teuchos::ParameterList coloring_params;
250 std::string matrixType =
"Jacobian";
251 bool symmetrize =
false;
253 coloring_params.set(
"MatrixType", matrixType);
254 coloring_params.set(
"symmetrize", symmetrize);
256 ok = testColorer.
run(
"Test Two", coloring_params);
261 Teuchos::ParameterList coloring_params;
262 std::string matrixType =
"Jacobian";
264 coloring_params.set(
"MatrixType", matrixType);
266 ok = testColorer.
run(
"Test Three", coloring_params);
275 Tpetra::ScopeGuard scope(&narg, &arg);
276 Teuchos::RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
285 Teuchos::reduceAll<int, int>(*comm, Teuchos::REDUCE_SUM, 1, &ierr, &gerr);
286 if (comm->getRank() == 0) {
288 std::cout <<
"TEST PASSED" << std::endl;
290 std::cout <<
"TEST FAILED" << std::endl;
Tpetra::CrsMatrix< zscalar_t > matrix_t
int doTheTest(const Teuchos::RCP< const Teuchos::Comm< int > > &comm, int multiple)
typename matrix_t::device_type::execution_space execution_space_t
int main(int narg, char **arg)
common code used by tests
Tpetra::MultiVector< zscalar_t > multivector_t
ColorerTest(const Teuchos::RCP< const Teuchos::Comm< int > > &comm, int multiple)
bool run(const char *testname, Teuchos::ParameterList ¶ms)
typename map_t::global_ordinal_type gno_t
bool buildAndCheckSeedMatrix(const char *testname, Teuchos::ParameterList ¶ms, const bool useBlock)
Tpetra::global_size_t global_size_t
Tpetra::CrsGraph<> graph_t