Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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_Vector.hpp>
58 #include <MatrixMarket_Tpetra.hpp>
59 
60 using 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 
76 typedef zlno_t z2TestLO;
77 typedef zgno_t z2TestGO;
79 
80 typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO> SparseMatrix;
81 typedef Tpetra::CrsGraph<z2TestLO, z2TestGO> SparseGraph;
82 typedef Tpetra::Vector<z2TestScalar, z2TestLO, z2TestGO> Vector;
83 typedef Vector::node_type Node;
84 
88 
89 
90 // Integer vector
91 typedef Tpetra::Vector<int, z2TestLO, z2TestGO> IntVector;
93 
94 #define epsilon 0.00000001
95 #define NNZ_IDX 1
96 
98 int 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 }
zgno_t z2TestGO
Definition: coloring1.cpp:76
zlno_t z2TestLO
Definition: coloring1.cpp:75
#define nParts
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
Definition: coloring1.cpp:79
#define NNZ_IDX
Provides access for Zoltan2 to Xpetra::CrsGraph data.
Zoltan2::XpetraCrsGraphAdapter< SparseGraph > SparseGraphAdapter
int main(int narg, char **arg)
Definition: coloring1.cpp:199
common code used by tests
typename InputTraits< User >::part_t part_t
Tpetra::Vector< z2TestScalar, z2TestLO, z2TestGO > Vector
Definition: coloring1.cpp:80
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
Definition: coloring1.cpp:85
Tpetra::Vector< int, z2TestLO, z2TestGO > IntVector
Defines the PartitioningProblem class.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
zscalar_t z2TestScalar
Definition: coloring1.cpp:77
float zscalar_t
#define epsilon
Zoltan2::XpetraMultiVectorAdapter< Vector > MultiVectorAdapter
Tpetra::Map::global_ordinal_type zgno_t
void solve(bool updateInputData=true)
Direct the problem to create a solution.
Vector::node_type Node
Definition: coloring1.cpp:81
std::string testDataFilePath(".")