47 #include "Teuchos_CommandLineProcessor.hpp"
48 #include "Tpetra_CrsMatrix.hpp"
49 #include "Tpetra_Core.hpp"
50 #include "Tpetra_KokkosCompat_DefaultNode.hpp"
56 #include "Teuchos_TimeMonitor.hpp"
57 #include "Teuchos_StackedTimer.hpp"
60 #include <Galeri_MultiVectorTraits.hpp>
61 #include <Galeri_XpetraProblemFactory.hpp>
62 #include <Galeri_XpetraParameters.hpp>
63 #include <MatrixMarket_Tpetra.hpp>
79 template <
typename lno_t,
typename gno_t,
typename scalar_t,
typename nod_t>
81 const Teuchos::RCP<
const Teuchos::Comm<int> > &comm,
82 Teuchos::RCP<Tpetra::CrsMatrix<scalar_t, lno_t, gno_t, nod_t>> &M_)
85 if (comm->getRank() == 0){
86 std::cout <<
"Create matrix with " << problemType;
87 std::cout <<
" (a " << xdim;
89 std::cout <<
" x " << ydim <<
" x " << zdim <<
" ";
91 std::cout <<
" x" << ydim <<
" x 1 ";
93 std::cout <<
"x 1 x 1 ";
94 std::cout <<
" mesh)" << std::endl;
97 Teuchos::CommandLineProcessor tclp;
98 Galeri::Xpetra::Parameters<gno_t> params(tclp, xdim, ydim, zdim, problemType);
100 Teuchos::RCP<const Tpetra::Map<lno_t, gno_t> > map =
101 Teuchos::rcp(
new Tpetra::Map<lno_t, gno_t>(params.GetNumGlobalElements(), 0, comm));
105 Teuchos::RCP<Galeri::Xpetra::Problem<Tpetra::Map<lno_t, gno_t>,
106 Tpetra::CrsMatrix<scalar_t, lno_t, gno_t, nod_t>,
107 Tpetra::MultiVector<scalar_t, lno_t, gno_t> > > Pr=
108 Galeri::Xpetra::BuildProblem<scalar_t,
lno_t,
gno_t,
109 Tpetra::Map<lno_t, gno_t>,
110 Tpetra::CrsMatrix<scalar_t, lno_t, gno_t, nod_t>,
111 Tpetra::MultiVector<scalar_t, lno_t, gno_t, nod_t> >
112 (params.GetMatrixType(), map, params.GetParameterList());
114 M_ = Pr->BuildMatrix();
116 catch (std::exception &e) {
117 if(comm->getRank() == 0) std::cout <<
"Error returned from Galeri " << e.what() << std::endl;
126 template <
typename adapter_type>
132 typedef typename graph_type::global_ordinal_type GO;
133 typedef typename graph_type::local_ordinal_type LO;
134 typedef typename graph_type::node_type NO;
137 using ordinal_view_t = Kokkos::View<GO*, typename NO::device_type>;
138 using part_view_t = Kokkos::View<PT*, typename NO::device_type>;
140 auto graph = adapter->getUserGraph();
141 auto rowMap = graph->getRowMap();
142 auto colMap = graph->getColMap();
144 size_t numLclRows = rowMap->getLocalNumElements();
145 size_t numGblRows = rowMap->getGlobalNumElements();
146 size_t numLclCols = colMap->getLocalNumElements();
148 ordinal_view_t colLocalToGlobal(Kokkos::view_alloc(
"colLocalToGlobal", Kokkos::WithoutInitializing), numLclCols);
149 auto colMapHost = Kokkos::create_mirror_view (Kokkos::HostSpace (), colLocalToGlobal);
150 for(
size_t i = 0; i < numLclCols; ++i)
151 colMapHost[i] = colMap->getGlobalElement(i);
152 Kokkos::deep_copy (colLocalToGlobal, colMapHost);
154 ordinal_view_t rowLocalToGlobal(Kokkos::view_alloc(
"rowLocalToGlobal", Kokkos::WithoutInitializing), numLclRows);
155 auto rowMapHost = Kokkos::create_mirror_view (Kokkos::HostSpace (), rowLocalToGlobal);
156 for(
size_t i = 0; i < numLclRows; ++i)
157 rowMapHost[i] = rowMap->getGlobalElement(i);
158 Kokkos::deep_copy (rowLocalToGlobal, rowMapHost);
160 part_view_t localParts(
"localParts", numGblRows);
161 part_view_t globalParts(
"globalParts", numGblRows);
162 auto localPartsHost = Kokkos::create_mirror_view(Kokkos::HostSpace(), localParts);
165 for(
size_t i = 0; i < numLclRows; i++){
166 GO gi = rowMap->getGlobalElement(i);
167 localPartsHost(gi) = parts[i];
169 Kokkos::deep_copy(localParts, localPartsHost);
171 auto comm = graph->getComm();
172 Teuchos::reduceAll<int, PT> (*comm, Teuchos::REDUCE_SUM, numGblRows, localParts.data(), globalParts.data());
174 auto rowPtr = graph->getLocalGraphHost().row_map;
175 auto colInd = graph->getLocalGraphHost().entries;
177 size_t localtotalcut = 0, totalcut = 0;
179 using execution_space =
typename NO::device_type::execution_space;
180 using range_policy = Kokkos::RangePolicy<execution_space, Kokkos::IndexType<LO>>;
181 Kokkos::parallel_reduce(
"Compute cut", range_policy(0, numLclRows),
182 KOKKOS_LAMBDA(
const LO i,
size_t &cut){
184 const GO gRid = rowLocalToGlobal(i);
185 const PT gi = globalParts(gRid);
187 const size_t start = rowPtr(i);
188 const size_t end = rowPtr(i+1);
189 for(
size_t j = start; j < end; ++j) {
191 const GO gCid = colLocalToGlobal(colInd(j));
192 PT gj = globalParts(gCid);
198 Teuchos::reduceAll (*comm, Teuchos::REDUCE_SUM, 1, &localtotalcut, &totalcut);
201 auto rowPtr_h = Kokkos::create_mirror_view(rowPtr);
202 Kokkos::deep_copy(rowPtr_h, rowPtr);
205 size_t *partw =
new size_t[nparts];
206 size_t *partc =
new size_t[nparts];
208 size_t *gpartw =
new size_t[nparts];
209 size_t *gpartc =
new size_t[nparts];
211 for(
int i = 0; i < nparts; i++){
212 partw[i] = 0; partc[i] = 0;
213 gpartw[i] = 0; gpartc[i] = 0;
216 for(
size_t i = 0; i < numLclRows; i++){
217 partw[parts[i]] += rowPtr_h(i+1) - rowPtr_h(i) - 1;
221 Teuchos::reduceAll (*comm, Teuchos::REDUCE_SUM, nparts, partw, gpartw);
222 Teuchos::reduceAll (*comm, Teuchos::REDUCE_SUM, nparts, partc, gpartc);
224 size_t maxc = 0, totc = 0;
225 size_t maxw = 0, totw = 0;
227 for(
int i = 0; i < nparts; i++){
236 double imbw = (double)maxw/((
double)totw/nparts);
237 double imbc = (double)maxc/((
double)totc/nparts);
239 if(comm->getRank() == 0) {
241 std::cout <<
"\n\n************************************************" << std::endl;
242 std::cout <<
" EDGECUT: " << totalcut << std::endl;
243 std::cout <<
" MAX/AVG WEIGHT: " << imbw << std::endl;
244 std::cout <<
" MAX/AVG COUNT: " << imbc << std::endl;
245 std::cout <<
"************************************************\n\n" << std::endl;
251 int main(
int narg,
char *arg[])
254 Tpetra::ScopeGuard tpetraScope (&narg, &arg);
257 const Teuchos::RCP<const Teuchos::Comm<int>> pComm= Tpetra::getDefaultComm();
259 int me = pComm->getRank();
263 int max_iters = 1000;
268 std::string matrix_file =
"";
269 std::string eigensolve =
"LOBPCG";
270 bool parmetis =
false;
275 std::string ptype =
"";
276 std::string prec =
"jacobi";
277 std::string init =
"random";
282 for (
int i = 0; i < narg; i++)
283 std::cout << arg[i] <<
" ";
284 std::cout << std::endl;
287 Teuchos::CommandLineProcessor cmdp(
false,
true);
288 cmdp.setOption(
"matrix_file",&matrix_file,
289 "Path and filename of the matrix to be read.");
290 cmdp.setOption(
"nparts",&nparts,
291 "Number of global parts desired in the resulting partition.");
292 cmdp.setOption(
"rand_seed",&rand_seed,
293 "Seed for the random multivector.");
294 cmdp.setOption(
"max_iters",&max_iters,
295 "Maximum iters for the eigensolver.");
296 cmdp.setOption(
"block_size",&block_size,
298 cmdp.setOption(
"verbosity", &verbosity,
300 cmdp.setOption(
"parmetis",
"sphynx", &parmetis,
301 "Whether to use parmetis.");
302 cmdp.setOption(
"pulp",
"sphynx", &pulp,
303 "Whether to use pulp.");
304 cmdp.setOption(
"prec", &prec,
305 "Prec type to use.");
306 cmdp.setOption(
"eigensolve", &eigensolve,
307 "Eigensolver to use: LOBPCG, BlockDavidson, GeneralizedDavidson, BlockKrylovSchur or randomized.");
308 cmdp.setOption(
"prob", &ptype,
309 "Problem type to use. Options are combinatorial, normalized or generalized.");
310 cmdp.setOption(
"tol", &tol,
311 "Tolerance to use.");
312 cmdp.setOption(
"init", &init,
313 "Sphynx Initial guess. Options: random or constants. Default: random if randomized solver is used.");
314 cmdp.setOption(
"resFreq", &resFreq,
315 "(For randomized) Specify how often to check the residuals. Orthogonalization of the basis is also done.");
316 cmdp.setOption(
"orthoFreq", &orthoFreq,
317 "(For randomized) Specify how often to orthogonalize the basis.");
319 if (cmdp.parse(narg,arg)!=Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL)
return -1;
323 std::cout << std::endl <<
"--------------------------------------------------" << matrix_file << std::endl;
324 std::cout <<
"| Sphynx Parameters |" << matrix_file << std::endl;
325 std::cout <<
"--------------------------------------------------" << matrix_file << std::endl;
326 std::cout <<
" matrix file = " << matrix_file << std::endl;
327 std::cout <<
" nparts = " << nparts << std::endl;
328 std::cout <<
" verbosity = " << verbosity << std::endl;
329 std::cout <<
" parmetis = " << parmetis << std::endl;
330 std::cout <<
" pulp = " << pulp << std::endl;
331 std::cout <<
" prec = " << prec << std::endl;
332 std::cout <<
" eigensolver = " << eigensolve << std::endl;
333 std::cout <<
" prob = " << ptype << std::endl;
334 std::cout <<
" tol = " << tol << std::endl;
335 std::cout <<
" init = " << init << std::endl;
336 std::cout <<
" resFreq = " << resFreq << std::endl;
337 std::cout <<
" orthoFreq = " << orthoFreq << std::endl;
338 std::cout <<
"--------------------------------------------------" << matrix_file << std::endl << std::endl;
341 using scalar_type = Tpetra::Details::DefaultTypes::scalar_type;
342 using local_ordinal_type = Tpetra::Details::DefaultTypes::local_ordinal_type;
343 using global_ordinal_type = Tpetra::Details::DefaultTypes::global_ordinal_type;
344 using node_type = Tpetra::Details::DefaultTypes::node_type;
346 using crs_matrix_type = Tpetra::CrsMatrix<scalar_type, local_ordinal_type, global_ordinal_type, node_type>;
351 Teuchos::RCP<adapter_type> adapter;
352 Teuchos::RCP<crs_matrix_type> tmatrix;
355 std::srand(rand_seed);
358 std::string mtx =
".mtx", mm =
".mm", lc =
".largestComp", lc2 =
".bin";
359 if(std::equal(lc.rbegin(), lc.rend(), matrix_file.rbegin()) || std::equal(lc2.rbegin(), lc2.rend(), matrix_file.rbegin())) {
360 tmatrix = readMatrixFromBinaryFile<crs_matrix_type>(matrix_file, pComm,
true, verbosity>0);
361 if(me == 0 && verbosity > 1) std::cout <<
"Used Seher's reader for Largest Comp." << std::endl;
363 else if(std::equal(mtx.rbegin(), mtx.rend(), matrix_file.rbegin()) || std::equal(mm.rbegin(), mm.rend(), matrix_file.rbegin())) {
364 typedef Tpetra::MatrixMarket::Reader<crs_matrix_type> reader_type;
366 tmatrix = r.readSparseFile(matrix_file, pComm);
367 if(me == 0 && verbosity > 1) std::cout <<
"Used standard Matrix Market reader." << std::endl;
374 if(matrix_file !=
"")
376 std::string::const_iterator it = matrix_file.begin();
377 while(it != matrix_file.end() && std::isdigit(*it)) ++it;
378 if(it == matrix_file.end())
380 meshdim = std::stoi(matrix_file);
383 if(me == 0) std::cout <<
"Invalid matrix file entered. Reverting to default matrix." << std::endl;
386 int ierr = buildCrsMatrix<local_ordinal_type, global_ordinal_type, scalar_type>
387 (meshdim, meshdim, meshdim,
"Brick3D", pComm, tmatrix);
388 if (ierr != 0) std::cout <<
"Error! Brick3D failed to build!" << std::endl;
391 if(me == 0 && verbosity > 0) std::cout <<
"Done with reading/creating the matrix." << std::endl;
393 Teuchos::RCP<const Tpetra::Map<> > map = tmatrix->getMap();
395 adapter = Teuchos::rcp(
new adapter_type(tmatrix->getCrsGraph(), 1));
396 adapter->setVertexWeightIsDegree(0);
399 Teuchos::RCP<Teuchos::ParameterList> params = Teuchos::rcp(
new Teuchos::ParameterList());
400 Teuchos::RCP<Teuchos::ParameterList> sphynxParams(
new Teuchos::ParameterList);
401 params->set(
"num_global_parts", nparts);
403 Teuchos::RCP<Teuchos::StackedTimer> stacked_timer;
404 stacked_timer = Teuchos::rcp(
new Teuchos::StackedTimer(
"SphynxDriver"));
405 Teuchos::TimeMonitor::setStackedTimer(stacked_timer);
406 if(parmetis || pulp) {
408 params->set(
"partitioning_approach",
"partition");
409 params->set(
"imbalance_tolerance", 1.01);
411 params->set(
"algorithm",
"parmetis");
412 params->set(
"imbalance_tolerance", 1.01);
415 params->set(
"algorithm",
"pulp");
416 params->set(
"pulp_vert_imbalance", 1.01);
420 Teuchos::RCP<problem_type> problem;
423 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer(
"Partitioning::All"));
425 Teuchos::TimeMonitor t2(*Teuchos::TimeMonitor::getNewTimer(
"Partitioning::Problem"));
426 problem = Teuchos::rcp(
new problem_type(adapter.getRawPtr(), params.getRawPtr(), Tpetra::getDefaultComm()));
429 Teuchos::TimeMonitor t3(*Teuchos::TimeMonitor::getNewTimer(
"Partitioning::Solve"));
435 solution_type solution = problem->getSolution();
436 compute_edgecut<adapter_type>(adapter, solution);
440 sphynxParams->set(
"sphynx_verbosity", verbosity);
441 sphynxParams->set(
"sphynx_max_iterations", max_iters);
442 sphynxParams->set(
"sphynx_ortho_freq", orthoFreq);
443 sphynxParams->set(
"sphynx_res_freq", resFreq);
444 sphynxParams->set(
"sphynx_skip_preprocessing",
true);
445 sphynxParams->set(
"sphynx_eigensolver", eigensolve);
446 if (block_size > 0) sphynxParams->set(
"sphynx_block_size", block_size);
447 if (ptype !=
"") sphynxParams->set(
"sphynx_problem_type", ptype);
448 if (init !=
"") sphynxParams->set(
"sphynx_initial_guess", init);
449 if (prec !=
"") sphynxParams->set(
"sphynx_preconditioner_type", prec);
450 if (tol != -1) sphynxParams->set(
"sphynx_tolerance", tol);
453 Teuchos::RCP<problem_type> problem;
456 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer(
"Partitioning::All"));
458 Teuchos::TimeMonitor t2(*Teuchos::TimeMonitor::getNewTimer(
"Partitioning::Problem"));
459 problem = Teuchos::rcp(
new problem_type(adapter.get(), params.get(), sphynxParams, Tpetra::getDefaultComm()));
462 if(me == 0) std::cout << eigensolve <<
" will be used to solve the partitioning problem." << std::endl;
467 solution_type solution = problem->getSolution();
468 compute_edgecut<adapter_type>(adapter, solution);
470 stacked_timer->stopBaseTimer();
472 Teuchos::RCP<Teuchos::FancyOStream> fancy2 = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
473 Teuchos::FancyOStream& out2 = *fancy2;
474 Teuchos::StackedTimer::OutputOptions options;
475 options.output_fraction = options.output_histogram = options.output_minmax =
true;
476 stacked_timer->report(out2, pComm, options);
478 Teuchos::TimeMonitor::summarize();
void compute_edgecut(Teuchos::RCP< adapter_type > &adapter, Zoltan2::PartitioningSolution< adapter_type > &solution)
map_t::global_ordinal_type gno_t
Provides access for Zoltan2 to Xpetra::CrsGraph data.
int main(int narg, char **arg)
int buildCrsMatrix(int xdim, int ydim, int zdim, std::string problemType, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, Teuchos::RCP< Tpetra::CrsMatrix< scalar_t, lno_t, gno_t, nod_t >> &M_)
Defines XpetraCrsGraphAdapter class.
SparseMatrixAdapter_t::part_t part_t
A PartitioningSolution is a solution to a partitioning problem.
const part_t * getPartListView() const
Returns the part list corresponding to the global ID list.
map_t::local_ordinal_type lno_t
PartitioningProblem sets up partitioning problems for the user.
Defines the PartitioningProblem class.
size_t getTargetGlobalNumberOfParts() const
Returns the global number of parts desired in the solution.
Zoltan2::BasicUserTypes< zscalar_t, zlno_t, zgno_t > user_t