Zoltan2
Loading...
Searching...
No Matches
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
63namespace Zoltan2{
64
65
79template <typename Adapter>
80class AlgMatrix : public Algorithm<Adapter>
81{
82
83public:
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
92 typedef typename Adapter::base_adapter_t base_adapter_t;
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
114private:
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
135template <typename Adapter>
136void AlgMatrix<Adapter>::partitionMatrix(const RCP<MatrixPartitioningSolution<Adapter> > &solution)
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
172 MatrixAdapter<user_t, userCoord_t> *adapterPtr =
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
Defines the PartitioningProblem class.
Defines the PartitioningSolution class.
AlgMatrix(const RCP< const Environment > &_env, const RCP< const Comm< int > > &_problemComm, const RCP< const MatrixAdapter< user_t, userCoord_t > > &_matrixAdapter)
Adapter::base_adapter_t base_adapter_t
void partitionMatrix(const RCP< MatrixPartitioningSolution< Adapter > > &solution)
Matrix Partitioning method.
Algorithm defines the base class for all algorithms.
virtual void getRowIDsView(const gno_t *&rowIds) const
Sets pointer to this process' rows' global IDs.
PartitioningProblem sets up partitioning problems for the user.
A PartitioningSolution is a solution to a partitioning problem.
Created by mbenlioglu on Aug 31, 2020.