106 using lno_t =
typename Adapter::lno_t;
107 using gno_t =
typename Adapter::gno_t;
113 using graph_t = Tpetra::CrsGraph<lno_t, gno_t, node_t>;
114 using matrix_t = Tpetra::CrsMatrix<scalar_t, lno_t, gno_t, node_t>;
115 using mvector_t = Tpetra::MultiVector<scalar_t, lno_t, gno_t, node_t>;
116 using op_t = Tpetra::Operator<scalar_t, lno_t, gno_t, node_t>;
125 Sphynx(
const RCP<const Environment> &env,
126 const RCP<Teuchos::ParameterList> ¶ms,
127 const RCP<Teuchos::ParameterList> &sphynxParams,
128 const RCP<
const Comm<int> > &comm,
129 const RCP<
const XpetraCrsGraphAdapter<graph_t> > &adapter) :
130 env_(env), params_(params), sphynxParams_(sphynxParams), comm_(comm), adapter_(adapter)
134 const Teuchos::ParameterEntry *pe = params_->getEntryPtr(
"num_global_parts");
135 numGlobalParts_ = pe->getValue<
int>(&numGlobalParts_);
136 if(numGlobalParts_ > 1){
138 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer(
"Sphynx::Laplacian"));
141 pe = sphynxParams_->getEntryPtr(
"sphynx_verbosity");
143 verbosity_ = pe->getValue<
int>(&verbosity_);
146 pe = sphynxParams_->getEntryPtr(
"sphynx_skip_preprocessing");
148 skipPrep_ = pe->getValue<
bool>(&skipPrep_);
153 graph_ = adapter_->getUserGraph();
155 throw std::runtime_error(
"\nSphynx Error: Preprocessing has not been implemented yet.\n");
177 void partition(
const RCP<PartitioningSolution<Adapter> > &solution);
182 template<
typename problem_t>
185 template<
typename problem_t>
188 template<
typename problem_t>
191 template<
typename problem_t>
201 std::vector<int> strides);
205 std::vector<const weight_t *>
weights,
206 std::vector<int> strides,
207 const Teuchos::RCP<PartitioningSolution<Adapter>> &solution);
222 auto rowOffsets = graph_->getLocalGraphDevice().row_map;
223 auto rowOffsets_h = Kokkos::create_mirror_view(rowOffsets);
224 Kokkos::deep_copy(rowOffsets_h, rowOffsets);
227 const size_t numGlobalEntries = graph_->getGlobalNumEntries();
228 const size_t numLocalRows = graph_->getLocalNumRows();
229 const size_t numGlobalRows = graph_->getGlobalNumRows();
233 for(
size_t i = 0; i < numLocalRows; i++){
234 if(rowOffsets_h(i+1) - rowOffsets_h(i) - 1 > localMax)
235 localMax = rowOffsets_h(i+1) - rowOffsets_h(i) - 1;
239 size_t globalMax = 0;
240 Teuchos::reduceAll<int, size_t> (*comm_, Teuchos::REDUCE_MAX, 1, &localMax, &globalMax);
242 double avg =
static_cast<double>(numGlobalEntries-numGlobalRows)/numGlobalRows;
243 double maxOverAvg =
static_cast<double>(globalMax)/avg;
247 if(maxOverAvg > 10) {
253 if(comm_->getRank() == 0) {
254 std::cout <<
"Determining Regularity -- max degree: " << globalMax
255 <<
", avg degree: " << avg <<
", max/avg: " << globalMax/avg << std::endl
256 <<
"Determined Regularity -- regular: " << !irregular_ << std::endl;
277 precType_ =
"jacobi";
278#ifdef HAVE_ZOLTAN2SPHYNX_MUELU
283 const Teuchos::ParameterEntry *pe = sphynxParams_->getEntryPtr(
"sphynx_preconditioner_type");
285 precType_ = pe->getValue<std::string>(&precType_);
286 if(precType_ !=
"muelu" && precType_ !=
"jacobi" && precType_ !=
"polynomial")
287 throw std::runtime_error(
"\nSphynx Error: " + precType_ +
" is not recognized as a preconditioner.\n"
288 +
" Possible values: muelu (if enabled), jacobi, and polynomial\n");
294 if(precType_ ==
"polynomial")
301 pe = sphynxParams_->getEntryPtr(
"sphynx_problem_type");
303 std::string probType =
"";
304 probType = pe->getValue<std::string>(&probType);
306 if(probType ==
"combinatorial")
308 else if(probType ==
"generalized")
310 else if(probType ==
"normalized")
313 throw std::runtime_error(
"\nSphynx Error: " + probType +
" is not recognized as a problem type.\n"
314 +
" Possible values: combinatorial, generalized, and normalized.\n");
320 if(!irregular_ && (precType_ ==
"jacobi" || precType_ ==
"polynomial"))
325 pe = sphynxParams_->getEntryPtr(
"sphynx_tolerance");
327 tolerance_ = pe->getValue<
scalar_t>(&tolerance_);
336 pe = sphynxParams_->getEntryPtr(
"sphynx_initial_guess");
338 std::string initialGuess =
"";
339 initialGuess = pe->getValue<std::string>(&initialGuess);
341 if(initialGuess ==
"random")
343 else if(initialGuess ==
"constants")
346 throw std::runtime_error(
"\nSphynx Error: " + initialGuess +
" is not recognized as an option for initial guess.\n"
347 +
" Possible values: random and constants.\n");
372 auto rowOffsets = graph_->getLocalGraphHost().row_map;
375 Teuchos::RCP<matrix_t> degMat(
new matrix_t (graph_->getRowMap(),
381 lno_t numRows =
static_cast<lno_t>(graph_->getLocalNumRows());
384 for (
lno_t i = 0; i < numRows; ++i) {
385 val[0] = rowOffsets(i+1) - rowOffsets(i) - 1;
387 degMat->insertLocalValues(i, 1, val, ind);
392 degMat->fillComplete(graph_->getDomainMap(), graph_->getRangeMap());
404 using range_policy = Kokkos::RangePolicy<
405 typename node_t::device_type::execution_space, Kokkos::IndexType<lno_t>>;
406 using values_view_t = Kokkos::View<scalar_t*, typename node_t::device_type>;
407 using offset_view_t = Kokkos::View<size_t*, typename node_t::device_type>;
409 const size_t numEnt = graph_->getLocalNumEntries();
410 const size_t numRows = graph_->getLocalNumRows();
413 values_view_t newVal (Kokkos::view_alloc(
"CombLapl::val", Kokkos::WithoutInitializing), numEnt);
414 Kokkos::deep_copy(newVal, -1);
417 offset_view_t diagOffsets(Kokkos::view_alloc(
"Diag Offsets", Kokkos::WithoutInitializing), numRows);
418 graph_->getLocalDiagOffsets(diagOffsets);
421 auto rowOffsets = graph_->getLocalGraphDevice().row_map;
424 Kokkos::parallel_for(
"Combinatorial Laplacian", range_policy(0, numRows),
425 KOKKOS_LAMBDA(
const lno_t i){
426 newVal(rowOffsets(i) + diagOffsets(i)) = rowOffsets(i+1) - rowOffsets(i) - 1;
432 Teuchos::RCP<matrix_t> laplacian (
new matrix_t(graph_, newVal));
433 laplacian->fillComplete (graph_->getDomainMap(), graph_->getRangeMap());
447 using range_policy = Kokkos::RangePolicy<
448 typename node_t::device_type::execution_space, Kokkos::IndexType<lno_t>>;
449 using values_view_t = Kokkos::View<scalar_t*, typename node_t::device_type>;
450 using offset_view_t = Kokkos::View<size_t*, typename node_t::device_type>;
451 using vector_t = Tpetra::Vector<scalar_t, lno_t, gno_t, node_t>;
452 using dual_view_t =
typename vector_t::dual_view_type;
453 using KAT = Kokkos::Details::ArithTraits<scalar_t>;
455 const size_t numEnt = graph_->getLocalNumEntries();
456 const size_t numRows = graph_->getLocalNumRows();
459 values_view_t newVal (Kokkos::view_alloc(
"NormLapl::val", Kokkos::WithoutInitializing), numEnt);
460 Kokkos::deep_copy(newVal, -1);
463 dual_view_t dv (
"MV::DualView", numRows, 1);
464 auto deginvsqrt = dv.d_view;
467 offset_view_t diagOffsets(Kokkos::view_alloc(
"Diag Offsets", Kokkos::WithoutInitializing), numRows);
468 graph_->getLocalDiagOffsets(diagOffsets);
471 auto rowOffsets = graph_->getLocalGraphDevice().row_map;
474 Kokkos::parallel_for(
"Combinatorial Laplacian", range_policy(0, numRows),
475 KOKKOS_LAMBDA(
const lno_t i){
476 newVal(rowOffsets(i) + diagOffsets(i)) = rowOffsets(i+1) - rowOffsets(i) - 1;
477 deginvsqrt(i,0) = 1.0/KAT::sqrt(rowOffsets(i+1) - rowOffsets(i) - 1);
483 Teuchos::RCP<matrix_t> laplacian (
new matrix_t(graph_, newVal));
484 laplacian->fillComplete (graph_->getDomainMap(), graph_->getRangeMap());
487 vector_t degInvSqrt(graph_->getRowMap(), dv);
490 laplacian->leftScale(degInvSqrt);
491 laplacian->rightScale(degInvSqrt);
503 const Teuchos::RCP<const Environment> env_;
504 Teuchos::RCP<Teuchos::ParameterList> params_;
505 Teuchos::RCP<Teuchos::ParameterList> sphynxParams_;
506 const Teuchos::RCP<const Teuchos::Comm<int> > comm_;
507 const Teuchos::RCP<const Adapter> adapter_;
511 Teuchos::RCP<const graph_t> graph_;
512 Teuchos::RCP<matrix_t> laplacian_;
513 Teuchos::RCP<const matrix_t> degMatrix_;
514 Teuchos::RCP<mvector_t> eigenVectors_;
517 std::string precType_;
541 if(numGlobalParts_ == 1) {
543 size_t numRows =adapter_->getUserGraph()->getLocalNumRows();
544 Teuchos::ArrayRCP<part_t> parts(numRows,0);
545 solution->setParts(parts);
552 int numEigenVectors = (int) log2(numGlobalParts_)+1;
557 if(computedNumEv <= 1) {
559 std::runtime_error(
"\nAnasazi Error: LOBPCGSolMgr::solve() returned unconverged.\n"
560 "Sphynx Error: LOBPCG could not compute any eigenvectors.\n"
561 " Increase either max iters or tolerance.\n");
570 std::vector<const weight_t *>
weights;
571 std::vector<int> wstrides;
587 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer(
"Sphynx::LOBPCG"));
590 std::string which =
"SR";
591 std::string ortho =
"SVQB";
592 bool relConvTol =
false;
593 int maxIterations = 1000;
594 bool isHermitian =
true;
595 bool relLockTol =
false;
597 bool useFullOrtho =
true;
602 double solvetime = 0;
605 const Teuchos::ParameterEntry *pe;
607 pe = sphynxParams_->getEntryPtr(
"sphynx_maxiterations");
609 maxIterations = pe->getValue<
int>(&maxIterations);
611 pe = sphynxParams_->getEntryPtr(
"sphynx_use_full_ortho");
613 useFullOrtho = pe->getValue<
bool>(&useFullOrtho);
617 int anasaziVerbosity = Anasazi::Errors + Anasazi::Warnings;
619 anasaziVerbosity += Anasazi::FinalSummary + Anasazi::TimingDetails;
621 anasaziVerbosity += Anasazi::IterationDetails;
623 anasaziVerbosity += Anasazi::StatusTestDetails
624 + Anasazi::OrthoDetails
629 Teuchos::ParameterList anasaziParams;
630 anasaziParams.set(
"Verbosity", anasaziVerbosity);
631 anasaziParams.set(
"Which", which);
632 anasaziParams.set(
"Convergence Tolerance", tolerance_);
633 anasaziParams.set(
"Maximum Iterations", maxIterations);
634 anasaziParams.set(
"Block Size", numEigenVectors);
635 anasaziParams.set(
"Relative Convergence Tolerance", relConvTol);
636 anasaziParams.set(
"Orthogonalization", ortho);
637 anasaziParams.set(
"Use Locking", lock);
638 anasaziParams.set(
"Relative Locking Tolerance", relLockTol);
639 anasaziParams.set(
"Full Ortho", useFullOrtho);
642 RCP<mvector_t> ivec(
new mvector_t(laplacian_->getRangeMap(), numEigenVectors));
647 Anasazi::MultiVecTraits<scalar_t, mvector_t>::MvRandom(*ivec);
648 ivec->getVectorNonConst(0)->putScalar(1.);
656 ivec->getVectorNonConst(0)->putScalar(1.);
657 for (
int j = 1; j < numEigenVectors; j++)
658 ivec->getVectorNonConst(j)->putScalar(0.);
660 auto map = laplacian_->getRangeMap();
661 gno_t blockSize = map->getGlobalNumElements() / numEigenVectors;
662 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::runtime_error,
"Blocksize too small for \"constants\" initial guess. Try \"random\".");
664 for (
size_t lid = 0; lid < ivec->getLocalLength(); lid++) {
665 gno_t gid = map->getGlobalElement(lid);
666 for (
int j = 1; j < numEigenVectors; j++){
667 if (((j-1)*blockSize <= gid) && (j*blockSize > gid))
668 ivec->replaceLocalValue(lid,j,1.);
674 using problem_t = Anasazi::BasicEigenproblem<scalar_t, mvector_t, op_t>;
675 Teuchos::RCP<problem_t> problem (
new problem_t(laplacian_, ivec));
676 problem->setHermitian(isHermitian);
677 problem->setNEV(numEigenVectors);
684 problem->setM(degMatrix_);
687 bool boolret = problem->setProblem();
688 if (boolret !=
true) {
689 throw std::runtime_error(
"\nAnasazi::BasicEigenproblem::setProblem() returned with error.\n");
693 using solver_t = Anasazi::LOBPCGSolMgr<scalar_t, mvector_t, op_t>;
694 solver_t solver(problem, anasaziParams);
696 if (verbosity_ > 0 && comm_->getRank() == 0)
697 anasaziParams.print(std::cout);
700 if (verbosity_ > 0 && comm_->getRank() == 0)
701 std::cout <<
"Beginning the LOBPCG solve..." << std::endl;
702 Anasazi::ReturnType returnCode = solver.solve();
705 if (returnCode != Anasazi::Converged) {
708 iter = solver.getNumIters();
709 solvetime = (solver.getTimers()[0])->totalElapsedTime();
713 using solution_t = Anasazi::Eigensolution<scalar_t, mvector_t>;
714 solution_t sol = problem->getSolution();
715 eigenVectors_ = sol.Evecs;
716 int numev = sol.numVecs;
719 if (verbosity_ > 0 && comm_->getRank() == 0) {
720 std::cout << std::endl;
721 std::cout <<
"LOBPCG SUMMARY" << std::endl;
722 std::cout <<
"Failed to converge: " << numfailed << std::endl;
723 std::cout <<
"No of iterations : " << iter << std::endl;
724 std::cout <<
"Solve time: " << solvetime << std::endl;
725 std::cout <<
"No of comp. vecs. : " << numev << std::endl;
759#ifdef HAVE_ZOLTAN2SPHYNX_MUELU
760 Teuchos::ParameterList paramList;
762 paramList.set(
"verbosity",
"none");
763 else if(verbosity_ == 1)
764 paramList.set(
"verbosity",
"low");
765 else if(verbosity_ == 2)
766 paramList.set(
"verbosity",
"medium");
767 else if(verbosity_ == 3)
768 paramList.set(
"verbosity",
"high");
769 else if(verbosity_ >= 4)
770 paramList.set(
"verbosity",
"extreme");
772 paramList.set(
"smoother: type",
"CHEBYSHEV");
773 Teuchos::ParameterList smootherParamList;
774 smootherParamList.set(
"chebyshev: degree", 3);
775 smootherParamList.set(
"chebyshev: ratio eigenvalue", 7.0);
776 smootherParamList.set(
"chebyshev: eigenvalue max iterations", irregular_ ? 100 : 10);
777 paramList.set(
"smoother: params", smootherParamList);
778 paramList.set(
"use kokkos refactor",
true);
779 paramList.set(
"transpose: use implicit",
true);
783 paramList.set(
"multigrid algorithm",
"unsmoothed");
785 paramList.set(
"coarse: type",
"CHEBYSHEV");
786 Teuchos::ParameterList coarseParamList;
787 coarseParamList.set(
"chebyshev: degree", 3);
788 coarseParamList.set(
"chebyshev: ratio eigenvalue", 7.0);
789 paramList.set(
"coarse: params", coarseParamList);
791 paramList.set(
"max levels", 5);
792 paramList.set(
"aggregation: drop tol", 0.40);
796 using prec_t = MueLu::TpetraOperator<scalar_t, lno_t, gno_t, node_t>;
797 Teuchos::RCP<prec_t> prec = MueLu::CreateTpetraPreconditioner<
800 problem->setPrec(prec);
803 throw std::runtime_error(
"\nSphynx Error: MueLu requested but not compiled into Trilinos.\n");
813 int verbosity2 = Belos::Errors;
815 verbosity2 += Belos::Warnings;
816 else if(verbosity_ == 2)
817 verbosity2 += Belos::Warnings + Belos::FinalSummary;
818 else if(verbosity_ == 3)
819 verbosity2 += Belos::Warnings + Belos::FinalSummary + Belos::TimingDetails;
820 else if(verbosity_ >= 4)
821 verbosity2 += Belos::Warnings + Belos::FinalSummary + Belos::TimingDetails
822 + Belos::StatusTestDetails;
824 Teuchos::ParameterList paramList;
825 paramList.set(
"Polynomial Type",
"Roots");
826 paramList.set(
"Orthogonalization",
"ICGS");
827 paramList.set(
"Maximum Degree", laplacian_->getGlobalNumRows() > 100 ? 25 : 5);
828 paramList.set(
"Polynomial Tolerance", 1.0e-6 );
829 paramList.set(
"Verbosity", verbosity2 );
830 paramList.set(
"Random RHS",
false );
831 paramList.set(
"Outer Solver",
"");
832 paramList.set(
"Timer Label",
"Belos Polynomial Solve" );
835 using lproblem_t = Belos::LinearProblem<scalar_t, mvector_t, op_t>;
836 Teuchos::RCP<lproblem_t> innerPolyProblem(
new lproblem_t());
837 innerPolyProblem->setOperator(laplacian_);
839 using btop_t = Belos::TpetraOperator<scalar_t>;
840 Teuchos::RCP<btop_t> polySolver(
new btop_t(innerPolyProblem,
841 Teuchos::rcpFromRef(paramList),
843 problem->setPrec(polySolver);
894 std::vector<int> strides)
897 int numWeights = adapter_->getNumWeightsPerVertex();
898 int numConstraints = numWeights > 0 ? numWeights : 1;
900 size_t myNumVertices = adapter_->getLocalNumVertices();
902 for(
int j = 0; j < numConstraints; j++)
907 const gno_t *columnId;
910 if(numWeights == 0) {
913 adapter_->getEdgesView(offset, columnId);
914 for (
size_t i = 0; i < myNumVertices; i++)
915 weights[0][i] = offset[i+1] - offset[i] - 1;
917 vecweights.push_back(
weights[0]);
918 strides.push_back(1);
923 for(
int j = 0; j < numConstraints; j++) {
925 if(adapter_->useDegreeAsVertexWeight(j)) {
927 adapter_->getEdgesView(offset, columnId);
928 for (
size_t i = 0; i < myNumVertices; i++)
929 weights[j][i] = offset[i+1] - offset[i];
934 adapter_->getVertexWeightsView(wgt, stride, j);
936 for (
size_t i = 0; i < myNumVertices; i++)
940 vecweights.push_back(
weights[j]);
941 strides.push_back(1);
953 std::vector<const weight_t *>
weights,
954 std::vector<int> strides,
955 const Teuchos::RCP<PartitioningSolution<Adapter>> &solution)
958 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer(
"Sphynx::MJ"));
961 using base_adapter_t =
typename mvector_adapter_t::base_adapter_t;
967 size_t myNumVertices =
coordinates->getLocalLength();
970 Teuchos::RCP<mvector_adapter_t> adapcoordinates(
new mvector_adapter_t(
coordinates,
973 Teuchos::RCP<const base_adapter_t> baseAdapter =
974 Teuchos::rcp(
dynamic_cast<const base_adapter_t *
>(adapcoordinates.get()),
false);
981 Teuchos::RCP<const Comm<int>> comm2 = comm_;
982 Teuchos::RCP<mj_t> mj(
new mj_t(env_, comm2, baseAdapter));
985 Teuchos::RCP<solution_t> vectorsolution(
new solution_t(env_, comm2, 1, mj));
986 mj->partition(vectorsolution);
989 Teuchos::ArrayRCP<part_t> parts(myNumVertices);
990 for(
size_t i = 0; i < myNumVertices; i++) {
991 parts[i] = vectorsolution->getPartListView()[i];
993 solution->setParts(parts);