Zoltan2
Loading...
Searching...
No Matches
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_CrsMatrix.hpp>
54#include <Tpetra_Vector.hpp>
55#include <Galeri_XpetraMaps.hpp>
56#include <Galeri_XpetraProblemFactory.hpp>
57
58using Teuchos::RCP;
59
61// Program to demonstrate use of Zoltan2 to partition a TPetra matrix
62// using graph partitioning via Scotch or ParMETIS.
64
65int 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;
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}
Defines the PartitioningProblem class.
Defines the XpetraCrsMatrixAdapter class.
Defines the XpetraMultiVectorAdapter.
int main()
PartitioningProblem sets up partitioning problems for the user.
const PartitioningSolution< Adapter > & getSolution()
Get the solution to the problem.
void solve(bool updateInputData=true)
Direct the problem to create a solution.
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
void applyPartitioningSolution(const User &in, User *&out, const PartitioningSolution< Adapter > &solution) const
#define epsilon
Definition: nd.cpp:82
Zoltan2::XpetraMultiVectorAdapter< tMVector_t > MultiVectorAdapter_t
Definition: nd.cpp:80