Zoltan2
Loading...
Searching...
No Matches
partitioning1.cpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Zoltan2: A package of combinatorial algorithms for scientific computing
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Karen Devine (kddevin@sandia.gov)
39// Erik Boman (egboman@sandia.gov)
40// Siva Rajamanickam (srajama@sandia.gov)
41//
42// ***********************************************************************
43//
44// @HEADER
50#include <iostream>
51#include <limits>
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>
59
60using Teuchos::RCP;
61
63// Program to demonstrate use of Zoltan2 to partition a TPetra matrix
64// (read from a MatrixMarket file or generated by Galeri::Xpetra).
65// Usage:
66// a.out [--inputFile=filename] [--outputFile=outfile] [--verbose]
67// [--x=#] [--y=#] [--z=#] [--matrix={Laplace1D,Laplace2D,Laplace3D}
68// Karen Devine, 2011
70
72// Eventually want to use Teuchos unit tests to vary z2TestLO and
73// GO. For now, we set them at compile time based on whether Tpetra
74// is built with explicit instantiation on. (in Zoltan2_TestHelpers.hpp)
75
79
80typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO> SparseMatrix;
81typedef Tpetra::CrsGraph<z2TestLO, z2TestGO> SparseGraph;
82typedef Tpetra::Vector<z2TestScalar, z2TestLO, z2TestGO> Vector;
83typedef Vector::node_type Node;
84
88
89
90// Integer vector
91typedef Tpetra::Vector<int, z2TestLO, z2TestGO> IntVector;
93
94#define epsilon 0.00000001
95#define NNZ_IDX 1
96
98int main(int narg, char** arg)
99{
100 std::string inputFile = ""; // Matrix Market or Zoltan file to read
101 std::string outputFile = ""; // Matrix Market or Zoltan file to write
102 std::string inputPath = testDataFilePath; // Directory with input file
103 std::string method = "scotch";
104 bool verbose = false; // Verbosity of output
105 bool distributeInput = true;
106 bool haveFailure = false;
107 int nParts = -1;
108 int nVwgts = 0;
109 int nEwgts = 0;
110 int testReturn = 0;
111
113 Tpetra::ScopeGuard tscope(&narg, &arg);
114 RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
115 int me = comm->getRank();
116
117 // Read run-time options.
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");
141
143 // Even with cmdp option "true", I get errors for having these
144 // arguments on the command line. (On redsky build)
145 // KDDKDD Should just be warnings, right? Code should still work with these
146 // KDDKDD params in the create-a-matrix file. Better to have them where
147 // KDDKDD they are used.
148 int xdim=10;
149 int ydim=10;
150 int zdim=10;
151 std::string matrixType("Laplace3D");
152
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");
164
166 // Quotient-specific parameters
167 int quotientThreshold = -1;
168 cmdp.setOption("qthreshold", &quotientThreshold,
169 "Threshold on the number of vertices for active MPI ranks to hold"
170 "after the migrating the communication graph to the active ranks.");
171
173
174 cmdp.parse(narg, arg);
175
176 RCP<UserInputForTests> uinput;
177
178 if (inputFile != "") // Input file specified; read a matrix
179 uinput = rcp(new UserInputForTests(inputPath, inputFile, comm,
180 true, distributeInput));
181
182 else // Let MueLu generate a default matrix
183 uinput = rcp(new UserInputForTests(xdim, ydim, zdim, string(""), comm,
184 true, distributeInput));
185
186 RCP<SparseMatrix> origMatrix = uinput->getUITpetraCrsMatrix();
187
188 if (origMatrix->getGlobalNumRows() < 40) {
189 Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,false));
190 origMatrix->describe(out, Teuchos::VERB_EXTREME);
191 }
192
193
194 if (outputFile != "") {
195 // Just a sanity check.
196 Tpetra::MatrixMarket::Writer<SparseMatrix>::writeSparseFile(outputFile,
197 origMatrix, verbose);
198 }
199
200 if (me == 0)
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;
205
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();
213
215 Teuchos::ParameterList params;
216
217 params.set("partitioning_approach", "partition");
218 params.set("algorithm", method);
219
221 if(nParts > 0) {
222 params.set("num_global_parts", nParts);
223 }
224
226 if(method == "quotient" && quotientThreshold > 0) {
227 params.set("quotient_threshold", quotientThreshold);
228 }
229
231 SparseGraphAdapter adapter(origMatrix->getCrsGraph(), nVwgts, nEwgts);
232
236
237 zscalar_t *vwgts = NULL, *ewgts = NULL;
238 if (nVwgts) {
239 // Test vertex weights with stride nVwgts.
240 size_t nrows = origMatrix->getLocalNumRows();
241 if (nrows) {
242 vwgts = new zscalar_t[nVwgts * nrows];
243 for (size_t i = 0; i < nrows; i++) {
244 size_t idx = i * nVwgts;
245 vwgts[idx] = zscalar_t(origMatrix->getRowMap()->getGlobalElement(i))
246 ;// + zscalar_t(0.5);
247 for (int j = 1; j < nVwgts; j++) vwgts[idx+j] = 1.;
248 }
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);
252 }
253 }
254 }
255
256 if (nEwgts) {
257 // Test edge weights with stride 1.
258 size_t nnz = origMatrix->getLocalNumEntries();
259 if (nnz) {
260 size_t nrows = origMatrix->getLocalNumRows();
261 size_t maxnzrow = origMatrix->getLocalMaxNumRowEntries();
262 ewgts = new zscalar_t[nEwgts * nnz];
263 size_t cnt = 0;
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++) {
267 size_t nnzinrow;
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.;
274 cnt++;
275 }
276 }
277 for (int j = 0; j < nEwgts; j++) {
278 adapter.setEdgeWeights(&ewgts[j*nnz], 1, j);
279 }
280 }
281 }
282
283
285 Zoltan2::PartitioningProblem<SparseGraphAdapter> problem(&adapter, &params);
286
287 try {
288 if (me == 0) std::cout << "Calling solve() " << std::endl;
289
290 problem.solve();
291
292 if (me == 0) std::cout << "Done solve() " << std::endl;
293 }
294 catch (std::runtime_error &e) {
295 delete [] vwgts;
296 delete [] ewgts;
297 std::cout << "Runtime exception returned from solve(): " << e.what();
298 if (!strncmp(e.what(), "BUILD ERROR", 11)) {
299 // Catching build errors as exceptions is OK in the tests
300 std::cout << " PASS" << std::endl;
301 return 0;
302 }
303 else {
304 // All other runtime_errors are failures
305 std::cout << " FAIL" << std::endl;
306 return -1;
307 }
308 }
309 catch (std::logic_error &e) {
310 delete [] vwgts;
311 delete [] ewgts;
312 std::cout << "Logic exception returned from solve(): " << e.what()
313 << " FAIL" << std::endl;
314 return -1;
315 }
316 catch (std::bad_alloc &e) {
317 delete [] vwgts;
318 delete [] ewgts;
319 std::cout << "Bad_alloc exception returned from solve(): " << e.what()
320 << " FAIL" << std::endl;
321 return -1;
322 }
323 catch (std::exception &e) {
324 delete [] vwgts;
325 delete [] ewgts;
326 std::cout << "Unknown exception returned from solve(). " << e.what()
327 << " FAIL" << std::endl;
328 return -1;
329 }
330
333 size_t checkNparts = comm->getSize();
334 if(nParts != -1) checkNparts = size_t(nParts);
335 size_t checkLength = origMatrix->getLocalNumRows();
336
337 const SparseGraphAdapter::part_t *checkParts = problem.getSolution().getPartListView();
338
339 // Check for load balance
340 size_t *countPerPart = new size_t[checkNparts];
341 size_t *globalCountPerPart = new size_t[checkNparts];
342 zscalar_t *wtPerPart = (nVwgts ? new zscalar_t[checkNparts*nVwgts] : NULL);
343 zscalar_t *globalWtPerPart = (nVwgts ? new zscalar_t[checkNparts*nVwgts] : NULL);
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.;
346
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++) {
352 if (j != NNZ_IDX)
353 wtPerPart[checkParts[i]*nVwgts+j] += vwgts[i*nVwgts+j];
354 else
355 wtPerPart[checkParts[i]*nVwgts+j] += origMatrix->getNumEntriesInLocalRow(i);
356 }
357 }
358
359 // Quotient algorithm should produce the same result for each local row
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;
366 }
367 }
368
369
370 Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, checkNparts,
371 countPerPart, globalCountPerPart);
372 Teuchos::reduceAll<int, zscalar_t>(*comm, Teuchos::REDUCE_SUM,
373 checkNparts*nVwgts,
374 wtPerPart, globalWtPerPart);
375
376 size_t min = std::numeric_limits<std::size_t>::max();
377 size_t max = 0;
378 size_t sum = 0;
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];
384 }
385
386 if (me == 0) {
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"
394 : " ")
395 << std::endl;
396 std::cout << "Imbalance: " << max / avg << std::endl;
397 if (nVwgts) {
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];
407 }
408 }
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;
416 }
417 }
418 }
419
420 delete [] countPerPart;
421 delete [] wtPerPart;
422 delete [] globalCountPerPart;
423 delete [] globalWtPerPart;
424 delete [] vwgts;
425 delete [] ewgts;
426
427
429 if (me == 0) std::cout << "Redistributing matrix..." << std::endl;
430 SparseMatrix *redistribMatrix;
431 SparseMatrixAdapter adapterMatrix(origMatrix);
432 adapterMatrix.applyPartitioningSolution(*origMatrix, redistribMatrix,
433 problem.getSolution());
434 if (redistribMatrix->getGlobalNumRows() < 40) {
435 Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,false));
436 redistribMatrix->describe(out, Teuchos::VERB_EXTREME);
437 }
438
439 if (me == 0) std::cout << "Redistributing vectors..." << std::endl;
440 Vector *redistribVector;
441// std::vector<const zscalar_t *> weights;
442// std::vector<int> weightStrides;
443 MultiVectorAdapter adapterVector(origVector); //, weights, weightStrides);
444 adapterVector.applyPartitioningSolution(*origVector, redistribVector,
445 problem.getSolution());
446
447 RCP<Vector> redistribProd;
448 redistribProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
449 redistribMatrix->getRangeMap());
450
451 // Test redistributing an integer vector with the same solution.
452 // This test is mostly to make sure compilation always works.
453 RCP<IntVector> origIntVec;
454 IntVector *redistIntVec;
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);
459
460 IntVectorAdapter int_vec_adapter(origIntVec);
461 int_vec_adapter.applyPartitioningSolution(*origIntVec, redistIntVec,
462 problem.getSolution());
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;
469 haveFailure = true;
470 }
471 else if (me == 0) std::cout << " OK" << std::endl;
472 delete redistIntVec;
473
476
477 if (me == 0) std::cout << "Matvec original..." << std::endl;
478 origMatrix->apply(*origVector, *origProd);
479 z2TestScalar origNorm = origProd->norm2();
480 if (me == 0)
481 std::cout << "Norm of Original matvec prod: " << origNorm << std::endl;
482
483 if (me == 0) std::cout << "Matvec redistributed..." << std::endl;
484 redistribMatrix->apply(*redistribVector, *redistribProd);
485 z2TestScalar redistribNorm = redistribProd->norm2();
486 if (me == 0)
487 std::cout << "Norm of Redistributed matvec prod: " << redistribNorm << std::endl;
488
489 if (redistribNorm > origNorm+epsilon || redistribNorm < origNorm-epsilon) {
490 testReturn = 1;
491 haveFailure = true;
492 }
493
494 delete redistribVector;
495 delete redistribMatrix;
496
497 if (me == 0) {
498 if (testReturn) {
499 std::cout << "Mat-Vec product changed; FAIL" << std::endl;
500 haveFailure = true;
501 }
502 if (!haveFailure)
503 std::cout << "PASS" << std::endl;
504 }
505
506 return testReturn;
507}
#define nParts
Defines the PartitioningProblem class.
common code used by tests
float zscalar_t
Tpetra::Map ::local_ordinal_type zlno_t
std::string testDataFilePath(".")
Tpetra::Map ::global_ordinal_type zgno_t
Defines XpetraCrsGraphAdapter class.
Defines the XpetraCrsMatrixAdapter class.
Defines the XpetraMultiVectorAdapter.
int main()
InputTraits< User >::part_t part_t
PartitioningProblem sets up partitioning problems for the user.
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
void solve(bool updateInputData=true)
Direct the problem to create a solution.
Provides access for Zoltan2 to Xpetra::CrsGraph data.
void setVertexWeightIsDegree(int idx)
Specify an index for which the vertex weight should be the degree of the vertex.
void setVertexWeights(const scalar_t *val, int stride, int idx)
Provide a pointer to vertex weights.
void setEdgeWeights(const scalar_t *val, int stride, int idx)
Provide a pointer to edge weights.
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
Definition: coloring1.cpp:79
Tpetra::Vector< z2TestScalar, z2TestLO, z2TestGO > Vector
Definition: coloring1.cpp:80
zgno_t z2TestGO
Definition: coloring1.cpp:76
zscalar_t z2TestScalar
Definition: coloring1.cpp:77
Tpetra::Vector< int, z2TestLO, z2TestGO > IntVector
#define epsilon
Vector::node_type Node
zlno_t z2TestLO
Tpetra::CrsGraph< z2TestLO, z2TestGO > SparseGraph
Zoltan2::XpetraMultiVectorAdapter< IntVector > IntVectorAdapter
#define NNZ_IDX
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
Zoltan2::XpetraMultiVectorAdapter< Vector > MultiVectorAdapter
Zoltan2::XpetraCrsMatrixAdapter< SparseMatrix > SparseMatrixAdapter
Tpetra::Vector< z2TestScalar, z2TestLO, z2TestGO > Vector
zgno_t z2TestGO
Tpetra::Vector< int, z2TestLO, z2TestGO > IntVector
zscalar_t z2TestScalar
Zoltan2::XpetraCrsGraphAdapter< SparseGraph > SparseGraphAdapter