Zoltan2
graph.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
45 #include <iostream>
46 #include <limits>
50 #include <Teuchos_ParameterList.hpp>
51 #include <Teuchos_RCP.hpp>
52 #include <Teuchos_CommandLineProcessor.hpp>
53 #include <Tpetra_DefaultPlatform.hpp>
54 #include <Tpetra_CrsMatrix.hpp>
55 #include <Tpetra_Vector.hpp>
56 #include <Galeri_XpetraMaps.hpp>
57 #include <Galeri_XpetraProblemFactory.hpp>
58 
59 using Teuchos::RCP;
60 using namespace std;
61 
63 // Program to demonstrate use of Zoltan2 to partition a TPetra matrix
64 // using graph partitioning via Scotch or ParMETIS.
66 
67 int main(int narg, char** arg)
68 {
69  // Establish session; works both for MPI and non-MPI builds
70  Teuchos::GlobalMPISession mpiSession(&narg, &arg, NULL);
71  RCP<const Teuchos::Comm<int> > comm =
72  Tpetra::DefaultPlatform::getDefaultPlatform().getComm();
73  int me = comm->getRank();
74 
75  // Useful typedefs: Tpetra types
76  // In this example, we'll use Tpetra defaults for local/global ID type
77  typedef double scalar_t;
78  typedef Tpetra::Map<> Map_t;
79  typedef Map_t::local_ordinal_type localId_t;
80  typedef Map_t::global_ordinal_type globalId_t;
81  typedef Tpetra::CrsMatrix<scalar_t, localId_t, globalId_t> Matrix_t;
82  typedef Tpetra::MultiVector<scalar_t, localId_t, globalId_t> MultiVector_t;
83  typedef Tpetra::Vector<scalar_t, localId_t, globalId_t> Vector_t;
84 
85  // Useful typedefs: Zoltan2 types
86  typedef Zoltan2::XpetraCrsMatrixAdapter<Matrix_t> MatrixAdapter_t;
87  typedef Zoltan2::XpetraMultiVectorAdapter<Vector_t> MultiVectorAdapter_t;
88 
89  // Input parameters with default values
90  std::string method = "scotch"; // Partitioning method
91  globalId_t nx = 50, ny = 40, nz = 30; // Dimensions of mesh corresponding to
92  // the matrix to be partitioned
93 
94  // Read run-time options.
95  Teuchos::CommandLineProcessor cmdp (false, false);
96  cmdp.setOption("method", &method,
97  "Partitioning method to use: scotch or parmetis.");
98  cmdp.setOption("nx", &nx,
99  "number of gridpoints in X dimension for "
100  "mesh used to generate matrix; must be >= 1.");
101  cmdp.setOption("ny", &ny,
102  "number of gridpoints in Y dimension for "
103  "mesh used to generate matrix; must be >= 1.");
104  cmdp.setOption("nz", &nz,
105  "number of gridpoints in Z dimension for "
106  "mesh used to generate matrix; must be >= 1.");
107  cmdp.parse(narg, arg);
108 
109  if ((nx < 1) || (ny < 1) || (nz < 1)) {
110  cout << "Input error: nx, ny and nz must be >= 1" << endl;
111  return -1;
112  }
113 
114  // For this example, generate a matrix using Galeri with the
115  // default Tpetra distribution:
116  // Laplace3D matrix corresponding to mesh with dimensions (nx X ny X nz)
117  // with block row-based distribution
118  Teuchos::ParameterList galeriList;
119  galeriList.set("nx", nx);
120  galeriList.set("ny", ny);
121  galeriList.set("nz", nz);
122  Tpetra::global_size_t nGlobalElements = nx * ny * nz;
123 
124  RCP<Matrix_t> origMatrix;
125 
126  try {
127  RCP<const Map_t> map = rcp(new Map_t(nGlobalElements, 0, comm));
128 
129  typedef Galeri::Xpetra::Problem<Map_t,Matrix_t,MultiVector_t> Galeri_t;
130  RCP<Galeri_t> galeriProblem =
131  Galeri::Xpetra::BuildProblem<scalar_t, localId_t, globalId_t,
132  Map_t, Matrix_t, MultiVector_t>
133  ("Laplace3D", map, galeriList);
134  origMatrix = galeriProblem->BuildMatrix();
135  }
136  catch (std::exception &e) {
137  cout << "Exception in Galeri matrix generation. " << e.what() << endl;
138  return -1;
139  }
140 
141  if (me == 0)
142  cout << "NumRows = " << origMatrix->getGlobalNumRows() << endl
143  << "NumNonzeros = " << origMatrix->getGlobalNumEntries() << endl
144  << "NumProcs = " << comm->getSize() << endl;
145 
146  // Create vectors to use with the matrix for sparse matvec.
147  RCP<Vector_t> origVector, origProd;
148  origProd = Tpetra::createVector<scalar_t,localId_t,globalId_t>(
149  origMatrix->getRangeMap());
150  origVector = Tpetra::createVector<scalar_t,localId_t,globalId_t>(
151  origMatrix->getDomainMap());
152  origVector->randomize();
153 
154  // Specify partitioning parameters
155  Teuchos::ParameterList param;
156  param.set("partitioning_approach", "partition");
157  param.set("algorithm", method);
158 
159  // Create an input adapter for the Tpetra matrix.
160  MatrixAdapter_t adapter(origMatrix);
161 
162  // Create and solve partitioning problem
163  Zoltan2::PartitioningProblem<MatrixAdapter_t> problem(&adapter, &param);
164 
165  try {
166  problem.solve();
167  }
168  catch (std::exception &e) {
169  cout << "Exception returned from solve(). " << e.what() << endl;
170  return -1;
171  }
172 
173  // Redistribute matrix and vector into new matrix and vector.
174  // Can use PartitioningSolution from matrix to redistribute the vectors, too.
175 
176  if (me == 0) cout << "Redistributing matrix..." << endl;
177  RCP<Matrix_t> redistribMatrix;
178  adapter.applyPartitioningSolution(*origMatrix, redistribMatrix,
179  problem.getSolution());
180 
181  if (me == 0) cout << "Redistributing vectors..." << endl;
182  RCP<Vector_t> redistribVector;
183  MultiVectorAdapter_t adapterVector(origVector);
184  adapterVector.applyPartitioningSolution(*origVector, redistribVector,
185  problem.getSolution());
186 
187  // Create a new product vector for sparse matvec
188  RCP<Vector_t> redistribProd;
189  redistribProd = Tpetra::createVector<scalar_t,localId_t,globalId_t>(
190  redistribMatrix->getRangeMap());
191 
192  // SANITY CHECK
193  // A little output for small problems
194  if (origMatrix->getGlobalNumRows() <= 50) {
195  cout << me << " ORIGINAL: ";
196  for (size_t i = 0; i < origVector->getLocalLength(); i++)
197  cout << origVector->getMap()->getGlobalElement(i) << " ";
198  cout << endl;
199  cout << me << " REDISTRIB: ";
200  for (size_t i = 0; i < redistribVector->getLocalLength(); i++)
201  cout << redistribVector->getMap()->getGlobalElement(i) << " ";
202  cout << endl;
203  }
204 
205  // SANITY CHECK
206  // check that redistribution is "correct"; perform matvec with
207  // original and redistributed matrices/vectors and compare norms.
208 
209  if (me == 0) cout << "Matvec original..." << endl;
210  origMatrix->apply(*origVector, *origProd);
211  scalar_t origNorm = origProd->norm2();
212  if (me == 0)
213  cout << "Norm of Original matvec prod: " << origNorm << endl;
214 
215  if (me == 0) cout << "Matvec redistributed..." << endl;
216  redistribMatrix->apply(*redistribVector, *redistribProd);
217  scalar_t redistribNorm = redistribProd->norm2();
218  if (me == 0)
219  cout << "Norm of Redistributed matvec prod: " << redistribNorm << endl;
220 
221  if (me == 0) {
222  const double epsilon = 0.00000001;
223  if (redistribNorm > origNorm+epsilon || redistribNorm < origNorm-epsilon)
224  cout << "Mat-Vec product changed; FAIL" << std::endl;
225  else
226  cout << "PASS" << endl;
227  }
228 
229  return 0;
230 }
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
size_t global_size_t
Defines the XpetraMultiVectorAdapter.
Defines the XpetraCrsMatrixAdapter class.
int main(int narg, char **arg)
Definition: graph.cpp:67
An adapter for Xpetra::MultiVector.
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
PartitioningProblem sets up partitioning problems for the user.
Defines the PartitioningProblem class.
#define epsilon
void solve(bool updateInputData=true)
Direct the problem to create a solution.