Zoltan2
Loading...
Searching...
No Matches
Test_Sphynx.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 Seher Acer (sacer@sandia.gov)
39// Erik Boman (egboman@sandia.gov)
40// Siva Rajamanickam (srajama@sandia.gov)
41// Karen Devine (kddevin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
51#include <iostream>
52#include <limits>
53#include <Teuchos_ParameterList.hpp>
54#include <Teuchos_RCP.hpp>
55#include <Teuchos_FancyOStream.hpp>
56#include <Teuchos_CommandLineProcessor.hpp>
57#include <Tpetra_CrsMatrix.hpp>
58#include <Tpetra_Vector.hpp>
59#include <MatrixMarket_Tpetra.hpp>
60
61using Teuchos::RCP;
62
64// This program is a modified version of partitioning1.cpp (Karen Devine, 2011)
65// which can be found in zoltan2/test/core/partition/.
66// This version demonstrates use of Sphynx to partition a Tpetra matrix
67// (read from a MatrixMarket file or generated by Galeri::Xpetra).
68// Usage:
69// Zoltan2_Sphynx.exe [--inputFile=filename] [--outputFile=outfile] [--verbose]
70// [--x=#] [--y=#] [--z=#] [--matrix={Laplace1D,Laplace2D,Laplace3D}
71// [--normalized] [--generalized] [--polynomial]
73
75// Eventually want to use Teuchos unit tests to vary z2TestLO and
76// GO. For now, we set them at compile time based on whether Tpetra
77// is built with explicit instantiation on. (in Zoltan2_TestHelpers.hpp)
78
82
83typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO> SparseMatrix;
84typedef Tpetra::CrsGraph<z2TestLO, z2TestGO> SparseGraph;
85typedef Tpetra::Vector<z2TestScalar, z2TestLO, z2TestGO> VectorT;
86typedef VectorT::node_type Node;
87
91
92
93// Integer vector
94typedef Tpetra::Vector<int, z2TestLO, z2TestGO> IntVector;
96
97#define epsilon 0.00000001
98#define NNZ_IDX 1
99
101int main(int narg, char** arg)
102{
103 std::string inputFile = ""; // Matrix Market or Zoltan file to read
104 std::string outputFile = ""; // Matrix Market or Zoltan file to write
105 std::string inputPath = testDataFilePath; // Directory with input file
106 bool verbose = false; // Verbosity of output
107 bool distributeInput = true;
108 bool haveFailure = false;
109 int nVwgts = 0;
110 int testReturn = 0;
111
112 // Sphynx-related parameters
113 bool isNormalized = false;
114 bool isGeneralized = false;
115 std::string precType = "jacobi";
116 std::string initialGuess = "random";
117 bool useFullOrtho = true;
118
120 Tpetra::ScopeGuard tscope(&narg, &arg);
121 RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
122 int me = comm->getRank();
123 int commsize = comm->getSize();
124
125 // Read run-time options.
126 Teuchos::CommandLineProcessor cmdp (false, false);
127 cmdp.setOption("inputPath", &inputPath,
128 "Path to the MatrixMarket or Zoltan file to be read; "
129 "if not specified, a default path will be used.");
130 cmdp.setOption("inputFile", &inputFile,
131 "Name of the Matrix Market or Zoltan file to read; "
132 "if not specified, a matrix will be generated by MueLu.");
133 cmdp.setOption("outputFile", &outputFile,
134 "Name of the Matrix Market sparse matrix file to write, "
135 "echoing the input/generated matrix.");
136 cmdp.setOption("vertexWeights", &nVwgts,
137 "Number of weights to generate for each vertex");
138 cmdp.setOption("verbose", "quiet", &verbose,
139 "Print messages and results.");
140 cmdp.setOption("distribute", "no-distribute", &distributeInput,
141 "indicate whether or not to distribute "
142 "input across the communicator");
143 // Sphynx-related parameters:
144 cmdp.setOption("normalized", "combinatorial", &isNormalized,
145 "indicate whether or not to use a normalized Laplacian.");
146 cmdp.setOption("generalized", "non-generalized", &isGeneralized,
147 "indicate whether or not to use a generalized Laplacian.");
148 cmdp.setOption("precond", &precType,
149 "indicate which preconditioner to use [muelu|jacobi|polynomial].");
150 cmdp.setOption("initialGuess", &initialGuess,
151 "initial guess for LOBPCG");
152 cmdp.setOption("useFullOrtho", "partialOrtho", &useFullOrtho,
153 "use full orthogonalization.");
154
156 // Even with cmdp option "true", I get errors for having these
157 // arguments on the command line. (On redsky build)
158 // KDDKDD Should just be warnings, right? Code should still work with these
159 // KDDKDD params in the create-a-matrix file. Better to have them where
160 // KDDKDD they are used.
161 int xdim=10;
162 int ydim=10;
163 int zdim=10;
164 std::string matrixType("Laplace3D");
165
166 cmdp.setOption("x", &xdim,
167 "number of gridpoints in X dimension for "
168 "mesh used to generate matrix.");
169 cmdp.setOption("y", &ydim,
170 "number of gridpoints in Y dimension for "
171 "mesh used to generate matrix.");
172 cmdp.setOption("z", &zdim,
173 "number of gridpoints in Z dimension for "
174 "mesh used to generate matrix.");
175 cmdp.setOption("matrix", &matrixType,
176 "Matrix type: Laplace1D, Laplace2D, or Laplace3D");
178
179 cmdp.parse(narg, arg);
180
181 RCP<UserInputForTests> uinput;
182
183 if (inputFile != "") // Input file specified; read a matrix
184 uinput = rcp(new UserInputForTests(inputPath, inputFile, comm,
185 true, distributeInput));
186
187 else // Let MueLu generate a default matrix
188 uinput = rcp(new UserInputForTests(xdim, ydim, zdim, string(""), comm,
189 true, distributeInput));
190
191 RCP<SparseMatrix> origMatrix = uinput->getUITpetraCrsMatrix();
192
193 if (origMatrix->getGlobalNumRows() < 40) {
194 Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,false));
195 origMatrix->describe(out, Teuchos::VERB_EXTREME);
196 }
197
198
199 if (outputFile != "") {
200 // Just a sanity check.
201 Tpetra::MatrixMarket::Writer<SparseMatrix>::writeSparseFile(outputFile,
202 origMatrix, verbose);
203 }
204
205 if (me == 0)
206 std::cout << "NumRows = " << origMatrix->getGlobalNumRows() << std::endl
207 << "NumNonzeros = " << origMatrix->getGlobalNumEntries() << std::endl
208 << "NumProcs = " << comm->getSize() << std::endl
209 << "NumLocalRows (rank 0) = " << origMatrix->getLocalNumRows() << std::endl;
210
212 RCP<VectorT> origVector, origProd;
213 origProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
214 origMatrix->getRangeMap());
215 origVector = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
216 origMatrix->getDomainMap());
217 origVector->randomize();
218
220 Teuchos::RCP<Teuchos::ParameterList> params(new Teuchos::ParameterList);
221 params->set("num_global_parts", commsize);
222 Teuchos::RCP<Teuchos::ParameterList> sphynxParams(new Teuchos::ParameterList);
223 sphynxParams->set("sphynx_skip_preprocessing", true); // Preprocessing has not been implemented yet.
224 sphynxParams->set("sphynx_preconditioner_type", precType);
225 sphynxParams->set("sphynx_verbosity", verbose ? 1 : 0);
226 sphynxParams->set("sphynx_initial_guess", initialGuess);
227 sphynxParams->set("sphynx_use_full_ortho", useFullOrtho);
228 std::string problemType = "combinatorial";
229 if(isNormalized)
230 problemType = "normalized";
231 else if(isGeneralized)
232 problemType = "generalized";
233 sphynxParams->set("sphynx_problem_type", problemType); // Type of the eigenvalue problem.
234
236 Teuchos::RCP<SparseGraphAdapter> adapter = Teuchos::rcp( new SparseGraphAdapter(origMatrix->getCrsGraph(), nVwgts));
237
241
242 zscalar_t *vwgts = NULL;
243 if (nVwgts) {
244 // Test vertex weights with stride nVwgts.
245 size_t nrows = origMatrix->getLocalNumRows();
246 if (nrows) {
247 vwgts = new zscalar_t[nVwgts * nrows];
248 for (size_t i = 0; i < nrows; i++) {
249 size_t idx = i * nVwgts;
250 vwgts[idx] = zscalar_t(origMatrix->getRowMap()->getGlobalElement(i));
251 for (int j = 1; j < nVwgts; j++) vwgts[idx+j] = 1.;
252 }
253 for (int j = 0; j < nVwgts; j++) {
254 if (j != NNZ_IDX) adapter->setVertexWeights(&vwgts[j], nVwgts, j);
255 else adapter->setVertexWeightIsDegree(NNZ_IDX);
256 }
257 }
258 }
259
261 Zoltan2::SphynxProblem<SparseGraphAdapter> problem(adapter.get(), params.get(), sphynxParams);
262
263 try {
264 if (me == 0) std::cout << "Calling solve() " << std::endl;
265
266 problem.solve();
267
268 if (me == 0) std::cout << "Done solve() " << std::endl;
269 }
270 catch (std::runtime_error &e) {
271 delete [] vwgts;
272 std::cout << "Runtime exception returned from solve(): " << e.what();
273 if (!strncmp(e.what(), "BUILD ERROR", 11)) {
274 // Catching build errors as exceptions is OK in the tests
275 std::cout << " PASS" << std::endl;
276 return 0;
277 }
278 else {
279 // All other runtime_errors are failures
280 std::cout << " FAIL" << std::endl;
281 return -1;
282 }
283 }
284 catch (std::logic_error &e) {
285 delete [] vwgts;
286 std::cout << "Logic exception returned from solve(): " << e.what()
287 << " FAIL" << std::endl;
288 return -1;
289 }
290 catch (std::bad_alloc &e) {
291 delete [] vwgts;
292 std::cout << "Bad_alloc exception returned from solve(): " << e.what()
293 << " FAIL" << std::endl;
294 return -1;
295 }
296 catch (std::exception &e) {
297 delete [] vwgts;
298 std::cout << "Unknown exception returned from solve(). " << e.what()
299 << " FAIL" << std::endl;
300 return -1;
301 }
302
305 size_t checkNparts = comm->getSize();
306
307 size_t checkLength = origMatrix->getLocalNumRows();
308 const SparseGraphAdapter::part_t *checkParts = problem.getSolution().getPartListView();
309
310 // Check for load balance
311 size_t *countPerPart = new size_t[checkNparts];
312 size_t *globalCountPerPart = new size_t[checkNparts];
313 zscalar_t *wtPerPart = (nVwgts ? new zscalar_t[checkNparts*nVwgts] : NULL);
314 zscalar_t *globalWtPerPart = (nVwgts ? new zscalar_t[checkNparts*nVwgts] : NULL);
315 for (size_t i = 0; i < checkNparts; i++) countPerPart[i] = 0;
316 for (size_t i = 0; i < checkNparts * nVwgts; i++) wtPerPart[i] = 0.;
317
318 for (size_t i = 0; i < checkLength; i++) {
319 if (size_t(checkParts[i]) >= checkNparts)
320 std::cout << "Invalid Part " << checkParts[i] << ": FAIL" << std::endl;
321 countPerPart[checkParts[i]]++;
322 for (int j = 0; j < nVwgts; j++) {
323 if (j != NNZ_IDX)
324 wtPerPart[checkParts[i]*nVwgts+j] += vwgts[i*nVwgts+j];
325 else
326 wtPerPart[checkParts[i]*nVwgts+j] += origMatrix->getNumEntriesInLocalRow(i);
327 }
328 }
329 Teuchos::reduceAll<int, size_t>(*comm, Teuchos::REDUCE_SUM, checkNparts,
330 countPerPart, globalCountPerPart);
331 Teuchos::reduceAll<int, zscalar_t>(*comm, Teuchos::REDUCE_SUM,
332 checkNparts*nVwgts,
333 wtPerPart, globalWtPerPart);
334
335 size_t min = std::numeric_limits<std::size_t>::max();
336 size_t max = 0;
337 size_t sum = 0;
338 size_t minrank = 0, maxrank = 0;
339 for (size_t i = 0; i < checkNparts; i++) {
340 if (globalCountPerPart[i] < min) {min = globalCountPerPart[i]; minrank = i;}
341 if (globalCountPerPart[i] > max) {max = globalCountPerPart[i]; maxrank = i;}
342 sum += globalCountPerPart[i];
343 }
344
345 if (me == 0) {
346 float avg = (float) sum / (float) checkNparts;
347 std::cout << "Minimum count: " << min << " on rank " << minrank << std::endl;
348 std::cout << "Maximum count: " << max << " on rank " << maxrank << std::endl;
349 std::cout << "Average count: " << avg << std::endl;
350 std::cout << "Total count: " << sum
351 << (sum != origMatrix->getGlobalNumRows()
352 ? "Work was lost; FAIL"
353 : " ")
354 << std::endl;
355 std::cout << "Imbalance: " << max / avg << std::endl;
356 if (nVwgts) {
357 std::vector<zscalar_t> minwt(nVwgts, std::numeric_limits<zscalar_t>::max());
358 std::vector<zscalar_t> maxwt(nVwgts, 0.);
359 std::vector<zscalar_t> sumwt(nVwgts, 0.);
360 for (size_t i = 0; i < checkNparts; i++) {
361 for (int j = 0; j < nVwgts; j++) {
362 size_t idx = i*nVwgts+j;
363 if (globalWtPerPart[idx] < minwt[j]) minwt[j] = globalWtPerPart[idx];
364 if (globalWtPerPart[idx] > maxwt[j]) maxwt[j] = globalWtPerPart[idx];
365 sumwt[j] += globalWtPerPart[idx];
366 }
367 }
368 for (int j = 0; j < nVwgts; j++) {
369 float avgwt = (float) sumwt[j] / (float) checkNparts;
370 std::cout << std::endl;
371 std::cout << "Minimum weight[" << j << "]: " << minwt[j] << std::endl;
372 std::cout << "Maximum weight[" << j << "]: " << maxwt[j] << std::endl;
373 std::cout << "Average weight[" << j << "]: " << avgwt << std::endl;
374 std::cout << "Imbalance: " << maxwt[j] / avgwt << std::endl;
375 }
376 }
377 }
378
379 delete [] countPerPart;
380 delete [] wtPerPart;
381 delete [] globalCountPerPart;
382 delete [] globalWtPerPart;
383 delete [] vwgts;
384
386 if (me == 0) std::cout << "Redistributing matrix..." << std::endl;
387 SparseMatrix *redistribMatrix;
388 SparseMatrixAdapter adapterMatrix(origMatrix);
389 adapterMatrix.applyPartitioningSolution(*origMatrix, redistribMatrix,
390 problem.getSolution());
391 if (redistribMatrix->getGlobalNumRows() < 40) {
392 Teuchos::FancyOStream out(Teuchos::rcp(&std::cout,false));
393 redistribMatrix->describe(out, Teuchos::VERB_EXTREME);
394 }
395
396 if (me == 0) std::cout << "Redistributing vectors..." << std::endl;
397 VectorT *redistribVector;
398 MultiVectorAdapter adapterVector(origVector); //, weights, weightStrides);
399 adapterVector.applyPartitioningSolution(*origVector, redistribVector,
400 problem.getSolution());
401
402 RCP<VectorT> redistribProd;
403 redistribProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
404 redistribMatrix->getRangeMap());
405
406 // Test redistributing an integer vector with the same solution.
407 // This test is mostly to make sure compilation always works.
408 RCP<IntVector> origIntVec;
409 IntVector *redistIntVec;
410 origIntVec = Tpetra::createVector<int,z2TestLO,z2TestGO>(
411 origMatrix->getRangeMap());
412 for (size_t i = 0; i < origIntVec->getLocalLength(); i++)
413 origIntVec->replaceLocalValue(i, me);
414
415 IntVectorAdapter int_vec_adapter(origIntVec);
416 int_vec_adapter.applyPartitioningSolution(*origIntVec, redistIntVec,
417 problem.getSolution());
418 int origIntNorm = origIntVec->norm1();
419 int redistIntNorm = redistIntVec->norm1();
420 if (me == 0) std::cout << "IntegerVectorTest: " << origIntNorm << " == "
421 << redistIntNorm << " ?";
422 if (origIntNorm != redistIntNorm) {
423 if (me == 0) std::cout << " FAIL" << std::endl;
424 haveFailure = true;
425 }
426 else if (me == 0) std::cout << " OK" << std::endl;
427 delete redistIntVec;
428
431
432 if (me == 0) std::cout << "Matvec original..." << std::endl;
433 origMatrix->apply(*origVector, *origProd);
434 z2TestScalar origNorm = origProd->norm2();
435 if (me == 0)
436 std::cout << "Norm of Original matvec prod: " << origNorm << std::endl;
437
438 if (me == 0) std::cout << "Matvec redistributed..." << std::endl;
439 redistribMatrix->apply(*redistribVector, *redistribProd);
440 z2TestScalar redistribNorm = redistribProd->norm2();
441 if (me == 0)
442 std::cout << "Norm of Redistributed matvec prod: " << redistribNorm << std::endl;
443
444 if (redistribNorm > origNorm+epsilon || redistribNorm < origNorm-epsilon) {
445 testReturn = 1;
446 haveFailure = true;
447 }
448
449 delete redistribVector;
450 delete redistribMatrix;
451
452 if (me == 0) {
453 if (testReturn) {
454 std::cout << "Mat-Vec product changed; FAIL" << std::endl;
455 haveFailure = true;
456 }
457 if (!haveFailure)
458 std::cout << "PASS" << std::endl;
459 }
460
461 return testReturn;
462}
VectorT::node_type Node
#define epsilon
zlno_t z2TestLO
Tpetra::CrsGraph< z2TestLO, z2TestGO > SparseGraph
Zoltan2::XpetraMultiVectorAdapter< IntVector > IntVectorAdapter
#define NNZ_IDX
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
Tpetra::Vector< z2TestScalar, z2TestLO, z2TestGO > VectorT
Zoltan2::XpetraCrsMatrixAdapter< SparseMatrix > SparseMatrixAdapter
Zoltan2::XpetraMultiVectorAdapter< VectorT > MultiVectorAdapter
zgno_t z2TestGO
Tpetra::Vector< int, z2TestLO, z2TestGO > IntVector
zscalar_t z2TestScalar
Zoltan2::XpetraCrsGraphAdapter< SparseGraph > SparseGraphAdapter
common code used by tests
float zscalar_t
Tpetra::Map ::local_ordinal_type zlno_t
std::string testDataFilePath(".")
Tpetra::Map ::global_ordinal_type zgno_t
Defines XpetraCrsGraphAdapter class.
Defines the XpetraCrsMatrixAdapter class.
Defines the XpetraMultiVectorAdapter.
int main()
InputTraits< User >::part_t part_t
void solve(bool updateInputData=true)
Direct the problem to create a solution.
const PartitioningSolution< Adapter > & getSolution()
Provides access for Zoltan2 to Xpetra::CrsGraph data.
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
Definition coloring1.cpp:79
zscalar_t z2TestScalar
Definition coloring1.cpp:77
Tpetra::Vector< int, z2TestLO, z2TestGO > IntVector