partitioning1.cpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//   Zoltan2: A package of combinatorial algorithms for scientific computing
//
// Copyright 2012 NTESS and the Zoltan2 contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#include <Zoltan2_PartitioningProblem.hpp>
#include <Zoltan2_XpetraCrsMatrixAdapter.hpp>
#include <Zoltan2_XpetraCrsGraphAdapter.hpp>
#include <Zoltan2_XpetraMultiVectorAdapter.hpp>
#include <Zoltan2_TestHelpers.hpp>
#include <iostream>
#include <limits>
#include <Teuchos_ParameterList.hpp>
#include <Teuchos_RCP.hpp>
#include <Teuchos_FancyOStream.hpp>
#include <Teuchos_CommandLineProcessor.hpp>
#include <Tpetra_CrsMatrix.hpp>
#include <Tpetra_Vector.hpp>
#include <MatrixMarket_Tpetra.hpp>

using Teuchos::RCP;

// Program to demonstrate use of Zoltan2 to partition a TPetra matrix
// (read from a MatrixMarket file or generated by Galeri::Xpetra).
// Usage:
//     a.out [--inputFile=filename] [--outputFile=outfile] [--verbose]
//           [--x=#] [--y=#] [--z=#] [--matrix={Laplace1D,Laplace2D,Laplace3D}
// Karen Devine, 2011

// Eventually want to use Teuchos unit tests to vary z2TestLO and
// GO.  For now, we set them at compile time based on whether Tpetra
// is built with explicit instantiation on.  (in Zoltan2_TestHelpers.hpp)

typedef zlno_t z2TestLO;
typedef zgno_t z2TestGO;
typedef zscalar_t z2TestScalar;

typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO> SparseMatrix;
typedef Tpetra::CrsGraph<z2TestLO, z2TestGO> SparseGraph;
typedef Tpetra::Vector<z2TestScalar, z2TestLO, z2TestGO> Vector;
typedef Vector::node_type Node;

typedef Zoltan2::XpetraCrsMatrixAdapter<SparseMatrix> SparseMatrixAdapter;
typedef Zoltan2::XpetraCrsGraphAdapter<SparseGraph> SparseGraphAdapter;
typedef Zoltan2::XpetraMultiVectorAdapter<Vector> MultiVectorAdapter;


// Integer vector
typedef Tpetra::Vector<int, z2TestLO, z2TestGO> IntVector;
typedef Zoltan2::XpetraMultiVectorAdapter<IntVector> IntVectorAdapter;

#define epsilon 0.00000001
#define NNZ_IDX 1

int main(int narg, char** arg)
{
  std::string inputFile = "";        // Matrix Market or Zoltan file to read
  std::string outputFile = "";       // Matrix Market or Zoltan file to write
  std::string inputPath = testDataFilePath;  // Directory with input file
  std::string method = "scotch";
  bool verbose = false;              // Verbosity of output
  bool distributeInput = true;
  bool haveFailure = false;
  int nParts = -1;
  int nVwgts = 0;
  int nEwgts = 0;
  int testReturn = 0;

  Tpetra::ScopeGuard tscope(&narg, &arg);
  RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
  int me = comm->getRank();

  // Read run-time options.
  Teuchos::CommandLineProcessor cmdp (false, false);
  cmdp.setOption("inputPath", &inputPath,
                 "Path to the MatrixMarket or Zoltan file to be read; "
                 "if not specified, a default path will be used.");
  cmdp.setOption("inputFile", &inputFile,
                 "Name of the Matrix Market or Zoltan file to read; "
                 "if not specified, a matrix will be generated by MueLu.");
  cmdp.setOption("outputFile", &outputFile,
                 "Name of the Matrix Market sparse matrix file to write, "
                 "echoing the input/generated matrix.");
  cmdp.setOption("method", &method,
                 "Partitioning method to use:  scotch or parmetis.");
  cmdp.setOption("nparts", &nParts,
                 "Number of parts being requested");
  cmdp.setOption("vertexWeights", &nVwgts,
                 "Number of weights to generate for each vertex");
  cmdp.setOption("edgeWeights", &nEwgts,
                 "Number of weights to generate for each edge");
  cmdp.setOption("verbose", "quiet", &verbose,
                 "Print messages and results.");
  cmdp.setOption("distribute", "no-distribute", &distributeInput,
                "indicate whether or not to distribute "
                "input across the communicator");

  // Even with cmdp option "true", I get errors for having these
  //   arguments on the command line.  (On redsky build)
  // KDDKDD Should just be warnings, right?  Code should still work with these
  // KDDKDD params in the create-a-matrix file.  Better to have them where
  // KDDKDD they are used.
  int xdim=10;
  int ydim=10;
  int zdim=10;
  std::string matrixType("Laplace3D");

  cmdp.setOption("x", &xdim,
                "number of gridpoints in X dimension for "
                "mesh used to generate matrix.");
  cmdp.setOption("y", &ydim,
                "number of gridpoints in Y dimension for "
                "mesh used to generate matrix.");
  cmdp.setOption("z", &zdim,
                "number of gridpoints in Z dimension for "
                "mesh used to generate matrix.");
  cmdp.setOption("matrix", &matrixType,
                "Matrix type: Laplace1D, Laplace2D, or Laplace3D");

  // Quotient-specific parameters
  int quotientThreshold = -1;
  cmdp.setOption("qthreshold", &quotientThreshold,
                "Threshold on the number of vertices for active MPI ranks to hold"
    "after the migrating the communication graph to the active ranks.");  


  cmdp.parse(narg, arg);

  RCP<UserInputForTests> uinput;

  if (inputFile != "")   // Input file specified; read a matrix
    uinput = rcp(new UserInputForTests(inputPath, inputFile, comm,
                                       true, distributeInput));

  else                  // Let MueLu generate a default matrix
    uinput = rcp(new UserInputForTests(xdim, ydim, zdim, string(""), comm,
                                       true, distributeInput));

  RCP<SparseMatrix> origMatrix = uinput->getUITpetraCrsMatrix();

  if (origMatrix->getGlobalNumRows() < 40) {
    Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,false));
    origMatrix->describe(out, Teuchos::VERB_EXTREME);
  }


  if (outputFile != "") {
    // Just a sanity check.
    Tpetra::MatrixMarket::Writer<SparseMatrix>::writeSparseFile(outputFile,
                                                origMatrix, verbose);
  }

  if (me == 0)
    std::cout << "NumRows     = " << origMatrix->getGlobalNumRows() << std::endl
         << "NumNonzeros = " << origMatrix->getGlobalNumEntries() << std::endl
         << "NumProcs = " << comm->getSize() << std::endl
         << "NumLocalRows (rank 0) = " << origMatrix->getLocalNumRows() << std::endl;

  RCP<Vector> origVector, origProd;
  origProd   = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
                                    origMatrix->getRangeMap());
  origVector = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
                                    origMatrix->getDomainMap());
  origVector->randomize();

  Teuchos::ParameterList params;

  params.set("partitioning_approach", "partition");
  params.set("algorithm", method);

  if(nParts > 0) {
    params.set("num_global_parts", nParts);
  }

  if(method == "quotient" && quotientThreshold > 0) {
    params.set("quotient_threshold", quotientThreshold);    
  }

  SparseGraphAdapter adapter(origMatrix->getCrsGraph(), nVwgts, nEwgts);


  zscalar_t *vwgts = NULL, *ewgts = NULL;
  if (nVwgts) {
    // Test vertex weights with stride nVwgts.
    size_t nrows = origMatrix->getLocalNumRows();
    if (nrows) {
      vwgts = new zscalar_t[nVwgts * nrows];
      for (size_t i = 0; i < nrows; i++) {
        size_t idx = i * nVwgts;
        vwgts[idx] = zscalar_t(origMatrix->getRowMap()->getGlobalElement(i))
  ;//                 + zscalar_t(0.5);
        for (int j = 1; j < nVwgts; j++) vwgts[idx+j] = 1.;
      }
      for (int j = 0; j < nVwgts; j++) {
        if (j != NNZ_IDX) adapter.setVertexWeights(&vwgts[j], nVwgts, j);
        else              adapter.setVertexWeightIsDegree(NNZ_IDX);
      }
    }
  }

  if (nEwgts) {
    // Test edge weights with stride 1.
    size_t nnz = origMatrix->getLocalNumEntries();
    if (nnz) {
      size_t nrows = origMatrix->getLocalNumRows();
      size_t maxnzrow = origMatrix->getLocalMaxNumRowEntries();
      ewgts = new zscalar_t[nEwgts * nnz];
      size_t cnt = 0;
      typename SparseMatrix::nonconst_global_inds_host_view_type  egids("egids", maxnzrow);
      typename SparseMatrix::nonconst_values_host_view_type evals("evals", maxnzrow);
      for (size_t i = 0; i < nrows; i++) {
        size_t nnzinrow;
        z2TestGO gid = origMatrix->getRowMap()->getGlobalElement(i);
        origMatrix->getGlobalRowCopy(gid, egids, evals, nnzinrow);
        for (size_t k = 0; k < nnzinrow; k++) {
          ewgts[cnt] = (gid < egids[k] ? gid : egids[k]);
          if (nEwgts > 1) ewgts[cnt+nnz] = (gid < egids[k] ? egids[k] : gid);
          for (int j = 2; j < nEwgts; j++) ewgts[cnt+nnz*j] = 1.;
          cnt++;
        }
      }
      for (int j = 0; j < nEwgts; j++) {
        adapter.setEdgeWeights(&ewgts[j*nnz], 1, j);
      }
    }
  }


  Zoltan2::PartitioningProblem<SparseGraphAdapter> problem(&adapter, &params);

  try {
    if (me == 0) std::cout << "Calling solve() " << std::endl;

    problem.solve();

    if (me == 0) std::cout << "Done solve() " << std::endl;
  }
  catch (std::runtime_error &e) {
    delete [] vwgts;
    delete [] ewgts;
    std::cout << "Runtime exception returned from solve(): " << e.what();
    if (!strncmp(e.what(), "BUILD ERROR", 11)) {
      // Catching build errors as exceptions is OK in the tests
      std::cout << " PASS" << std::endl;
      return 0;
    }
    else {
      // All other runtime_errors are failures
      std::cout << " FAIL" << std::endl;
      return -1;
    }
  }
  catch (std::logic_error &e) {
    delete [] vwgts;
    delete [] ewgts;
    std::cout << "Logic exception returned from solve(): " << e.what()
         << " FAIL" << std::endl;
    return -1;
  }
  catch (std::bad_alloc &e) {
    delete [] vwgts;
    delete [] ewgts;
    std::cout << "Bad_alloc exception returned from solve(): " << e.what()
         << " FAIL" << std::endl;
    return -1;
  }
  catch (std::exception &e) {
    delete [] vwgts;
    delete [] ewgts;
    std::cout << "Unknown exception returned from solve(). " << e.what()
         << " FAIL" << std::endl;
    return -1;
  }
  
  size_t checkNparts = comm->getSize();
  if(nParts != -1) checkNparts = size_t(nParts);
  size_t checkLength = origMatrix->getLocalNumRows();

  const SparseGraphAdapter::part_t *checkParts = problem.getSolution().getPartListView();

  // Check for load balance
  size_t *countPerPart = new size_t[checkNparts];
  size_t *globalCountPerPart = new size_t[checkNparts];
  zscalar_t *wtPerPart = (nVwgts ? new zscalar_t[checkNparts*nVwgts] : NULL);
  zscalar_t *globalWtPerPart = (nVwgts ? new zscalar_t[checkNparts*nVwgts] : NULL);
  for (size_t i = 0; i < checkNparts; i++) countPerPart[i] = 0;
  for (size_t i = 0; i < checkNparts * nVwgts; i++) wtPerPart[i] = 0.;

  for (size_t i = 0; i < checkLength; i++) {
    if (size_t(checkParts[i]) >= checkNparts)
      std::cout << "Invalid Part " << checkParts[i] << ": FAIL" << std::endl;
    countPerPart[checkParts[i]]++;
    for (int j = 0; j < nVwgts; j++) {
      if (j != NNZ_IDX)
        wtPerPart[checkParts[i]*nVwgts+j] += vwgts[i*nVwgts+j];
      else
        wtPerPart[checkParts[i]*nVwgts+j] += origMatrix->getNumEntriesInLocalRow(i);
    }
  }

  // Quotient algorithm should produce the same result for each local row
  if(method == "quotient") {
    size_t result = size_t(checkParts[0]);    
    for (size_t i = 1; i < checkLength; i++) {
      if (size_t(checkParts[i]) != result)
  std::cout << "Different parts in the quotient algorithm: " 
      << result << "!=" << checkParts[i] << ": FAIL" << std::endl;
    }
  }


  Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, checkNparts,
                                  countPerPart, globalCountPerPart);
  Teuchos::reduceAll<int, zscalar_t>(*comm, Teuchos::REDUCE_SUM,
                                    checkNparts*nVwgts,
                                    wtPerPart, globalWtPerPart);

  size_t min = std::numeric_limits<std::size_t>::max();
  size_t max = 0;
  size_t sum = 0;
  size_t minrank = 0, maxrank = 0;
  for (size_t i = 0; i < checkNparts; i++) {
    if (globalCountPerPart[i] < min) {min = globalCountPerPart[i]; minrank = i;}
    if (globalCountPerPart[i] > max) {max = globalCountPerPart[i]; maxrank = i;}
    sum += globalCountPerPart[i];
  }
  
  if (me == 0) {
    float avg = (float) sum / (float) checkNparts;
    std::cout << "Minimum count:  " << min << " on rank " << minrank << std::endl;
    std::cout << "Maximum count:  " << max << " on rank " << maxrank << std::endl;
    std::cout << "Average count:  " << avg << std::endl;
    std::cout << "Total count:    " << sum
         << (sum != origMatrix->getGlobalNumRows()
                 ? "Work was lost; FAIL"
                 : " ")
         << std::endl;
    std::cout << "Imbalance:     " << max / avg << std::endl;
    if (nVwgts) {
      std::vector<zscalar_t> minwt(nVwgts, std::numeric_limits<zscalar_t>::max());
      std::vector<zscalar_t> maxwt(nVwgts, 0.);
      std::vector<zscalar_t> sumwt(nVwgts, 0.);
      for (size_t i = 0; i < checkNparts; i++) {
        for (int j = 0; j < nVwgts; j++) {
          size_t idx = i*nVwgts+j;
          if (globalWtPerPart[idx] < minwt[j]) minwt[j] = globalWtPerPart[idx];
          if (globalWtPerPart[idx] > maxwt[j]) maxwt[j] = globalWtPerPart[idx];
          sumwt[j] += globalWtPerPart[idx];
        }
      }
      for (int j = 0; j < nVwgts; j++) {
        float avgwt = (float) sumwt[j] / (float) checkNparts;
        std::cout << std::endl;
        std::cout << "Minimum weight[" << j << "]:  " << minwt[j] << std::endl;
        std::cout << "Maximum weight[" << j << "]:  " << maxwt[j] << std::endl;
        std::cout << "Average weight[" << j << "]:  " << avgwt << std::endl;
        std::cout << "Imbalance:       " << maxwt[j] / avgwt << std::endl;
      }
    }
  }

  delete [] countPerPart;
  delete [] wtPerPart;
  delete [] globalCountPerPart;
  delete [] globalWtPerPart;
  delete [] vwgts;
  delete [] ewgts;


  if (me == 0) std::cout << "Redistributing matrix..." << std::endl;
  SparseMatrix *redistribMatrix;
  SparseMatrixAdapter adapterMatrix(origMatrix);
  adapterMatrix.applyPartitioningSolution(*origMatrix, redistribMatrix,
                                          problem.getSolution());
  if (redistribMatrix->getGlobalNumRows() < 40) {
    Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,false));
    redistribMatrix->describe(out, Teuchos::VERB_EXTREME);
  }

  if (me == 0) std::cout << "Redistributing vectors..." << std::endl;
  Vector *redistribVector;
//  std::vector<const zscalar_t *> weights;
//  std::vector<int> weightStrides;
  MultiVectorAdapter adapterVector(origVector); //, weights, weightStrides);
  adapterVector.applyPartitioningSolution(*origVector, redistribVector,
                                          problem.getSolution());

  RCP<Vector> redistribProd;
  redistribProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
                                       redistribMatrix->getRangeMap());

  // Test redistributing an integer vector with the same solution.
  // This test is mostly to make sure compilation always works.
  RCP<IntVector> origIntVec;
  IntVector *redistIntVec;
  origIntVec = Tpetra::createVector<int,z2TestLO,z2TestGO>(
                                        origMatrix->getRangeMap());
  for (size_t i = 0; i < origIntVec->getLocalLength(); i++)
    origIntVec->replaceLocalValue(i, me);

  IntVectorAdapter int_vec_adapter(origIntVec);
  int_vec_adapter.applyPartitioningSolution(*origIntVec, redistIntVec,
                                             problem.getSolution());
  int origIntNorm = origIntVec->norm1();
  int redistIntNorm = redistIntVec->norm1();
  if (me == 0) std::cout << "IntegerVectorTest:  " << origIntNorm << " == "
                         << redistIntNorm << " ?";
  if (origIntNorm != redistIntNorm) {
    if (me == 0) std::cout << " FAIL" << std::endl;
    haveFailure = true;
  }
  else if (me == 0) std::cout << " OK" << std::endl;
  delete redistIntVec;


  if (me == 0) std::cout << "Matvec original..." << std::endl;
  origMatrix->apply(*origVector, *origProd);
  z2TestScalar origNorm = origProd->norm2();
  if (me == 0)
    std::cout << "Norm of Original matvec prod:       " << origNorm << std::endl;

  if (me == 0) std::cout << "Matvec redistributed..." << std::endl;
  redistribMatrix->apply(*redistribVector, *redistribProd);
  z2TestScalar redistribNorm = redistribProd->norm2();
  if (me == 0)
    std::cout << "Norm of Redistributed matvec prod:  " << redistribNorm << std::endl;

  if (redistribNorm > origNorm+epsilon || redistribNorm < origNorm-epsilon) {
    testReturn = 1;
    haveFailure = true;
  }

  delete redistribVector;
  delete redistribMatrix;

  if (me == 0) {
    if (testReturn) {
      std::cout << "Mat-Vec product changed; FAIL" << std::endl;
      haveFailure = true;
    }
    if (!haveFailure)
      std::cout << "PASS" << std::endl;
  }
  
  return testReturn;
}