Zoltan2
Loading...
Searching...
No Matches
Zoltan2_Sphynx.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Sphynx
6// Copyright 2020 National Technology & Engineering
7// Solutions of Sandia, LLC (NTESS).
8//
9// Under the terms of Contract DE-NA0003525 with NTESS,
10// the U.S. Government retains certain rights in this software.
11//
12// Redistribution and use in source and binary forms, with or without
13// modification, are permitted provided that the following conditions are
14// met:
15//
16// 1. Redistributions of source code must retain the above copyright
17// notice, this list of conditions and the following disclaimer.
18//
19// 2. Redistributions in binary form must reproduce the above copyright
20// notice, this list of conditions and the following disclaimer in the
21// documentation and/or other materials provided with the distribution.
22//
23// 3. Neither the name of the Corporation nor the names of the
24// contributors may be used to endorse or promote products derived from
25// this software without specific prior written permission.
26//
27// THIS SOFTWARE IS PROVIDED BY NTESS "AS IS" AND ANY
28// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
29// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
30// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL NTESS OR THE
31// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
32// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
33// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
34// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
35// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
36// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
37// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38//
39// Questions? Contact Seher Acer (sacer@sandia.gov)
40// Erik Boman (egboman@sandia.gov)
41// Siva Rajamanickam (srajama@sandia.gov)
42// Karen Devine (kddevin@sandia.gov)
43//
44// ***********************************************************************
45//
46// @HEADER
47
48#ifndef _ZOLTAN2_SPHYNXALGORITHM_HPP_
49#define _ZOLTAN2_SPHYNXALGORITHM_HPP_
50
51
53// This file contains the Sphynx algorithm.
54//
55// Sphynx is a graph partitioning algorithm that is based on a spectral method.
56// It has three major steps:
57// 1) compute the Laplacian matrix of the input graph,
58// 2) compute logK+1 eigenvectors on the Laplacian matrix,
59// 3) use eigenvector entries as coordinates and compute a K-way partition on
60// them using a geometric method.
61//
62// Step1: Sphynx provides three eigenvalue problems and hence Laplacian matrix:
63// i) combinatorial (Lx = \lambdax, where L = A - D)
64// ii) generalized (Lx = \lambdaDx, where L = A - D)
65// iii) normalized (L_nx, \lambdax, where Ln = D^{-1/2}LD^{-1/2}
66// and L = A - D)
67//
68// Step2: Sphynx calls the LOBPCG algorithm provided in Anasazi to obtain
69// logK+1 eigenvectors.
70// Step3: Sphynx calls the MJ algorithm provided in Zoltan2Core to compute the
71// partition.
73
74#include "Zoltan2Sphynx_config.h"
75
78
81
82#include "AnasaziLOBPCGSolMgr.hpp"
83#include "AnasaziBasicEigenproblem.hpp"
84#include "AnasaziTpetraAdapter.hpp"
85
86#include "BelosLinearProblem.hpp"
87#include "BelosTpetraOperator.hpp"
88
89#include "Ifpack2_Factory.hpp"
90
91#include "Teuchos_TimeMonitor.hpp"
92
93#ifdef HAVE_ZOLTAN2SPHYNX_MUELU
94#include "MueLu_CreateTpetraPreconditioner.hpp"
95#endif
96
97namespace Zoltan2 {
98
99 template <typename Adapter>
100 class Sphynx : public Algorithm<Adapter>
101 {
102
103 public:
104
105 using scalar_t = double; // Sphynx with scalar_t=double obtains better cutsize
106 using lno_t = typename Adapter::lno_t;
107 using gno_t = typename Adapter::gno_t;
108 using node_t = typename Adapter::node_t;
109 using offset_t = typename Adapter::offset_t;
110 using part_t = typename Adapter::part_t;
111 using weight_t = typename Adapter::scalar_t;
112
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>;
117
119
123
124 // Takes the graph from the input adapter and computes the Laplacian matrix
125 Sphynx(const RCP<const Environment> &env,
126 const RCP<Teuchos::ParameterList> &params,
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)
131 {
132
133 // Don't compute the Laplacian if the number of requested parts is 1
134 const Teuchos::ParameterEntry *pe = params_->getEntryPtr("num_global_parts");
135 numGlobalParts_ = pe->getValue<int>(&numGlobalParts_);
136 if(numGlobalParts_ > 1){
137
138 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer("Sphynx::Laplacian"));
139
140 // The verbosity is common throughout the algorithm, better to check and set now
141 pe = sphynxParams_->getEntryPtr("sphynx_verbosity");
142 if (pe)
143 verbosity_ = pe->getValue<int>(&verbosity_);
144
145 // Do we need to pre-process the input?
146 pe = sphynxParams_->getEntryPtr("sphynx_skip_preprocessing");
147 if (pe)
148 skipPrep_ = pe->getValue<bool>(&skipPrep_);
149
150 // Get the graph from XpetraCrsGraphAdapter if skipPrep_ is true
151 // We assume the graph has all of the symmetric and diagonal entries
152 if(skipPrep_)
153 graph_ = adapter_->getUserGraph();
154 else {
155 throw std::runtime_error("\nSphynx Error: Preprocessing has not been implemented yet.\n");
156 }
157
158 // Check if the graph is regular
160
161 // Set Sphynx defaults: preconditioner, problem type, tolerance, initial vectors.
162 setDefaults();
163
164 // Compute the Laplacian matrix
166
167 if(problemType_ == GENERALIZED)
169
170 }
171 }
172
176
177 void partition(const RCP<PartitioningSolution<Adapter> > &solution);
178
179
180 int LOBPCGwrapper(const int numEigenVectors);
181
182 template<typename problem_t>
183 void setPreconditioner(Teuchos::RCP<problem_t> &problem);
184
185 template<typename problem_t>
186 void setMueLuPreconditioner(Teuchos::RCP<problem_t> &problem);
187
188 template<typename problem_t>
189 void setJacobiPreconditioner(Teuchos::RCP<problem_t> &problem);
190
191 template<typename problem_t>
192 void setPolynomialPreconditioner(Teuchos::RCP<problem_t> &problem);
193
194
195 void eigenvecsToCoords(Teuchos::RCP<mvector_t> &eigenVectors,
196 int computedNumEv,
197 Teuchos::RCP<mvector_t> &coordinates);
198
199
200 void computeWeights(std::vector<const weight_t *> vecweights,
201 std::vector<int> strides);
202
203
204 void MJwrapper(const Teuchos::RCP<const mvector_t> &coordinates,
205 std::vector<const weight_t *> weights,
206 std::vector<int> strides,
207 const Teuchos::RCP<PartitioningSolution<Adapter>> &solution);
208
209
213
215 // Determine input graph's regularity = maxDegree/AvgDegree < 10.
216 // If Laplacian type is not specified, then use combinatorial for regular
217 // and generalized for irregular.
218 // MueLu settings will differ depending on the regularity, too.
220 {
221 // Get the row pointers in the host
222 auto rowOffsets = graph_->getLocalGraphDevice().row_map;
223 auto rowOffsets_h = Kokkos::create_mirror_view(rowOffsets);
224 Kokkos::deep_copy(rowOffsets_h, rowOffsets);
225
226 // Get size information
227 const size_t numGlobalEntries = graph_->getGlobalNumEntries();
228 const size_t numLocalRows = graph_->getLocalNumRows();
229 const size_t numGlobalRows = graph_->getGlobalNumRows();
230
231 // Compute local maximum degree
232 size_t localMax = 0;
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;
236 }
237
238 // Compute global maximum degree
239 size_t globalMax = 0;
240 Teuchos::reduceAll<int, size_t> (*comm_, Teuchos::REDUCE_MAX, 1, &localMax, &globalMax);
241
242 double avg = static_cast<double>(numGlobalEntries-numGlobalRows)/numGlobalRows;
243 double maxOverAvg = static_cast<double>(globalMax)/avg;
244
245 // Use generalized Laplacian on irregular graphs
246 irregular_ = false;
247 if(maxOverAvg > 10) {
248 irregular_ = true;
249 }
250
251 // Let the user know about what we determined if verbose
252 if(verbosity_ > 0) {
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;
257
258 }
259 }
260 }
261
263 // If preconditioner type is not specified:
264 // use muelu if it is enabled, and jacobi otherwise.
265 // If eigenvalue problem type is not specified:
266 // use combinatorial for regular and
267 // normalized for irregular with polynomial preconditioner,
268 // generalized for irregular with other preconditioners.
269 // If convergence tolerance is not specified:
270 // use 1e-3 for regular with jacobi and polynomial, and 1e-2 otherwise.
271 // If how to decide the initial vectors is not specified:
272 // use random for regular and constant for irregular
274 {
275
276 // Set the default preconditioner to muelu if it is enabled, jacobi otherwise.
277 precType_ = "jacobi";
278#ifdef HAVE_ZOLTAN2SPHYNX_MUELU
279 precType_ = "muelu";
280#endif
281
282 // Override the preconditioner type with the user's preference
283 const Teuchos::ParameterEntry *pe = sphynxParams_->getEntryPtr("sphynx_preconditioner_type");
284 if (pe) {
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");
289 }
290
291 // Set the default problem type
292 problemType_ = COMBINATORIAL;
293 if(irregular_) {
294 if(precType_ == "polynomial")
295 problemType_ = NORMALIZED;
296 else
297 problemType_ = GENERALIZED;
298 }
299
300 // Override the problem type with the user's preference
301 pe = sphynxParams_->getEntryPtr("sphynx_problem_type");
302 if (pe) {
303 std::string probType = "";
304 probType = pe->getValue<std::string>(&probType);
305
306 if(probType == "combinatorial")
307 problemType_ = COMBINATORIAL;
308 else if(probType == "generalized")
309 problemType_ = GENERALIZED;
310 else if(probType == "normalized")
311 problemType_ = NORMALIZED;
312 else
313 throw std::runtime_error("\nSphynx Error: " + probType + " is not recognized as a problem type.\n"
314 + " Possible values: combinatorial, generalized, and normalized.\n");
315 }
316
317
318 // Set the default for tolerance
319 tolerance_ = 1.0e-2;
320 if(!irregular_ && (precType_ == "jacobi" || precType_ == "polynomial"))
321 tolerance_ = 1.0e-3;
322
323
324 // Override the tolerance with the user's preference
325 pe = sphynxParams_->getEntryPtr("sphynx_tolerance");
326 if (pe)
327 tolerance_ = pe->getValue<scalar_t>(&tolerance_);
328
329
330 // Set the default for initial vectors
331 randomInit_ = true;
332 if(irregular_)
333 randomInit_ = false;
334
335 // Override the initialization method with the user's preference
336 pe = sphynxParams_->getEntryPtr("sphynx_initial_guess");
337 if (pe) {
338 std::string initialGuess = "";
339 initialGuess = pe->getValue<std::string>(&initialGuess);
340
341 if(initialGuess == "random")
342 randomInit_ = true;
343 else if(initialGuess == "constants")
344 randomInit_ = false;
345 else
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");
348 }
349
350 }
351
353 // Compute the Laplacian matrix depending on the eigenvalue problem type.
354 // There are 3 options for the type: combinatorial, generalized, and normalized.
355 // Combinatorial and generalized share the same Laplacian but generalized
356 // also needs a degree matrix.
358 {
359
360 if(problemType_ == NORMALIZED)
361 laplacian_ = computeNormalizedLaplacian();
362 else
363 laplacian_ = computeCombinatorialLaplacian();
364 }
365
367 // Compute a diagonal matrix with the vertex degrees in the input graph
369 {
370
371 // Get the row pointers in the host
372 auto rowOffsets = graph_->getLocalGraphHost().row_map;
373
374 // Create the degree matrix with max row size set to 1
375 Teuchos::RCP<matrix_t> degMat(new matrix_t (graph_->getRowMap(),
376 graph_->getRowMap(),
377 1));
378
379 scalar_t *val = new scalar_t[1];
380 lno_t *ind = new lno_t[1];
381 lno_t numRows = static_cast<lno_t>(graph_->getLocalNumRows());
382
383 // Insert the diagonal entries as the degrees
384 for (lno_t i = 0; i < numRows; ++i) {
385 val[0] = rowOffsets(i+1) - rowOffsets(i) - 1;
386 ind[0] = i;
387 degMat->insertLocalValues(i, 1, val, ind);
388 }
389 delete [] val;
390 delete [] ind;
391
392 degMat->fillComplete(graph_->getDomainMap(), graph_->getRangeMap());
393
394 degMatrix_ = degMat;
395 }
396
398 // Compute the combinatorial Laplacian: L = D - A.
399 // l_ij = degree of vertex i if i = j
400 // l_ij = -1 if i != j and a_ij != 0
401 // l_ij = 0 if i != j and a_ij = 0
402 Teuchos::RCP<matrix_t> computeCombinatorialLaplacian()
403 {
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>;
408
409 const size_t numEnt = graph_->getLocalNumEntries();
410 const size_t numRows = graph_->getLocalNumRows();
411
412 // Create new values for the Laplacian, initialize to -1
413 values_view_t newVal (Kokkos::view_alloc("CombLapl::val", Kokkos::WithoutInitializing), numEnt);
414 Kokkos::deep_copy(newVal, -1);
415
416 // Get the diagonal offsets
417 offset_view_t diagOffsets(Kokkos::view_alloc("Diag Offsets", Kokkos::WithoutInitializing), numRows);
418 graph_->getLocalDiagOffsets(diagOffsets);
419
420 // Get the row pointers in the host
421 auto rowOffsets = graph_->getLocalGraphDevice().row_map;
422
423 // Compute the diagonal entries as the vertex degrees
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;
427 }
428 );
429 Kokkos::fence ();
430
431 // Create the Laplacian maatrix using the input graph and with the new values
432 Teuchos::RCP<matrix_t> laplacian (new matrix_t(graph_, newVal));
433 laplacian->fillComplete (graph_->getDomainMap(), graph_->getRangeMap());
434
435 return laplacian;
436
437 }
438
439
441 // Compute the normalized Laplacian: L_N = D^{-1/2} L D^{-1/2}, where L = D - A.
442 // l_ij = 1
443 // l_ij = -1/(sqrt(deg(v_i))sqrt(deg(v_j)) if i != j and a_ij != 0
444 // l_ij = 0 if i != j and a_ij = 0
445 Teuchos::RCP<matrix_t> computeNormalizedLaplacian()
446 {
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>;
454
455 const size_t numEnt = graph_->getLocalNumEntries();
456 const size_t numRows = graph_->getLocalNumRows();
457
458 // Create new values for the Laplacian, initialize to -1
459 values_view_t newVal (Kokkos::view_alloc("NormLapl::val", Kokkos::WithoutInitializing), numEnt);
460 Kokkos::deep_copy(newVal, -1);
461
462 // D^{-1/2}
463 dual_view_t dv ("MV::DualView", numRows, 1);
464 auto deginvsqrt = dv.d_view;
465
466 // Get the diagonal offsets
467 offset_view_t diagOffsets(Kokkos::view_alloc("Diag Offsets", Kokkos::WithoutInitializing), numRows);
468 graph_->getLocalDiagOffsets(diagOffsets);
469
470 // Get the row pointers
471 auto rowOffsets = graph_->getLocalGraphDevice().row_map;
472
473 // Compute the diagonal entries as the vertex degrees
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);
478 }
479 );
480 Kokkos::fence ();
481
482 // Create the Laplacian graph using the same graph structure with the new values
483 Teuchos::RCP<matrix_t> laplacian (new matrix_t(graph_, newVal));
484 laplacian->fillComplete (graph_->getDomainMap(), graph_->getRangeMap());
485
486 // Create the vector for D^{-1/2}
487 vector_t degInvSqrt(graph_->getRowMap(), dv);
488
489 // Normalize the laplacian matrix as D^{-1/2} L D^{-1/2}
490 laplacian->leftScale(degInvSqrt);
491 laplacian->rightScale(degInvSqrt);
492
493 return laplacian;
494 }
495
499
500 private:
501
502 // User-provided members
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_;
508
509 // Internal members
510 int numGlobalParts_;
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_;
515
516 bool irregular_; // decided internally
517 std::string precType_; // obtained from user params or decided internally
518 problemType problemType_; // obtained from user params or decided internally
519 double tolerance_; // obtained from user params or decided internally
520 bool randomInit_; // obtained from user params or decided internally
521 int verbosity_; // obtained from user params
522 bool skipPrep_; // obtained from user params
523 };
524
525
526
530
531
533 // Compute a partition using the Laplacian matrix (and possibly the degree
534 // matrix as well). First call LOBPCG to compute logK+1 eigenvectors, then
535 // transform the eigenvectors to coordinates, and finally call MJ to compute
536 // a partition on the coordinates.
537 template <typename Adapter>
538 void Sphynx<Adapter>::partition(const RCP<PartitioningSolution<Adapter>> &solution)
539 {
540 // Return a trivial solution if only one part is requested
541 if(numGlobalParts_ == 1) {
542
543 size_t numRows =adapter_->getUserGraph()->getLocalNumRows();
544 Teuchos::ArrayRCP<part_t> parts(numRows,0);
545 solution->setParts(parts);
546
547 return;
548
549 }
550
551 // The number of eigenvectors to be computed
552 int numEigenVectors = (int) log2(numGlobalParts_)+1;
553
554 // Compute the eigenvectors using LOBPCG
555 int computedNumEv = Sphynx::LOBPCGwrapper(numEigenVectors);
556
557 if(computedNumEv <= 1) {
558 throw
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");
562
563 }
564
565 // Transform the eigenvectors into coordinates
566 Teuchos::RCP<mvector_t> coordinates;
567 Sphynx::eigenvecsToCoords(eigenVectors_, computedNumEv, coordinates);
568
569 // Get the weights from the adapter
570 std::vector<const weight_t *> weights;
571 std::vector<int> wstrides;
573
574
575 // Compute the partition using MJ on coordinates
576 Sphynx::MJwrapper(coordinates, weights, wstrides, solution);
577
578 }
579
580
582 // Call LOBPCG on the Laplacian matrix.
583 template <typename Adapter>
584 int Sphynx<Adapter>::LOBPCGwrapper(const int numEigenVectors)
585 {
586
587 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer("Sphynx::LOBPCG"));
588
589 // Set defaults for the parameters
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;
596 bool lock = false;
597 bool useFullOrtho = true;
598
599 // Information to output in a verbose run
600 int numfailed = 0;
601 int iter = 0;
602 double solvetime = 0;
603
604 // Get the user values for the parameters
605 const Teuchos::ParameterEntry *pe;
606
607 pe = sphynxParams_->getEntryPtr("sphynx_maxiterations");
608 if (pe)
609 maxIterations = pe->getValue<int>(&maxIterations);
610
611 pe = sphynxParams_->getEntryPtr("sphynx_use_full_ortho");
612 if (pe)
613 useFullOrtho = pe->getValue<bool>(&useFullOrtho);
614
615
616 // Set Anasazi verbosity level
617 int anasaziVerbosity = Anasazi::Errors + Anasazi::Warnings;
618 if (verbosity_ >= 1) // low
619 anasaziVerbosity += Anasazi::FinalSummary + Anasazi::TimingDetails;
620 if (verbosity_ >= 2) // medium
621 anasaziVerbosity += Anasazi::IterationDetails;
622 if (verbosity_ >= 3) // high
623 anasaziVerbosity += Anasazi::StatusTestDetails
624 + Anasazi::OrthoDetails
625 + Anasazi::Debug;
626
627
628 // Create the parameter list to pass into solver
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);
640
641 // Create and set initial vectors
642 RCP<mvector_t> ivec( new mvector_t(laplacian_->getRangeMap(), numEigenVectors));
643
644 if (randomInit_) {
645
646 // 0-th vector constant 1, rest random
647 Anasazi::MultiVecTraits<scalar_t, mvector_t>::MvRandom(*ivec);
648 ivec->getVectorNonConst(0)->putScalar(1.);
649
650 }
651 else { // This implies we will use constant initial vectors.
652
653 // 0-th vector constant 1, other vectors constant per block
654 // Create numEigenVectors blocks, but only use numEigenVectors-1 of them.
655 // This assures orthogonality.
656 ivec->getVectorNonConst(0)->putScalar(1.);
657 for (int j = 1; j < numEigenVectors; j++)
658 ivec->getVectorNonConst(j)->putScalar(0.);
659
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\".");
663
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.);
669 }
670 }
671 }
672
673 // Create the eigenproblem to be solved
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);
678
679
680 // Set preconditioner
682
683 if(problemType_ == Sphynx::GENERALIZED)
684 problem->setM(degMatrix_);
685
686 // Inform the eigenproblem that you are finished passing it information
687 bool boolret = problem->setProblem();
688 if (boolret != true) {
689 throw std::runtime_error("\nAnasazi::BasicEigenproblem::setProblem() returned with error.\n");
690 }
691
692 // Set LOBPCG
693 using solver_t = Anasazi::LOBPCGSolMgr<scalar_t, mvector_t, op_t>;
694 solver_t solver(problem, anasaziParams);
695
696 if (verbosity_ > 0 && comm_->getRank() == 0)
697 anasaziParams.print(std::cout);
698
699 // Solve the problem
700 if (verbosity_ > 0 && comm_->getRank() == 0)
701 std::cout << "Beginning the LOBPCG solve..." << std::endl;
702 Anasazi::ReturnType returnCode = solver.solve();
703
704 // Check convergence, niters, and solvetime
705 if (returnCode != Anasazi::Converged) {
706 ++numfailed;
707 }
708 iter = solver.getNumIters();
709 solvetime = (solver.getTimers()[0])->totalElapsedTime();
710
711
712 // Retrieve the solution
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;
717
718 // Summarize iteration counts and solve time
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;
726 }
727
728 return numev;
729
730 }
731
732
733
735 template <typename Adapter>
736 template <typename problem_t>
737 void Sphynx<Adapter>::setPreconditioner(Teuchos::RCP<problem_t> &problem)
738 {
739
740 // Set the preconditioner
741 if(precType_ == "muelu") {
743 }
744 else if(precType_ == "polynomial") {
746 }
747 else if(precType_ == "jacobi") {
749 }
750 // else: do we want a case where no preconditioning is applied?
751 }
752
754 template <typename Adapter>
755 template <typename problem_t>
756 void Sphynx<Adapter>::setMueLuPreconditioner(Teuchos::RCP<problem_t> &problem)
757 {
758
759#ifdef HAVE_ZOLTAN2SPHYNX_MUELU
760 Teuchos::ParameterList paramList;
761 if(verbosity_ == 0)
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");
771
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);
780
781 if(irregular_) {
782
783 paramList.set("multigrid algorithm", "unsmoothed");
784
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);
790
791 paramList.set("max levels", 5);
792 paramList.set("aggregation: drop tol", 0.40);
793
794 }
795
796 using prec_t = MueLu::TpetraOperator<scalar_t, lno_t, gno_t, node_t>;
797 Teuchos::RCP<prec_t> prec = MueLu::CreateTpetraPreconditioner<
798 scalar_t, lno_t, gno_t, node_t>(laplacian_, paramList);
799
800 problem->setPrec(prec);
801
802#else
803 throw std::runtime_error("\nSphynx Error: MueLu requested but not compiled into Trilinos.\n");
804#endif
805
806 }
807
809 template <typename Adapter>
810 template <typename problem_t>
811 void Sphynx<Adapter>::setPolynomialPreconditioner(Teuchos::RCP<problem_t> &problem)
812 {
813 int verbosity2 = Belos::Errors;
814 if(verbosity_ == 1)
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;
823
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" );
833
834 // Construct a linear problem for the polynomial solver manager
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_);
838
839 using btop_t = Belos::TpetraOperator<scalar_t>;
840 Teuchos::RCP<btop_t> polySolver(new btop_t(innerPolyProblem,
841 Teuchos::rcpFromRef(paramList),
842 "GmresPoly", true));
843 problem->setPrec(polySolver);
844 }
845
847 template <typename Adapter>
848 template <typename problem_t>
849 void Sphynx<Adapter>::setJacobiPreconditioner(Teuchos::RCP<problem_t> &problem)
850 {
851
852 Teuchos::RCP<Ifpack2::Preconditioner<scalar_t, lno_t, gno_t, node_t>> prec;
853 std::string precType = "RELAXATION";
854
855 prec = Ifpack2::Factory::create<matrix_t> (precType, laplacian_);
856
857 Teuchos::ParameterList precParams;
858 precParams.set("relaxation: type", "Jacobi");
859 precParams.set("relaxation: fix tiny diagonal entries", true);
860 precParams.set("relaxation: min diagonal value", 1.0e-8);
861
862 prec->setParameters(precParams);
863 prec->initialize();
864 prec->compute();
865
866 problem->setPrec(prec);
867
868 }
869
871 // Transform the computed eigenvectors into coordinates
872 template <typename Adapter>
873 void Sphynx<Adapter>::eigenvecsToCoords(Teuchos::RCP<mvector_t> &eigenVectors,
874 int computedNumEv,
875 Teuchos::RCP<mvector_t> &coordinates)
876 {
877 // Extract the meaningful eigenvectors by getting rid of the first one
878 Teuchos::Array<size_t> columns (computedNumEv-1);
879 for (int j = 0; j < computedNumEv-1; ++j) {
880 columns[j] = j+1;
881 }
882 coordinates = eigenVectors->subCopy (columns());
883 coordinates->setCopyOrView (Teuchos::View);
884
885 }
886
887
889 // If user provided some weights, use them by getting them from the adapter.
890 // If user didn't provide weights but told us to use degree as weight, do so.
891 // If user neither provided weights nor told us what to do, use degree as weight.
892 template <typename Adapter>
893 void Sphynx<Adapter>::computeWeights(std::vector<const weight_t *> vecweights,
894 std::vector<int> strides)
895 {
896
897 int numWeights = adapter_->getNumWeightsPerVertex();
898 int numConstraints = numWeights > 0 ? numWeights : 1;
899
900 size_t myNumVertices = adapter_->getLocalNumVertices();
901 weight_t ** weights = new weight_t*[numConstraints];
902 for(int j = 0; j < numConstraints; j++)
903 weights[j] = new weight_t[myNumVertices];
904
905 // Will be needed if we use degree as weight
906 const offset_t *offset;
907 const gno_t *columnId;
908
909 // If user hasn't set any weights, use vertex degrees as weights
910 if(numWeights == 0) {
911
912 // Compute the weight of vertex i as the number of nonzeros in row i.
913 adapter_->getEdgesView(offset, columnId);
914 for (size_t i = 0; i < myNumVertices; i++)
915 weights[0][i] = offset[i+1] - offset[i] - 1;
916
917 vecweights.push_back(weights[0]);
918 strides.push_back(1);
919 }
920 else {
921
922 // Use the weights if there are any already set in the adapter
923 for(int j = 0; j < numConstraints; j++) {
924
925 if(adapter_->useDegreeAsVertexWeight(j)) {
926 // Compute the weight of vertex i as the number of nonzeros in row i.
927 adapter_->getEdgesView(offset, columnId);
928 for (size_t i = 0; i < myNumVertices; i++)
929 weights[j][i] = offset[i+1] - offset[i];
930 }
931 else{
932 int stride;
933 const weight_t *wgt = NULL;
934 adapter_->getVertexWeightsView(wgt, stride, j);
935
936 for (size_t i = 0; i < myNumVertices; i++)
937 weights[j][i] = wgt[i];
938 }
939
940 vecweights.push_back(weights[j]);
941 strides.push_back(1);
942
943 }
944 }
945
946 }
947
948
950 // Compute a partition by calling MJ on given coordinates with given weights
951 template <typename Adapter>
952 void Sphynx<Adapter>::MJwrapper(const Teuchos::RCP<const mvector_t> &coordinates,
953 std::vector<const weight_t *> weights,
954 std::vector<int> strides,
955 const Teuchos::RCP<PartitioningSolution<Adapter>> &solution)
956 {
957
958 Teuchos::TimeMonitor t(*Teuchos::TimeMonitor::getNewTimer("Sphynx::MJ"));
959
960 using mvector_adapter_t = Zoltan2::XpetraMultiVectorAdapter<mvector_t>;
961 using base_adapter_t = typename mvector_adapter_t::base_adapter_t;
965
966
967 size_t myNumVertices = coordinates->getLocalLength();
968
969 // Create the base adapter for the multivector adapter
970 Teuchos::RCP<mvector_adapter_t> adapcoordinates(new mvector_adapter_t(coordinates,
971 weights,
972 strides));
973 Teuchos::RCP<const base_adapter_t> baseAdapter =
974 Teuchos::rcp(dynamic_cast<const base_adapter_t *>(adapcoordinates.get()), false);
975
976 // Create the coordinate model using the base multivector adapter
978 // Teuchos::RCP<const cmodel_t> coordModel (new cmodel_t(baseAdapter, env_, comm_, flags));
979
980 // Create the MJ object
981 Teuchos::RCP<const Comm<int>> comm2 = comm_;
982 Teuchos::RCP<mj_t> mj(new mj_t(env_, comm2, baseAdapter));
983
984 // Partition with MJ
985 Teuchos::RCP<solution_t> vectorsolution( new solution_t(env_, comm2, 1, mj));
986 mj->partition(vectorsolution);
987
988 // Transform the solution
989 Teuchos::ArrayRCP<part_t> parts(myNumVertices);
990 for(size_t i = 0; i < myNumVertices; i++) {
991 parts[i] = vectorsolution->getPartListView()[i];
992 }
993 solution->setParts(parts);
994
995 }
996
997} // namespace Zoltan2
998
999#endif
Contains the Multi-jagged algorthm.
Defines the CoordinateModel classes.
Defines XpetraCrsGraphAdapter class.
Defines the XpetraMultiVectorAdapter.
Algorithm defines the base class for all algorithms.
This class provides geometric coordinates with optional weights to the Zoltan2 algorithm.
A PartitioningSolution is a solution to a partitioning problem.
Sphynx(const RCP< const Environment > &env, const RCP< Teuchos::ParameterList > &params, const RCP< Teuchos::ParameterList > &sphynxParams, const RCP< const Comm< int > > &comm, const RCP< const XpetraCrsGraphAdapter< graph_t > > &adapter)
void setMueLuPreconditioner(Teuchos::RCP< problem_t > &problem)
Teuchos::RCP< matrix_t > computeNormalizedLaplacian()
typename Adapter::offset_t offset_t
int LOBPCGwrapper(const int numEigenVectors)
void MJwrapper(const Teuchos::RCP< const mvector_t > &coordinates, std::vector< const weight_t * > weights, std::vector< int > strides, const Teuchos::RCP< PartitioningSolution< Adapter > > &solution)
Tpetra::CrsGraph< lno_t, gno_t, node_t > graph_t
void partition(const RCP< PartitioningSolution< Adapter > > &solution)
Partitioning method.
void eigenvecsToCoords(Teuchos::RCP< mvector_t > &eigenVectors, int computedNumEv, Teuchos::RCP< mvector_t > &coordinates)
Tpetra::MultiVector< scalar_t, lno_t, gno_t, node_t > mvector_t
typename Adapter::node_t node_t
Tpetra::Operator< scalar_t, lno_t, gno_t, node_t > op_t
typename Adapter::lno_t lno_t
typename Adapter::scalar_t weight_t
Teuchos::RCP< matrix_t > computeCombinatorialLaplacian()
typename Adapter::gno_t gno_t
void setPreconditioner(Teuchos::RCP< problem_t > &problem)
void setPolynomialPreconditioner(Teuchos::RCP< problem_t > &problem)
void setJacobiPreconditioner(Teuchos::RCP< problem_t > &problem)
Tpetra::CrsMatrix< scalar_t, lno_t, gno_t, node_t > matrix_t
typename Adapter::part_t part_t
void computeWeights(std::vector< const weight_t * > vecweights, std::vector< int > strides)
Multi Jagged coordinate partitioning algorithm.
map_t::local_ordinal_type lno_t
map_t::global_ordinal_type gno_t
Created by mbenlioglu on Aug 31, 2020.
std::bitset< NUM_MODEL_FLAGS > modelFlag_t
static RCP< tMVector_t > coordinates
static ArrayRCP< ArrayRCP< zscalar_t > > weights