52 #include <Teuchos_ParameterList.hpp>
53 #include <Teuchos_RCP.hpp>
54 #include <Teuchos_FancyOStream.hpp>
55 #include <Teuchos_CommandLineProcessor.hpp>
56 #include <Tpetra_CrsMatrix.hpp>
57 #include <Tpetra_Vector.hpp>
58 #include <MatrixMarket_Tpetra.hpp>
80 typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO>
SparseMatrix;
82 typedef Tpetra::Vector<z2TestScalar, z2TestLO, z2TestGO>
Vector;
83 typedef Vector::node_type
Node;
91 typedef Tpetra::Vector<int, z2TestLO, z2TestGO>
IntVector;
94 #define epsilon 0.00000001
98 int main(
int narg,
char** arg)
100 std::string inputFile =
"";
101 std::string outputFile =
"";
103 std::string method =
"scotch";
104 bool verbose =
false;
105 bool distributeInput =
true;
106 bool haveFailure =
false;
113 Tpetra::ScopeGuard tscope(&narg, &arg);
114 RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
115 int me = comm->getRank();
118 Teuchos::CommandLineProcessor cmdp (
false,
false);
119 cmdp.setOption(
"inputPath", &inputPath,
120 "Path to the MatrixMarket or Zoltan file to be read; "
121 "if not specified, a default path will be used.");
122 cmdp.setOption(
"inputFile", &inputFile,
123 "Name of the Matrix Market or Zoltan file to read; "
124 "if not specified, a matrix will be generated by MueLu.");
125 cmdp.setOption(
"outputFile", &outputFile,
126 "Name of the Matrix Market sparse matrix file to write, "
127 "echoing the input/generated matrix.");
128 cmdp.setOption(
"method", &method,
129 "Partitioning method to use: scotch or parmetis.");
130 cmdp.setOption(
"nparts", &nParts,
131 "Number of parts being requested");
132 cmdp.setOption(
"vertexWeights", &nVwgts,
133 "Number of weights to generate for each vertex");
134 cmdp.setOption(
"edgeWeights", &nEwgts,
135 "Number of weights to generate for each edge");
136 cmdp.setOption(
"verbose",
"quiet", &verbose,
137 "Print messages and results.");
138 cmdp.setOption(
"distribute",
"no-distribute", &distributeInput,
139 "indicate whether or not to distribute "
140 "input across the communicator");
151 std::string matrixType(
"Laplace3D");
153 cmdp.setOption(
"x", &xdim,
154 "number of gridpoints in X dimension for "
155 "mesh used to generate matrix.");
156 cmdp.setOption(
"y", &ydim,
157 "number of gridpoints in Y dimension for "
158 "mesh used to generate matrix.");
159 cmdp.setOption(
"z", &zdim,
160 "number of gridpoints in Z dimension for "
161 "mesh used to generate matrix.");
162 cmdp.setOption(
"matrix", &matrixType,
163 "Matrix type: Laplace1D, Laplace2D, or Laplace3D");
167 int quotientThreshold = -1;
168 cmdp.setOption(
"qthreshold", "ientThreshold,
169 "Threshold on the number of vertices for active MPI ranks to hold"
170 "after the migrating the communication graph to the active ranks.");
174 cmdp.parse(narg, arg);
176 RCP<UserInputForTests> uinput;
180 true, distributeInput));
184 true, distributeInput));
186 RCP<SparseMatrix> origMatrix = uinput->getUITpetraCrsMatrix();
188 if (origMatrix->getGlobalNumRows() < 40) {
189 Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,
false));
190 origMatrix->describe(out, Teuchos::VERB_EXTREME);
194 if (outputFile !=
"") {
196 Tpetra::MatrixMarket::Writer<SparseMatrix>::writeSparseFile(outputFile,
197 origMatrix, verbose);
201 std::cout <<
"NumRows = " << origMatrix->getGlobalNumRows() << std::endl
202 <<
"NumNonzeros = " << origMatrix->getGlobalNumEntries() << std::endl
203 <<
"NumProcs = " << comm->getSize() << std::endl
204 <<
"NumLocalRows (rank 0) = " << origMatrix->getLocalNumRows() << std::endl;
207 RCP<Vector> origVector, origProd;
208 origProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
209 origMatrix->getRangeMap());
210 origVector = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
211 origMatrix->getDomainMap());
212 origVector->randomize();
215 Teuchos::ParameterList params;
217 params.set(
"partitioning_approach",
"partition");
218 params.set(
"algorithm", method);
222 params.set(
"num_global_parts", nParts);
226 if(method ==
"quotient" && quotientThreshold > 0) {
227 params.set(
"quotient_threshold", quotientThreshold);
240 size_t nrows = origMatrix->getLocalNumRows();
243 for (
size_t i = 0; i < nrows; i++) {
244 size_t idx = i * nVwgts;
245 vwgts[idx] =
zscalar_t(origMatrix->getRowMap()->getGlobalElement(i))
247 for (
int j = 1; j < nVwgts; j++) vwgts[idx+j] = 1.;
249 for (
int j = 0; j < nVwgts; j++) {
250 if (j !=
NNZ_IDX) adapter.setVertexWeights(&vwgts[j], nVwgts, j);
251 else adapter.setVertexWeightIsDegree(
NNZ_IDX);
258 size_t nnz = origMatrix->getLocalNumEntries();
260 size_t nrows = origMatrix->getLocalNumRows();
261 size_t maxnzrow = origMatrix->getLocalMaxNumRowEntries();
264 typename SparseMatrix::nonconst_global_inds_host_view_type egids(
"egids", maxnzrow);
265 typename SparseMatrix::nonconst_values_host_view_type evals(
"evals", maxnzrow);
266 for (
size_t i = 0; i < nrows; i++) {
268 z2TestGO gid = origMatrix->getRowMap()->getGlobalElement(i);
269 origMatrix->getGlobalRowCopy(gid, egids, evals, nnzinrow);
270 for (
size_t k = 0; k < nnzinrow; k++) {
271 ewgts[cnt] = (gid < egids[k] ? gid : egids[k]);
272 if (nEwgts > 1) ewgts[cnt+nnz] = (gid < egids[k] ? egids[k] : gid);
273 for (
int j = 2; j < nEwgts; j++) ewgts[cnt+nnz*j] = 1.;
277 for (
int j = 0; j < nEwgts; j++) {
278 adapter.setEdgeWeights(&ewgts[j*nnz], 1, j);
288 if (me == 0) std::cout <<
"Calling solve() " << std::endl;
292 if (me == 0) std::cout <<
"Done solve() " << std::endl;
294 catch (std::runtime_error &e) {
297 std::cout <<
"Runtime exception returned from solve(): " << e.what();
298 if (!strncmp(e.what(),
"BUILD ERROR", 11)) {
300 std::cout <<
" PASS" << std::endl;
305 std::cout <<
" FAIL" << std::endl;
309 catch (std::logic_error &e) {
312 std::cout <<
"Logic exception returned from solve(): " << e.what()
313 <<
" FAIL" << std::endl;
316 catch (std::bad_alloc &e) {
319 std::cout <<
"Bad_alloc exception returned from solve(): " << e.what()
320 <<
" FAIL" << std::endl;
323 catch (std::exception &e) {
326 std::cout <<
"Unknown exception returned from solve(). " << e.what()
327 <<
" FAIL" << std::endl;
333 size_t checkNparts = comm->getSize();
334 if(nParts != -1) checkNparts = size_t(nParts);
335 size_t checkLength = origMatrix->getLocalNumRows();
340 size_t *countPerPart =
new size_t[checkNparts];
341 size_t *globalCountPerPart =
new size_t[checkNparts];
344 for (
size_t i = 0; i < checkNparts; i++) countPerPart[i] = 0;
345 for (
size_t i = 0; i < checkNparts * nVwgts; i++) wtPerPart[i] = 0.;
347 for (
size_t i = 0; i < checkLength; i++) {
348 if (
size_t(checkParts[i]) >= checkNparts)
349 std::cout <<
"Invalid Part " << checkParts[i] <<
": FAIL" << std::endl;
350 countPerPart[checkParts[i]]++;
351 for (
int j = 0; j < nVwgts; j++) {
353 wtPerPart[checkParts[i]*nVwgts+j] += vwgts[i*nVwgts+j];
355 wtPerPart[checkParts[i]*nVwgts+j] += origMatrix->getNumEntriesInLocalRow(i);
360 if(method ==
"quotient") {
361 size_t result = size_t(checkParts[0]);
362 for (
size_t i = 1; i < checkLength; i++) {
363 if (
size_t(checkParts[i]) != result)
364 std::cout <<
"Different parts in the quotient algorithm: "
365 << result <<
"!=" << checkParts[i] <<
": FAIL" << std::endl;
370 Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, checkNparts,
371 countPerPart, globalCountPerPart);
372 Teuchos::reduceAll<int, zscalar_t>(*comm, Teuchos::REDUCE_SUM,
374 wtPerPart, globalWtPerPart);
376 size_t min = std::numeric_limits<std::size_t>::max();
379 size_t minrank = 0, maxrank = 0;
380 for (
size_t i = 0; i < checkNparts; i++) {
381 if (globalCountPerPart[i] < min) {min = globalCountPerPart[i]; minrank = i;}
382 if (globalCountPerPart[i] > max) {max = globalCountPerPart[i]; maxrank = i;}
383 sum += globalCountPerPart[i];
387 float avg = (float) sum / (
float) checkNparts;
388 std::cout <<
"Minimum count: " << min <<
" on rank " << minrank << std::endl;
389 std::cout <<
"Maximum count: " << max <<
" on rank " << maxrank << std::endl;
390 std::cout <<
"Average count: " << avg << std::endl;
391 std::cout <<
"Total count: " << sum
392 << (sum != origMatrix->getGlobalNumRows()
393 ?
"Work was lost; FAIL"
396 std::cout <<
"Imbalance: " << max / avg << std::endl;
398 std::vector<zscalar_t> minwt(nVwgts, std::numeric_limits<zscalar_t>::max());
399 std::vector<zscalar_t> maxwt(nVwgts, 0.);
400 std::vector<zscalar_t> sumwt(nVwgts, 0.);
401 for (
size_t i = 0; i < checkNparts; i++) {
402 for (
int j = 0; j < nVwgts; j++) {
403 size_t idx = i*nVwgts+j;
404 if (globalWtPerPart[idx] < minwt[j]) minwt[j] = globalWtPerPart[idx];
405 if (globalWtPerPart[idx] > maxwt[j]) maxwt[j] = globalWtPerPart[idx];
406 sumwt[j] += globalWtPerPart[idx];
409 for (
int j = 0; j < nVwgts; j++) {
410 float avgwt = (float) sumwt[j] / (
float) checkNparts;
411 std::cout << std::endl;
412 std::cout <<
"Minimum weight[" << j <<
"]: " << minwt[j] << std::endl;
413 std::cout <<
"Maximum weight[" << j <<
"]: " << maxwt[j] << std::endl;
414 std::cout <<
"Average weight[" << j <<
"]: " << avgwt << std::endl;
415 std::cout <<
"Imbalance: " << maxwt[j] / avgwt << std::endl;
420 delete [] countPerPart;
422 delete [] globalCountPerPart;
423 delete [] globalWtPerPart;
429 if (me == 0) std::cout <<
"Redistributing matrix..." << std::endl;
434 if (redistribMatrix->getGlobalNumRows() < 40) {
435 Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,
false));
436 redistribMatrix->describe(out, Teuchos::VERB_EXTREME);
439 if (me == 0) std::cout <<
"Redistributing vectors..." << std::endl;
447 RCP<Vector> redistribProd;
448 redistribProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
449 redistribMatrix->getRangeMap());
453 RCP<IntVector> origIntVec;
455 origIntVec = Tpetra::createVector<int,z2TestLO,z2TestGO>(
456 origMatrix->getRangeMap());
457 for (
size_t i = 0; i < origIntVec->getLocalLength(); i++)
458 origIntVec->replaceLocalValue(i, me);
463 int origIntNorm = origIntVec->norm1();
464 int redistIntNorm = redistIntVec->norm1();
465 if (me == 0) std::cout <<
"IntegerVectorTest: " << origIntNorm <<
" == "
466 << redistIntNorm <<
" ?";
467 if (origIntNorm != redistIntNorm) {
468 if (me == 0) std::cout <<
" FAIL" << std::endl;
471 else if (me == 0) std::cout <<
" OK" << std::endl;
477 if (me == 0) std::cout <<
"Matvec original..." << std::endl;
478 origMatrix->apply(*origVector, *origProd);
481 std::cout <<
"Norm of Original matvec prod: " << origNorm << std::endl;
483 if (me == 0) std::cout <<
"Matvec redistributed..." << std::endl;
484 redistribMatrix->apply(*redistribVector, *redistribProd);
487 std::cout <<
"Norm of Redistributed matvec prod: " << redistribNorm << std::endl;
489 if (redistribNorm > origNorm+
epsilon || redistribNorm < origNorm-
epsilon) {
494 delete redistribVector;
495 delete redistribMatrix;
499 std::cout <<
"Mat-Vec product changed; FAIL" << std::endl;
503 std::cout <<
"PASS" << std::endl;
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
Provides access for Zoltan2 to Xpetra::CrsGraph data.
Zoltan2::XpetraCrsGraphAdapter< SparseGraph > SparseGraphAdapter
int main(int narg, char **arg)
common code used by tests
typename InputTraits< User >::part_t part_t
Tpetra::Vector< z2TestScalar, z2TestLO, z2TestGO > Vector
Defines the XpetraMultiVectorAdapter.
Defines XpetraCrsGraphAdapter class.
Defines the XpetraCrsMatrixAdapter class.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
An adapter for Xpetra::MultiVector.
Tpetra::Map::local_ordinal_type zlno_t
Zoltan2::XpetraMultiVectorAdapter< IntVector > IntVectorAdapter
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
Tpetra::CrsGraph< z2TestLO, z2TestGO > SparseGraph
PartitioningProblem sets up partitioning problems for the user.
Zoltan2::XpetraCrsMatrixAdapter< SparseMatrix > SparseMatrixAdapter
Tpetra::Vector< int, z2TestLO, z2TestGO > IntVector
Defines the PartitioningProblem class.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
Zoltan2::XpetraMultiVectorAdapter< Vector > MultiVectorAdapter
Tpetra::Map::global_ordinal_type zgno_t
void solve(bool updateInputData=true)
Direct the problem to create a solution.
std::string testDataFilePath(".")