Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
Zoltan2_MatrixPartitioningAlgs.hpp
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 #ifndef _ZOLTAN2_ALGMATRIX_HPP_
46 #define _ZOLTAN2_ALGMATRIX_HPP_
47 
48 // #include <Zoltan2_IdentifierModel.hpp>
51 // #include <Zoltan2_Algorithm.hpp>
52 
54 
55 // #include <sstream>
56 // #include <string>
57 // #include <bitset>
58 
63 namespace Zoltan2{
64 
65 
79 template <typename Adapter>
80 class AlgMatrix : public Algorithm<Adapter>
81 {
82 
83 public:
84  typedef typename Adapter::lno_t lno_t; // local ids
85  typedef typename Adapter::gno_t gno_t; // global ids
86  typedef typename Adapter::scalar_t scalar_t; // scalars
87  typedef typename Adapter::part_t part_t; // part numbers
88 
89  typedef typename Adapter::user_t user_t;
90  typedef typename Adapter::userCoord_t userCoord_t;
91 
93 
94 
95 
96  // Constructor
97  AlgMatrix(const RCP<const Environment> &_env,
98  const RCP<const Comm<int> > &_problemComm,
99  const RCP<const MatrixAdapter<user_t, userCoord_t> > &_matrixAdapter)
100  : mEnv(_env), mProblemComm(_problemComm), mMatrixAdapter(_matrixAdapter)
101  {}
102 
103 
104  // AlgMatrix(const RCP<const Environment> &_env,
105  // const RCP<const Comm<int> > &_problemComm,
106  // const RCP<const Adapter> &_matrixAdapter)
107  // : mEnv(_env), mProblemComm(_problemComm), mMatrixAdapter(_matrixAdapter)
108  // {}
109 
110  // Partitioning method
111  void partitionMatrix(const RCP<MatrixPartitioningSolution<Adapter> > &solution);
112 
113 
114 private:
115  const RCP<const Environment> mEnv;
116  const RCP<const Comm<int> > mProblemComm;
117 
118  const RCP<const MatrixAdapter<user_t, userCoord_t> > mMatrixAdapter;
119  //const RCP<const Adapter > mMatrixAdapter;
120 
121  //typedef Tpetra::CrsMatrix<z2TestScalar, z2TestLO, z2TestGO> SparseMatrix;
122  // const RCP<const XpetraCrsMatrixAdapter<user_t, userCoord_t> > mMatrixAdapter;
123 
124 
125 
126 
127 };
128 
129 
131 // Partitioning method
132 //
133 // -- For now implementing 2D Random Cartesian
135 template <typename Adapter>
137 {
138  mEnv->debug(DETAILED_STATUS, std::string("Entering AlgBlock"));
139 
140  int myrank = mEnv->myRank_;
141  int nprocs = mEnv->numProcs_;
142 
144  // Determine processor dimension d for 2D block partitioning d = sqrt(p)
145  // Determine processor row and processor column for this rank for 2D partitioning
147  int procDim = sqrt(nprocs);
148  assert(procDim * procDim == nprocs); // TODO: Should check this earlier and produce more useful error
149 
150  int myProcRow = myrank / procDim;
151  int myProcCol = myrank % procDim;
153 
155  // Create 1D Random partitioning problem to create partitioning for the 2D partitioning
157 
159  // Create parameters for 1D partitioning
161  Teuchos::ParameterList params1D("Params for 1D partitioning");
162  // params1D.set("algorithm", "random"); //TODO add support for random
163  params1D.set("algorithm", "block");
164  params1D.set("imbalance_tolerance", 1.1);
165  params1D.set("num_global_parts", procDim);
167 
169  // Create Zoltan2 partitioning problem
170  // -- Alternatively we could simply call algorithm directly
173  const_cast<MatrixAdapter<user_t, userCoord_t>*>(mMatrixAdapter.getRawPtr());
174 
175 
177  new Zoltan2::PartitioningProblem<base_adapter_t>(adapterPtr, &params1D);
178 
180 
182  // Solve the problem, get the solution
184  problem->solve();
185 
186  // Zoltan2::PartitioningSolution<SparseMatrixAdapter_t> solution = problem.getSolution();
187  const Zoltan2::PartitioningSolution<base_adapter_t> &solution1D = problem->getSolution();
188 
190 
192  // Get part assignments for each local_id
194  // size_t numGlobalParts = solution1D->getTargetGlobalNumberOfParts();
195 
196  const part_t *parts = solution1D.getPartListView();
197 
199 
201 
202 
204  // Create column Ids ArrayRCP colIDs to store in solution
205  // This will define which column IDs are allowed in a particular processor
206  // column.
208 
210  // Group gids corresponding to local ids based on process column
212  const gno_t *rowIds;
213  mMatrixAdapter->getRowIDsView(rowIds);
214 
215  size_t nLocIDs = mMatrixAdapter->getLocalNumRows();
216 
217  std::vector<std::vector<gno_t> > idsInProcColI(procDim);
218  Teuchos::ArrayRCP<gno_t> colIDs;
219 
220  for(size_t i=0; i<nLocIDs; i++)
221  {
222  // Nonzeros with columns of index rowIds[i] belong to some processor
223  // in processor column parts[i]
224  idsInProcColI[ parts[i] ].push_back(rowIds[i]);
225  }
227 
228  delete problem; // delete 1D partitioning problem
229 
230 
232  // Communicate gids to process col roots
234 
235  // Loop over each of the processor columns
236  for(int procColI=0; procColI<procDim; procColI++)
237  {
238 
240  // This rank is root of procColI, need to receive
242  if(myProcCol==procColI && myProcRow==0)
243  {
244  assert(myrank==procColI); // Could make this above conditional?
245 
246 
247  std::set<gno_t> gidSet; // to get sorted list of gids
248 
250  // Copy local gids to sorted gid set
252  for(int i=0;i<idsInProcColI[procColI].size();i++)
253  {
254  gidSet.insert(idsInProcColI[procColI][i]);
255  }
257 
259  // Receive gids from remote processors, insert into set
261  std::vector<gno_t> recvBuf;
262  for(int src=0;src<nprocs;src++)
263  {
264  if(src!=myrank)
265  {
266  int buffSize;
267 
268  Teuchos::receive<int,int>(*mProblemComm, src, 1, &buffSize);
269 
270  if(buffSize>0)
271  {
272  recvBuf.resize(buffSize);
273 
274  Teuchos::receive<int,gno_t>(*mProblemComm, src, buffSize, recvBuf.data());
275 
276  for(int i=0;i<buffSize;i++)
277  {
278  gidSet.insert(recvBuf[i]);
279  }
280  }
281  }
282  }
284 
286  //Copy data to std::vector
288  colIDs.resize(gidSet.size());
289 
290  typename std::set<gno_t>::const_iterator iter;
291  int indx=0;
292  for (iter=gidSet.begin(); iter!=gidSet.end(); ++iter)
293  {
294  colIDs[indx] = *iter;
295  indx++;
296  }
298  }
300 
302  // This rank is not root of procColI, need to senddata block to root
304  else
305  {
306  int sizeToSend = idsInProcColI[procColI].size();
307 
308  //is the dst proc info correct here?
309 
310  // Root of procColI is rank procColI
311  Teuchos::send<int,int>(*mProblemComm,1,&sizeToSend,procColI);
312 
313  Teuchos::send<int,gno_t>(*mProblemComm,sizeToSend,idsInProcColI[procColI].data(),procColI);
314  }
316 
317  }
319 
320  // Free memory if possible
321 
322 
324  // Communicate gids from process col root down process column to other procs
326 
328  // This rank is root of its processor column, need to send
330  if(myProcRow==0)
331  {
333  // Send data to remote processors in processor column
335  for(int procColRank=0; procColRank<procDim; procColRank++)
336  {
337  // Convert from proc column rank to global rank
338  int dstRank = procColRank*procDim + myProcCol;
339 
340  if(dstRank!=myrank)
341  {
342  int sizeToSend = colIDs.size();
343 
344  Teuchos::send<int,int>(*mProblemComm,1,&sizeToSend,dstRank);
345  Teuchos::send<int,gno_t>(*mProblemComm,sizeToSend,colIDs.get(),dstRank);
346 
347  }
348  }
350 
351  }
353 
355  // This rank is not root of processor, need to recv data from root
357  else
358  {
359  // Src rank is root of processor column = myProcCol
360  int srcRank = myProcCol;
361 
362  int buffSize;
363  Teuchos::receive<int,int>(*mProblemComm, srcRank, 1, &buffSize);
364 
365  colIDs.resize(buffSize);
366  Teuchos::receive<int,gno_t>(*mProblemComm, srcRank, buffSize, colIDs.get());
367  }
369 
371  // Create domain/range IDs (same for both now) Array RCP to store in solution
372  // Created from equal division of column IDs
374 
375  // Split processor column ids into nDim parts
376  // Processor column i ids split between ranks 0+i, nDim+i, 2*nDim+i, ...
377 
378  //
379 
380  ArrayRCP<gno_t> domRangeIDs; // Domain/Range vector Ids assigned to this process
381 
382  size_t nCols = colIDs.size();
383  lno_t nColsDivDim = nCols / procDim;
384  lno_t nColsModDim = nCols % procDim;
385 
386  size_t sColIndx;
387  size_t domRangeSize;
388 
389  // This proc will have nColsDivDim+1 domain/range ids
390  if(myProcRow < nColsModDim)
391  {
392  sColIndx = myProcRow * (nColsDivDim+1);
393  domRangeSize = nColsDivDim+1;
394  }
395  // This proc will have nColsDivDim domain/range ids
396  else
397  {
398  sColIndx = nColsModDim*(nColsDivDim+1) + (myProcRow-nColsModDim) * nColsDivDim;
399  domRangeSize = nColsDivDim;
400  }
401 
402  domRangeIDs.resize(domRangeSize);
403 
404  for(size_t i=0;i<domRangeSize;i++)
405  {
406  domRangeIDs[i] = colIDs[sColIndx+i];
407  }
409 
411  // Create row IDs Array RCP to store in solution
412  // Created from union of domRangeIDs
413  //
414  // Procs 0, 1, ... nDim-1 share the same set of rowIDs
415  // Procs nDim, nDim+1, ..., 2*nDim-1 share the same set of rowIDs
417 
419  // Create subcommunicator for processor row
421  std::vector<int> subcommRanks(procDim);
422 
423  for(unsigned int i=0; i<procDim; i++)
424  {
425  subcommRanks[i] = myProcRow*procDim + i;
426  }
427 
428  ArrayView<const int> subcommRanksView = Teuchos::arrayViewFromVector (subcommRanks);
429 
430  RCP<Teuchos::Comm<int> > rowcomm = mProblemComm->createSubcommunicator(subcommRanksView);
432 
434  // Determine max number of columns in this process row
436  size_t maxProcRowNCols=0;
437 
438  Teuchos::reduceAll<int, size_t>(*rowcomm,Teuchos::REDUCE_MAX,1,&domRangeSize,&maxProcRowNCols);
440 
442  // Communicate all domRangeIDs to all processes in procRow
444  gno_t MAXVAL = std::numeric_limits<gno_t>::max();
445 
446  std::vector<gno_t> locRowIDs(maxProcRowNCols,MAXVAL);
447 
448  Teuchos::ArrayRCP<gno_t> rowIDs(maxProcRowNCols * procDim);
449 
450  // Copy local domRangeIDs into local row IDs, "extra elements" will have
451  // value MAXVAL
452  for(size_t i=0;i<domRangeIDs.size();i++)
453  {
454  locRowIDs[i] = domRangeIDs[i];
455  }
456 
457  // Gather all ids onto all processes in procRow
458  Teuchos::gatherAll<int,gno_t>(*rowcomm,maxProcRowNCols, locRowIDs.data(),
459  maxProcRowNCols*procDim, rowIDs.get());
461 
462  // Free local data
463  std::vector<int>().swap(locRowIDs);
464 
465 
467  // Insert local domRangeIDs into row IDs set, filter out extra items
469  std::set<gno_t> setRowIDs;
470 
471  for(size_t i=0;i<rowIDs.size();i++)
472  {
473  if(rowIDs[i]!=MAXVAL)
474  {
475  setRowIDs.insert(rowIDs[i]);
476  }
477  }
479 
481  // Resize rowIDs array and copy data to array (now sorted)
483  rowIDs.resize(setRowIDs.size());
484 
485  typename std::set<gno_t>::const_iterator iter;
486  size_t indx=0;
487 
488  for(iter=setRowIDs.begin(); iter!=setRowIDs.end(); ++iter)
489  {
490  rowIDs[indx] = *iter;
491  indx++;
492  }
494 
496 
498  // Finished. Store partition in solution.
500  solution->setIDLists(rowIDs,colIDs,domRangeIDs,domRangeIDs);
502 
503  mEnv->debug(DETAILED_STATUS, std::string("Exiting AlgMatrix"));
504 }
506 
507 
508 
509 } // namespace Zoltan2
510 
511 #endif
Zoltan2::BaseAdapter< userTypes_t > base_adapter_t
MatrixAdapter defines the adapter interface for matrices.
Defines the PartitioningSolution class.
sub-steps, each method&#39;s entry and exit
SparseMatrixAdapter_t::part_t part_t
AlgMatrix(const RCP< const Environment > &_env, const RCP< const Comm< int > > &_problemComm, const RCP< const MatrixAdapter< user_t, userCoord_t > > &_matrixAdapter)
A PartitioningSolution is a solution to a partitioning problem.
const part_t * getPartListView() const
Returns the part list corresponding to the global ID list.
Algorithm defines the base class for all algorithms.
void partitionMatrix(const RCP< MatrixPartitioningSolution< Adapter > > &solution)
Matrix Partitioning method.
PartitioningProblem sets up partitioning problems for the user.
Defines the PartitioningProblem class.
Adapter::base_adapter_t base_adapter_t
A PartitioningSolution is a solution to a partitioning problem.
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t
Definition: Metric.cpp:74