Zoltan2
Loading...
Searching...
No Matches
orderingScotch.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
48#include <iostream>
49#include <fstream>
50#include <limits>
51#include <vector>
52#include <Teuchos_ParameterList.hpp>
53#include <Teuchos_RCP.hpp>
54#include <Teuchos_CommandLineProcessor.hpp>
55#include <Tpetra_CrsMatrix.hpp>
56#include <Tpetra_Vector.hpp>
57#include <MatrixMarket_Tpetra.hpp>
58
59using Teuchos::RCP;
60
62// Program to demonstrate use of Zoltan2 to order a TPetra matrix
63// (read from a MatrixMarket file or generated by Galeri::Xpetra).
64// Usage:
65// a.out [--inputFile=filename] [--outputFile=outfile] [--verbose]
66// [--x=#] [--y=#] [--z=#] [--matrix={Laplace1D,Laplace2D,Laplace3D}
67// Karen Devine, 2011
69
71// Eventually want to use Teuchos unit tests to vary z2TestLO and
72// GO. For now, we set them at compile time based on whether Tpetra
73// is built with explicit instantiation on. (in Zoltan2_TestHelpers.hpp)
74
78
79typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO> SparseMatrix;
80typedef Tpetra::Vector<z2TestScalar, z2TestLO, z2TestGO> Vector;
81typedef Vector::node_type Node;
83
84
85// Using perm
86size_t computeBandwidth(RCP<SparseMatrix> A, z2TestLO *perm)
87{
88 z2TestLO ii, i, j, k;
89 typename SparseMatrix::local_inds_host_view_type indices;
90 typename SparseMatrix::values_host_view_type values;
91
92 z2TestLO bw_left = 0;
93 z2TestLO bw_right = 0;
94
95 z2TestLO n = A->getLocalNumRows();
96
97 // Loop over rows of matrix
98 for (ii=0; ii<n; ii++) {
99 A->getLocalRowView (ii, indices, values);
100 for (k=0; k< static_cast<z2TestLO>(indices.size()); k++){
101 if (indices[k] < n){ // locally owned
102 if (perm){
103 i = perm[ii];
104 j = perm[indices[k]];
105 } else {
106 i = ii;
107 j = indices[k];
108 }
109 if (j-i > bw_right)
110 bw_right = j-i;
111 if (i-j > bw_left)
112 bw_left = i-j;
113 }
114 }
115 }
116 // Total bandwidth is the sum of left and right + 1
117 return (bw_left + bw_right + 1);
118}
119
120#define MDM
121#ifdef MDM
122// This is all temp code to be deleted when refactoring is compelte.
124 RCP<SparseMatrix> origMatrix, Zoltan2::LocalOrderingSolution<z2TestLO> *soln)
125{
126 typedef typename SparseMatrixAdapter::lno_t lno_t;
127
128 lno_t * perm = soln->getPermutationView();
129 lno_t * iperm = soln->getPermutationView(true);
130
131 lno_t numRows = origMatrix->getLocalNumRows();
132
133 std::cout << "origMatrix->getLocalNumRows(): " << numRows << std::endl;
134
135 if (numRows == 0) {
136 std::cout << "Skipping analysis - matrix is empty" << std::endl;
137 return;
138 }
139
140 Array<lno_t> oldMatrix(numRows*numRows);
141 Array<lno_t> newMatrix(numRows*numRows);
142
143 // print the solution permutation
144 lno_t showSize = numRows;
145 if(showSize > 30) {
146 showSize = 30;
147 }
148
149 std::cout << std::endl << "perm: ";
150 for(lno_t n = 0; n < showSize; ++n) {
151 std::cout << " " << perm[n] << " ";
152 }
153 if(showSize < numRows) {
154 std::cout << " ..."; // partial print...
155 }
156 std::cout << std::endl << "iperm: ";
157 for(lno_t n = 0; n < showSize; ++n) {
158 std::cout << " " << iperm[n] << " ";
159 }
160 if(showSize < numRows) {
161 std::cout << " ..."; // partial print...
162 }
163
164 std::cout << std::endl;
165
166 // write 1's in our matrix
167 for (lno_t j = 0; j < numRows; ++j) {
168 typename SparseMatrix::local_inds_host_view_type indices;
169 typename SparseMatrix::values_host_view_type wgts;
170 origMatrix->getLocalRowView( j, indices, wgts );
171 for (lno_t n = 0; n < static_cast<lno_t>(indices.size()); ++n) {
172 lno_t i = indices[n];
173 if (i < numRows) { // locally owned
174 oldMatrix[i + j*numRows] = 1;
175 newMatrix[perm[i] + perm[j]*numRows] = 1;
176 }
177 }
178 }
179
180 // print oldMatrix
181 std::cout << std::endl << "original graph in matrix form:" << std::endl;
182 for(lno_t y = 0; y < showSize; ++y) {
183 for(lno_t x = 0; x < showSize; ++x) {
184 std::cout << " " << oldMatrix[x + y*numRows];
185 }
186 if(showSize < numRows) {
187 std::cout << " ..."; // partial print...
188 }
189 std::cout << std::endl;
190 }
191 if(showSize < numRows) {
192 for(lno_t i = 0; i < showSize; ++i) {
193 std::cout << " ."; // partial print...
194 }
195 std::cout << " ..."; // partial print...
196 }
197 std::cout << std::endl;
198
199 // print newMatrix
200 std::cout << std::endl << "reordered graph in matrix form:" << std::endl;
201 for(lno_t y = 0; y < showSize; ++y) {
202 for(lno_t x = 0; x < showSize; ++x) {
203 std::cout << " " << newMatrix[x + y*numRows];
204 }
205 if(showSize < numRows) {
206 std::cout << " ..."; // partial print...
207 }
208 std::cout << std::endl;
209 }
210 if(showSize < numRows) {
211 for(lno_t i = 0; i < showSize; ++i) {
212 std::cout << " ."; // partial print...
213 }
214 std::cout << " ..."; // partial print...
215 }
216 std::cout << std::endl;
217}
218#endif
219
221int mainExecute(int narg, char** arg, RCP<const Teuchos::Comm<int> > comm)
222{
223 std::string inputFile = ""; // Matrix Market file to read
224 std::string outputFile = ""; // Output file to write
225 bool verbose = false; // Verbosity of output
226 int testReturn = 0;
227
228 int rank = comm->getRank();
229
230 // Read run-time options.
231 Teuchos::CommandLineProcessor cmdp (false, false);
232 cmdp.setOption("inputFile", &inputFile,
233 "Name of a Matrix Market file in the data directory; "
234 "if not specified, a matrix will be generated by Galeri.");
235 cmdp.setOption("outputFile", &outputFile,
236 "Name of file to write the permutation");
237 cmdp.setOption("verbose", "quiet", &verbose,
238 "Print messages and results.");
239
240 if (rank == 0 ) {
241 std::cout << "Starting everything" << std::endl;
242 }
243
245 // Even with cmdp option "true", I get errors for having these
246 // arguments on the command line. (On redsky build)
247 // KDDKDD Should just be warnings, right? Code should still work with these
248 // KDDKDD params in the create-a-matrix file. Better to have them where
249 // KDDKDD they are used.
250 int xdim=10;
251 int ydim=10;
252 int zdim=10;
253 std::string matrixType("Laplace3D");
254
255 cmdp.setOption("x", &xdim,
256 "number of gridpoints in X dimension for "
257 "mesh used to generate matrix.");
258 cmdp.setOption("y", &ydim,
259 "number of gridpoints in Y dimension for "
260 "mesh used to generate matrix.");
261 cmdp.setOption("z", &zdim,
262 "number of gridpoints in Z dimension for "
263 "mesh used to generate matrix.");
264 cmdp.setOption("matrix", &matrixType,
265 "Matrix type: Laplace1D, Laplace2D, or Laplace3D");
266
268 // Ordering options to test.
270 std::string orderMethod("scotch"); // NYI
271 cmdp.setOption("order_method", &orderMethod,
272 "order_method: natural, random, rcm, scotch");
273
274 std::string orderMethodType("local");
275 cmdp.setOption("order_method_type", &orderMethodType,
276 "local or global or both" );
277
279 cmdp.parse(narg, arg);
280
281
282 RCP<UserInputForTests> uinput;
283
284 // MDM - temp testing code
285 // testDataFilePath = "/Users/micheldemessieres/Documents/trilinos-build/trilinosrepo/parZoltan2/packages/zoltan2/test/driver/driverinputs/ordering";
286 // inputFile = "orderingTest";
287
288 if (inputFile != ""){ // Input file specified; read a matrix
289 uinput = rcp(new UserInputForTests(testDataFilePath, inputFile, comm, true));
290 }
291 else { // Let Galeri generate a matrix
292 uinput = rcp(
293 new UserInputForTests(xdim, ydim, zdim, matrixType, comm, true, true));
294 }
295
296 RCP<SparseMatrix> origMatrix = uinput->getUITpetraCrsMatrix();
297
298 if (rank == 0) {
299 std::cout << "NumRows = " << origMatrix->getGlobalNumRows() << std::endl
300 << "NumNonzeros = " << origMatrix->getGlobalNumEntries() << std::endl
301 << "NumProcs = " << comm->getSize() << std::endl;
302 }
303
305 // Currently Not Used
306 /*
307 RCP<Vector> origVector, origProd;
308 origProd = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
309 origMatrix->getRangeMap());
310 origVector = Tpetra::createVector<z2TestScalar,z2TestLO,z2TestGO>(
311 origMatrix->getDomainMap());
312 origVector->randomize();
313 */
314
316 Teuchos::ParameterList params;
317 params.set("order_method", orderMethod);
318 params.set("order_method_type", orderMethodType);
319
321 SparseMatrixAdapter adapter(origMatrix);
322
324
325 try {
326 Zoltan2::OrderingProblem<SparseMatrixAdapter> problem(&adapter, &params,
327 comm);
328
329 if( rank == 0 ) {
330 std::cout << "Going to solve" << std::endl;
331 }
332 problem.solve();
333
335 size_t checkLength;
336 z2TestLO *checkPerm, *checkInvPerm;
338 problem.getLocalOrderingSolution();
339
340 if( rank == 0 ) {
341 std::cout << "Going to get results" << std::endl;
342 }
343
344 #ifdef MDM // Temp debugging code all of which can be removed for final
345 for( int checkRank = 0; checkRank < comm->getSize(); ++checkRank ) {
346 comm->barrier();
347 if( checkRank == comm->getRank() ) {
348 std::cout << "Debug output matrix for rank: " << checkRank << std::endl;
349 tempDebugTest(origMatrix, soln);
350 }
351 comm->barrier();
352 }
353 #endif
354
355 // Permutation
356 checkLength = soln->getPermutationSize();
357 checkPerm = soln->getPermutationView();
358 checkInvPerm = soln->getPermutationView(true); // get permutation inverse
359
360 // Separators.
361 // The following methods needs to be supported:
362 // haveSeparators: true if Scotch Nested Dissection was called.
363 // getCBlkPtr: *CBlkPtr from Scotch_graphOrder
364 // getRangeTab: RangeTab from Scotch_graphOrder
365 // getTreeTab: TreeTab from Scotch_graphOrder
366 z2TestLO NumBlocks = 0;
367 z2TestLO *RangeTab;
368 z2TestLO *TreeTab;
369 if (soln->haveSeparators()) {
370 NumBlocks = soln->getNumSeparatorBlocks(); // BDD
371 RangeTab = soln->getSeparatorRangeView(); // BDD
372 TreeTab = soln->getSeparatorTreeView(); // BDD
373 }
374 else {
375 // TODO FAIL with error
376 }
377
378 if (outputFile != "") {
379 std::ofstream permFile;
380
381 // Write permutation (0-based) to file
382 // each process writes local perm to a separate file
383 //std::string fname = outputFile + "." + me;
384 std::stringstream fname;
385 fname << outputFile << "." << comm->getSize() << "." << rank;
386 permFile.open(fname.str().c_str());
387 for (size_t i=0; i<checkLength; i++){
388 permFile << " " << checkPerm[i] << std::endl;
389 }
390 permFile.close();
391 }
392
393 // Validate that checkPerm is a permutation
394 if (rank == 0 ) {
395 std::cout << "Checking permutation" << std::endl;
396 }
397
398 testReturn = soln->validatePerm();
399 if (testReturn) return testReturn;
400
401 // Validate the inverse permutation.
402 if (rank == 0 ) {
403 std::cout << "Checking inverse permutation" << std::endl;
404 }
405 for (size_t i=0; i< checkLength; i++){
406 testReturn = (checkInvPerm[checkPerm[i]] != z2TestLO(i));
407 if (testReturn) return testReturn;
408 }
409
410 // Validate NumBlocks
411 if (rank == 0 ) {
412 std::cout << "Checking num blocks" << std::endl;
413 }
414 testReturn = !((NumBlocks > 0) && (NumBlocks<z2TestLO(checkLength)));
415 if (testReturn) return testReturn;
416
417 // Validate RangeTab.
418 // Should be monitonically increasing, RT[0] = 0; RT[NumBlocks+1]=nVtx;
419 if (rank == 0 ) {
420 std::cout << "Checking range" << std::endl;
421 }
422 testReturn = RangeTab[0];
423 if (testReturn) return testReturn;
424
425 for (z2TestLO i = 0; i < NumBlocks; i++){
426 testReturn = !(RangeTab[i] < RangeTab[i+1]);
427 if (testReturn) return testReturn;
428 }
429
430 // How do we validate TreeTab?
431 // TreeTab root has -1, other values < NumBlocks
432 if (rank == 0 ) {
433 std::cout << "Checking Separator Tree" << std::endl;
434 }
435
436 if (checkLength) {
437 testReturn = (TreeTab[0] != -1);
438 }
439
440 if (testReturn) {
441 std::cout << "TreeTab[0] = " << TreeTab[0] << " != -1" << std::endl;
442 return testReturn;
443 }
444
445 for (size_t i=1; i< checkLength; i++){
446 testReturn = !(TreeTab[i] < NumBlocks);
447 if (testReturn) {
448 std::cout << "TreeTab[" << i << "] = " << TreeTab[i] << " >= NumBlocks "
449 << " = " << NumBlocks << std::endl;
450 return testReturn;
451 }
452 }
453
454 if (rank == 0 ) {
455 std::cout << "Going to compute the bandwidth" << std::endl;
456 }
457
458 // Compute original and permuted bandwidth
459 z2TestLO originalBandwidth = computeBandwidth(origMatrix, nullptr);
460 z2TestLO permutedBandwidth = computeBandwidth(origMatrix, checkPerm);
461
462 if (rank == 0 ) {
463 std::cout << "Original Bandwidth: " << originalBandwidth << std::endl;
464 std::cout << "Permuted Bandwidth: " << permutedBandwidth << std::endl;
465 }
466
467 if(permutedBandwidth >= originalBandwidth) {
468 if (rank == 0 ) {
469 std::cout << "Test failed: bandwidth was not improved!" << std::endl;
470
471 std::cout << "Test in progress - returning OK for now..." << std::endl;
472 }
473
474 // return 1; // TO DO - we need test to have proper ordering
475 }
476 else {
477 if (rank == 0) {
478 std::cout << "The bandwidth was improved. Good!" << std::endl;
479 }
480 }
481 }
482 catch (std::exception &e) {
483 std::cout << "Exception caught in ordering" << std::endl;
484 std::cout << e.what() << std::endl;
485 return 1;
486 }
487
488 return 0;
489}
490
491int main(int narg, char** arg)
492{
493 Tpetra::ScopeGuard tscope(&narg, &arg);
494 Teuchos::RCP<const Teuchos::Comm<int> > comm = Tpetra::getDefaultComm();
495
496 int result = mainExecute(narg, arg, comm);
497
498 // get reduced return code form all procsses
499 comm->barrier();
500 int resultReduced;
501 reduceAll<int,int>(*comm, Teuchos::EReductionType::REDUCE_MAX, 1,
502 &result, &resultReduced);
503
504 // provide a final message which guarantees that the test will fail
505 // if any of the processes failed
506 if (comm->getRank() == 0) {
507 std::cout << "Scotch Ordering test with " << comm->getSize()
508 << " processes has test return code " << resultReduced
509 << " and is exiting in the "
510 << ((resultReduced == 0 ) ? "PASSED" : "FAILED") << " state."
511 << std::endl;
512 }
513}
514
Defines the OrderingProblem class.
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 the XpetraCrsMatrixAdapter class.
int main()
InputTraits< User >::lno_t lno_t
OrderingProblem sets up ordering problems for the user.
void solve(bool updateInputData=true)
Direct the problem to create a solution.
LocalOrderingSolution< lno_t > * getLocalOrderingSolution()
Get the local ordering solution to the problem.
index_t * getSeparatorRangeView() const
Get pointer to separator range.
bool haveSeparators() const
Do we have the separators?
index_t getNumSeparatorBlocks() const
Get number of separator column blocks.
int validatePerm()
returns 0 if permutation is valid, negative if invalid.
index_t * getPermutationView(bool inverse=false) const
Get pointer to permutation. If inverse = true, return inverse permutation. By default,...
size_t getPermutationSize() const
Get (local) size of permutation.
index_t * getSeparatorTreeView() const
Get pointer to separator tree.
Provides access for Zoltan2 to Xpetra::CrsMatrix data.
zlno_t z2TestLO
Definition: coloring1.cpp:75
map_t::local_ordinal_type lno_t
Definition: mapRemotes.cpp:17
Vector::node_type Node
zlno_t z2TestLO
void tempDebugTest(RCP< SparseMatrix > origMatrix, Zoltan2::LocalOrderingSolution< z2TestLO > *soln)
int mainExecute(int narg, char **arg, RCP< const Teuchos::Comm< int > > comm)
Tpetra::CrsMatrix< z2TestScalar, z2TestLO, z2TestGO > SparseMatrix
size_t computeBandwidth(RCP< SparseMatrix > A, z2TestLO *perm)
Zoltan2::XpetraCrsMatrixAdapter< SparseMatrix > SparseMatrixAdapter
Tpetra::Vector< z2TestScalar, z2TestLO, z2TestGO > Vector
zgno_t z2TestGO
zscalar_t z2TestScalar