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