Zoltan2
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
49 #include <Zoltan2_TestHelpers.hpp>
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_DefaultPlatform.hpp>
58 #include <Tpetra_Vector.hpp>
59 #include <MatrixMarket_Tpetra.hpp>
60 
61 using Teuchos::RCP;
62 using namespace std;
63 
65 // Program to demonstrate use of Zoltan2 to partition a TPetra matrix
66 // (read from a MatrixMarket file or generated by Galeri::Xpetra).
67 // Usage:
68 // a.out [--inputFile=filename] [--outputFile=outfile] [--verbose]
69 // [--x=#] [--y=#] [--z=#] [--matrix={Laplace1D,Laplace2D,Laplace3D}
70 // Karen Devine, 2011
72 
74 // Eventually want to use Teuchos unit tests to vary z2TestLO and
75 // GO. For now, we set them at compile time based on whether Tpetra
76 // is built with explicit instantiation on. (in Zoltan2_TestHelpers.hpp)
77 
78 typedef zlno_t z2TestLO;
79 typedef zgno_t z2TestGO;
81 
82 typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO> SparseMatrix;
83 typedef Tpetra::CrsGraph<z2TestLO, z2TestGO> SparseGraph;
84 typedef Tpetra::Vector<z2TestScalar, z2TestLO, z2TestGO> Vector;
85 typedef Vector::node_type Node;
86 
90 
91 
92 // Integer vector
93 typedef Tpetra::Vector<int, z2TestLO, z2TestGO> IntVector;
95 
96 #define epsilon 0.00000001
97 #define NNZ_IDX 1
98 
100 int main(int narg, char** arg)
101 {
102  std::string inputFile = ""; // Matrix Market or Zoltan file to read
103  std::string outputFile = ""; // Matrix Market or Zoltan file to write
104  std::string inputPath = testDataFilePath; // Directory with input file
105  std::string method = "scotch";
106  bool verbose = false; // Verbosity of output
107  bool distributeInput = true;
108  bool haveFailure = false;
109  int nVwgts = 0;
110  int nEwgts = 0;
111  int testReturn = 0;
112 
114  Teuchos::GlobalMPISession mpiSession(&narg, &arg, NULL);
115  RCP<const Teuchos::Comm<int> > comm =
116  Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
117  int me = comm->getRank();
118 
119  // Read run-time options.
120  Teuchos::CommandLineProcessor cmdp (false, false);
121  cmdp.setOption("inputPath", &inputPath,
122  "Path to the MatrixMarket or Zoltan file to be read; "
123  "if not specified, a default path will be used.");
124  cmdp.setOption("inputFile", &inputFile,
125  "Name of the Matrix Market or Zoltan file to read; "
126  "if not specified, a matrix will be generated by MueLu.");
127  cmdp.setOption("outputFile", &outputFile,
128  "Name of the Matrix Market sparse matrix file to write, "
129  "echoing the input/generated matrix.");
130  cmdp.setOption("method", &method,
131  "Partitioning method to use: scotch or parmetis.");
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");
165 
166  cmdp.parse(narg, arg);
167 
168  RCP<UserInputForTests> uinput;
169 
170  if (inputFile != "") // Input file specified; read a matrix
171  uinput = rcp(new UserInputForTests(inputPath, inputFile, comm,
172  true, distributeInput));
173 
174  else // Let MueLu generate a default matrix
175  uinput = rcp(new UserInputForTests(xdim, ydim, zdim, string(""), comm,
176  true, distributeInput));
177 
178  RCP<SparseMatrix> origMatrix = uinput->getUITpetraCrsMatrix();
179 
180  if (origMatrix->getGlobalNumRows() < 40) {
181  Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,false));
182  origMatrix->describe(out, Teuchos::VERB_EXTREME);
183  }
184 
185 
186  if (outputFile != "") {
187  // Just a sanity check.
188  Tpetra::MatrixMarket::Writer<SparseMatrix>::writeSparseFile(outputFile,
189  origMatrix, verbose);
190  }
191 
192  if (me == 0)
193  cout << "NumRows = " << origMatrix->getGlobalNumRows() << endl
194  << "NumNonzeros = " << origMatrix->getGlobalNumEntries() << endl
195  << "NumProcs = " << comm->getSize() << endl
196  << "NumLocalRows (rank 0) = " << origMatrix->getNodeNumRows() << endl;
197 
199  RCP<Vector> origVector, origProd;
200  origProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
201  origMatrix->getRangeMap());
202  origVector = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
203  origMatrix->getDomainMap());
204  origVector->randomize();
205 
207  Teuchos::ParameterList params;
208 
209  params.set("partitioning_approach", "partition");
210  params.set("algorithm", method);
211 
213  SparseGraphAdapter adapter(origMatrix->getCrsGraph(), nVwgts, nEwgts);
214 
218 
219  zscalar_t *vwgts = NULL, *ewgts = NULL;
220  if (nVwgts) {
221  // Test vertex weights with stride nVwgts.
222  size_t nrows = origMatrix->getNodeNumRows();
223  if (nrows) {
224  vwgts = new zscalar_t[nVwgts * nrows];
225  for (size_t i = 0; i < nrows; i++) {
226  size_t idx = i * nVwgts;
227  vwgts[idx] = zscalar_t(origMatrix->getRowMap()->getGlobalElement(i))
228  ;// + zscalar_t(0.5);
229  for (int j = 1; j < nVwgts; j++) vwgts[idx+j] = 1.;
230  }
231  for (int j = 0; j < nVwgts; j++) {
232  if (j != NNZ_IDX) adapter.setVertexWeights(&vwgts[j], nVwgts, j);
233  else adapter.setVertexWeightIsDegree(NNZ_IDX);
234  }
235  }
236  }
237 
238  if (nEwgts) {
239  // Test edge weights with stride 1.
240  size_t nnz = origMatrix->getNodeNumEntries();
241  if (nnz) {
242  size_t nrows = origMatrix->getNodeNumRows();
243  size_t maxnzrow = origMatrix->getNodeMaxNumRowEntries();
244  ewgts = new zscalar_t[nEwgts * nnz];
245  size_t cnt = 0;
246  Array<z2TestGO> egids(maxnzrow);
247  Array<zscalar_t> evals(maxnzrow);
248  for (size_t i = 0; i < nrows; i++) {
249  size_t nnzinrow;
250  z2TestGO gid = origMatrix->getRowMap()->getGlobalElement(i);
251  origMatrix->getGlobalRowCopy(gid, egids(), evals(), nnzinrow);
252  for (size_t k = 0; k < nnzinrow; k++) {
253  ewgts[cnt] = (gid < egids[k] ? gid : egids[k]);
254  if (nEwgts > 1) ewgts[cnt+nnz] = (gid < egids[k] ? egids[k] : gid);
255  for (int j = 2; j < nEwgts; j++) ewgts[cnt+nnz*j] = 1.;
256  cnt++;
257  }
258  }
259  for (int j = 0; j < nEwgts; j++) {
260  adapter.setEdgeWeights(&ewgts[j*nnz], 1, j);
261  }
262  }
263  }
264 
265 
267  Zoltan2::PartitioningProblem<SparseGraphAdapter> problem(&adapter, &params);
268 
269  try {
270  if (me == 0) cout << "Calling solve() " << endl;
271 
272  problem.solve();
273 
274  if (me == 0) cout << "Done solve() " << endl;
275  }
276  catch (std::runtime_error &e) {
277  delete [] vwgts;
278  delete [] ewgts;
279  cout << "Runtime exception returned from solve(): " << e.what();
280  if (!strncmp(e.what(), "BUILD ERROR", 11)) {
281  // Catching build errors as exceptions is OK in the tests
282  cout << " PASS" << endl;
283  return 0;
284  }
285  else {
286  // All other runtime_errors are failures
287  cout << " FAIL" << endl;
288  return -1;
289  }
290  }
291  catch (std::logic_error &e) {
292  delete [] vwgts;
293  delete [] ewgts;
294  cout << "Logic exception returned from solve(): " << e.what()
295  << " FAIL" << endl;
296  return -1;
297  }
298  catch (std::bad_alloc &e) {
299  delete [] vwgts;
300  delete [] ewgts;
301  cout << "Bad_alloc exception returned from solve(): " << e.what()
302  << " FAIL" << endl;
303  return -1;
304  }
305  catch (std::exception &e) {
306  delete [] vwgts;
307  delete [] ewgts;
308  cout << "Unknown exception returned from solve(). " << e.what()
309  << " FAIL" << endl;
310  return -1;
311  }
312 
315  size_t checkNparts = comm->getSize();
316 
317  size_t checkLength = origMatrix->getNodeNumRows();
318  const SparseGraphAdapter::part_t *checkParts = problem.getSolution().getPartListView();
319 
320  // Check for load balance
321  size_t *countPerPart = new size_t[checkNparts];
322  size_t *globalCountPerPart = new size_t[checkNparts];
323  zscalar_t *wtPerPart = (nVwgts ? new zscalar_t[checkNparts*nVwgts] : NULL);
324  zscalar_t *globalWtPerPart = (nVwgts ? new zscalar_t[checkNparts*nVwgts] : NULL);
325  for (size_t i = 0; i < checkNparts; i++) countPerPart[i] = 0;
326  for (size_t i = 0; i < checkNparts * nVwgts; i++) wtPerPart[i] = 0.;
327 
328  for (size_t i = 0; i < checkLength; i++) {
329  if (size_t(checkParts[i]) >= checkNparts)
330  cout << "Invalid Part " << checkParts[i] << ": FAIL" << endl;
331  countPerPart[checkParts[i]]++;
332  for (int j = 0; j < nVwgts; j++) {
333  if (j != NNZ_IDX)
334  wtPerPart[checkParts[i]*nVwgts+j] += vwgts[i*nVwgts+j];
335  else
336  wtPerPart[checkParts[i]*nVwgts+j] += origMatrix->getNumEntriesInLocalRow(i);
337  }
338  }
339  Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, checkNparts,
340  countPerPart, globalCountPerPart);
341  Teuchos::reduceAll<int, zscalar_t>(*comm, Teuchos::REDUCE_SUM,
342  checkNparts*nVwgts,
343  wtPerPart, globalWtPerPart);
344 
345  size_t min = std::numeric_limits<std::size_t>::max();
346  size_t max = 0;
347  size_t sum = 0;
348  size_t minrank = 0, maxrank = 0;
349  for (size_t i = 0; i < checkNparts; i++) {
350  if (globalCountPerPart[i] < min) {min = globalCountPerPart[i]; minrank = i;}
351  if (globalCountPerPart[i] > max) {max = globalCountPerPart[i]; maxrank = i;}
352  sum += globalCountPerPart[i];
353  }
354 
355  if (me == 0) {
356  float avg = (float) sum / (float) checkNparts;
357  cout << "Minimum count: " << min << " on rank " << minrank << endl;
358  cout << "Maximum count: " << max << " on rank " << maxrank << endl;
359  cout << "Average count: " << avg << endl;
360  cout << "Total count: " << sum
361  << (sum != origMatrix->getGlobalNumRows()
362  ? "Work was lost; FAIL"
363  : " ")
364  << endl;
365  cout << "Imbalance: " << max / avg << endl;
366  if (nVwgts) {
367  std::vector<zscalar_t> minwt(nVwgts, std::numeric_limits<zscalar_t>::max());
368  std::vector<zscalar_t> maxwt(nVwgts, 0.);
369  std::vector<zscalar_t> sumwt(nVwgts, 0.);
370  for (size_t i = 0; i < checkNparts; i++) {
371  for (int j = 0; j < nVwgts; j++) {
372  size_t idx = i*nVwgts+j;
373  if (globalWtPerPart[idx] < minwt[j]) minwt[j] = globalWtPerPart[idx];
374  if (globalWtPerPart[idx] > maxwt[j]) maxwt[j] = globalWtPerPart[idx];
375  sumwt[j] += globalWtPerPart[idx];
376  }
377  }
378  for (int j = 0; j < nVwgts; j++) {
379  float avgwt = (float) sumwt[j] / (float) checkNparts;
380  cout << endl;
381  cout << "Minimum weight[" << j << "]: " << minwt[j] << endl;
382  cout << "Maximum weight[" << j << "]: " << maxwt[j] << endl;
383  cout << "Average weight[" << j << "]: " << avgwt << endl;
384  cout << "Imbalance: " << maxwt[j] / avgwt << endl;
385  }
386  }
387  }
388 
389  delete [] countPerPart;
390  delete [] wtPerPart;
391  delete [] globalCountPerPart;
392  delete [] globalWtPerPart;
393  delete [] vwgts;
394  delete [] ewgts;
395 
396 
398  if (me == 0) cout << "Redistributing matrix..." << endl;
399  SparseMatrix *redistribMatrix;
400  SparseMatrixAdapter adapterMatrix(origMatrix);
401  adapterMatrix.applyPartitioningSolution(*origMatrix, redistribMatrix,
402  problem.getSolution());
403  if (redistribMatrix->getGlobalNumRows() < 40) {
404  Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,false));
405  redistribMatrix->describe(out, Teuchos::VERB_EXTREME);
406  }
407 
408  if (me == 0) cout << "Redistributing vectors..." << endl;
409  Vector *redistribVector;
410 // std::vector<const zscalar_t *> weights;
411 // std::vector<int> weightStrides;
412  MultiVectorAdapter adapterVector(origVector); //, weights, weightStrides);
413  adapterVector.applyPartitioningSolution(*origVector, redistribVector,
414  problem.getSolution());
415 
416  RCP<Vector> redistribProd;
417  redistribProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
418  redistribMatrix->getRangeMap());
419 
420  // Test redistributing an integer vector with the same solution.
421  // This test is mostly to make sure compilation always works.
422  RCP<IntVector> origIntVec;
423  IntVector *redistIntVec;
424  origIntVec = Tpetra::createVector<int,z2TestLO,z2TestGO>(
425  origMatrix->getRangeMap());
426  for (size_t i = 0; i < origIntVec->getLocalLength(); i++)
427  origIntVec->replaceLocalValue(i, me);
428 
429  IntVectorAdapter int_vec_adapter(origIntVec);
430  int_vec_adapter.applyPartitioningSolution(*origIntVec, redistIntVec,
431  problem.getSolution());
432  int origIntNorm = origIntVec->norm1();
433  int redistIntNorm = redistIntVec->norm1();
434  if (me == 0) std::cout << "IntegerVectorTest: " << origIntNorm << " == "
435  << redistIntNorm << " ?";
436  if (origIntNorm != redistIntNorm) {
437  if (me == 0) std::cout << " FAIL" << std::endl;
438  haveFailure = true;
439  }
440  else if (me == 0) std::cout << " OK" << std::endl;
441  delete redistIntVec;
442 
445 
446  if (me == 0) cout << "Matvec original..." << endl;
447  origMatrix->apply(*origVector, *origProd);
448  z2TestScalar origNorm = origProd->norm2();
449  if (me == 0)
450  cout << "Norm of Original matvec prod: " << origNorm << endl;
451 
452  if (me == 0) cout << "Matvec redistributed..." << endl;
453  redistribMatrix->apply(*redistribVector, *redistribProd);
454  z2TestScalar redistribNorm = redistribProd->norm2();
455  if (me == 0)
456  cout << "Norm of Redistributed matvec prod: " << redistribNorm << endl;
457 
458  if (redistribNorm > origNorm+epsilon || redistribNorm < origNorm-epsilon) {
459  testReturn = 1;
460  haveFailure = true;
461  }
462 
463  delete redistribVector;
464  delete redistribMatrix;
465 
466  if (me == 0) {
467  if (testReturn) {
468  std::cout << "Mat-Vec product changed; FAIL" << std::endl;
469  haveFailure = true;
470  }
471  if (!haveFailure)
472  std::cout << "PASS" << std::endl;
473  }
474 
475  return testReturn;
476 }
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
double zscalar_t
zgno_t z2TestGO
#define NNZ_IDX
Provides access for Zoltan2 to Xpetra::CrsGraph data.
Zoltan2::XpetraCrsGraphAdapter< SparseGraph > SparseGraphAdapter
int zlno_t
zlno_t z2TestLO
common code used by tests
Tpetra::CrsGraph< z2TestLO, z2TestGO > SparseGraph
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
Zoltan2::XpetraMultiVectorAdapter< IntVector > IntVectorAdapter
Defines the XpetraMultiVectorAdapter.
Defines XpetraCrsGraphAdapter class.
Defines the XpetraCrsMatrixAdapter class.
Tpetra::Vector< int, z2TestLO, z2TestGO > IntVector
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
InputTraits< User >::part_t part_t
An adapter for Xpetra::MultiVector.
int zgno_t
Tpetra::Vector< z2TestScalar, z2TestLO, z2TestGO > Vector
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
PartitioningProblem sets up partitioning problems for the user.
zscalar_t z2TestScalar
Zoltan2::XpetraMultiVectorAdapter< Vector > MultiVectorAdapter
Zoltan2::XpetraCrsMatrixAdapter< SparseMatrix > SparseMatrixAdapter
Vector::node_type Node
Defines the PartitioningProblem class.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
int main(int narg, char **arg)
#define epsilon
void solve(bool updateInputData=true)
Direct the problem to create a solution.
std::string testDataFilePath(".")