All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MultiLevelSetup.h
Go to the documentation of this file.
1 /*
2  * MultiLevelSetup.h
3  *
4  * Created on: 10.01.2015
5  * Author: Michael Wegner (michael.wegner@student.kit.edu)
6  */
7 
8 #ifndef MULTILEVELSETUP_H_
9 #define MULTILEVELSETUP_H_
10 
11 #include "LevelHierarchy.h"
12 #include "../Smoother.h"
13 #include "../../algebraic/CSRMatrix.h"
14 
15 #include <limits>
16 
17 namespace NetworKit {
18 
23 template<class Matrix>
25 
26 #define UNDECIDED std::numeric_limits<index>::max()
27 
28 private:
29  const Smoother<Matrix> &smoother;
30 
39  bool coarseningElimination(Matrix& matrix, LevelHierarchy<Matrix>& hierarchy) const;
40 
50  count lowDegreeSweep(const Matrix& matrix, std::vector<bool>& fNode, index stage) const;
51 
61  void eliminationOperators(const Matrix& matrix, const std::vector<index>& fSet, const std::vector<index>& coarseIndex, Matrix& P, Vector& q) const;
62 
72  void coarseningAggregation(Matrix& matrix, LevelHierarchy<Matrix>& hierarchy, Vector& tv, count numTVVectors) const;
73 
82  std::vector<Vector> generateTVs(const Matrix& matrix, Vector& tv, const count numVectors) const;
83 
89  void addHighDegreeSeedNodes(const Matrix& matrix, std::vector<index>& status) const;
90 
97  void aggregateLooseNodes(const Matrix& strongAdjMatrix, std::vector<index>& status, count& nc) const;
98 
104  void computeStrongAdjacencyMatrix(const Matrix& matrix, Matrix& strongAdjMatrix) const;
105 
112  void computeAffinityMatrix(const Matrix& matrix, const std::vector<Vector>& tVs, Matrix& affinityMatrix) const;
113 
123  void aggregationStage(const Matrix& matrix, count& nc, const Matrix& strongAdjMatrix, const Matrix& affinityMatrix, std::vector<Vector>& tVs, std::vector<index>& status) const;
124 
131  void computeStrongNeighbors(const Matrix& affinityMatrix, const std::vector<index>& status, std::vector<std::vector<index>>& bins) const;
132 
144  bool findBestSeedEnergyCorrected(const Matrix& strongAdjMatrix, const Matrix& affinityMatrix, const std::vector<double>& diag, const std::vector<Vector>& tVs, const std::vector<index>& status, const index u, index& s) const;
145 
146 
152  inline bool canCoarsen(const Matrix& A) const {
153  return A.numberOfRows() > MAX_DIRECT_SOLVE_SIZE;
154  }
155 
163  bool isRelaxationFast(const Matrix& A, index lvlIndex, Vector& tv) const;
164 
173  void galerkinOperator(const Matrix& P, const Matrix& A, const std::vector<index>& PColIndex, const std::vector<std::vector<index>>& PRowIndex, Matrix& B) const;
174 
180  void setupForMatrix(Matrix& matrix, LevelHierarchy<Matrix>& hierarchy) const;
181 
182 public:
187  MultiLevelSetup(const Smoother<Matrix>& smoother) : smoother(smoother) {}
188 
194  void setup(const Graph& G, LevelHierarchy<Matrix>& hierarchy) const {
195  setup(Matrix::laplacianMatrix(G), hierarchy);
196  }
197 
203  void setup(const Matrix& matrix, LevelHierarchy<Matrix>& hierarchy) const;
204 };
205 
206 template<class Matrix>
207 void MultiLevelSetup<Matrix>::setup(const Matrix& matrix, LevelHierarchy<Matrix>& hierarchy) const {
208  CSRMatrix A = matrix;
209  setupForMatrix(A, hierarchy);
210 }
211 
212 template<class Matrix>
213 void MultiLevelSetup<Matrix>::setupForMatrix(Matrix& A, LevelHierarchy<Matrix>& hierarchy) const {
214  hierarchy.addFinestLevel(A);
215 
216 #ifndef NDEBUG
217  DEBUG("FINEST\t", A.numberOfRows(), "\t", A.nnz());
218 #endif
219 
220  bool doneCoarsening = false;
221  count numTVs = TV_NUM;
222  index level = 0;
223 
224  while (!doneCoarsening) {
225  // ELIMINATION
226  if (coarseningElimination(A, hierarchy)) {
227  if (!canCoarsen(A)) doneCoarsening = true;
228  level++;
229 #ifndef NDEBUG
230  DEBUG(level, " ELIM\t\t", A.numberOfRows(), "\t", A.nnz() / 2);
231 #endif
232  }
233 
234  // AGGREGATION
235  Vector tv;
236  if (doneCoarsening || isRelaxationFast(A, level, tv)) {
237  doneCoarsening = true;
238  } else {
239  coarseningAggregation(A, hierarchy, tv, numTVs);
240  level++;
241 #ifndef NDEBUG
242  DEBUG(level, " AGG\t\t", A.numberOfRows(), "\t", A.nnz() / 2);
243 #endif
244  if (numTVs < TV_MAX) {
245  numTVs += TV_INC;
246  }
247  }
248 
249  if (!canCoarsen(A)) doneCoarsening = true;
250  }
251 
252  hierarchy.setLastAsCoarsest();
253 }
254 
255 
256 template<class Matrix>
257 bool MultiLevelSetup<Matrix>::coarseningElimination(Matrix& matrix, LevelHierarchy<Matrix>& hierarchy) const {
258  std::vector<EliminationStage<Matrix>> coarseningStages;
259  count stageNum = 0;
260  while (stageNum < SETUP_ELIMINATION_MAX_STAGES) {
261  if (matrix.numberOfRows() <= MAX_DIRECT_SOLVE_SIZE) break; // we do not need to coarsen the matrix any further
262  std::vector<bool> fNode;
263  count nf = lowDegreeSweep(matrix, fNode, stageNum);
264  count nc = matrix.numberOfRows() - nf;
265 
266  if (nc == 0) { // do not eliminate all nodes -> leave one entry in c
267  nc = 1;
268  nf--;
269  }
270 
271  // add f nodes to fSet and c nodes to cSet
272  std::vector<index> fSet(nf);
273  std::vector<index> cSet(nc);
274 
275  std::vector<index> coarseIndex(matrix.numberOfRows());
276  count numFNodes = 0;
277  for (index i = 0, fIndex = 0, cIndex = 0; i < matrix.numberOfRows(); ++i) {
278  if (fNode[i] && fIndex < nf) {
279  coarseIndex[i] = fIndex;
280  fSet[fIndex++] = i;
281  numFNodes++;
282  } else {
283  coarseIndex[i] = cIndex;
284  cSet[cIndex++] = i;
285  }
286  }
287 
288  if (nf <= SETUP_ELIMINATION_MIN_ELIM_FRACTION * matrix.numberOfRows()) {
289  break;
290  }
291 
292  Matrix P;
293  Vector q;
294  eliminationOperators(matrix, fSet, coarseIndex, P, q);
295  coarseningStages.push_back(EliminationStage<Matrix>(P, q, fSet, cSet));
296 
297  Matrix Acc = matrix.extract(cSet, cSet); // Schur complement
298  Matrix Acf = matrix.extract(cSet, fSet); // Schur complement
299 
300  matrix = Acc + Acf * P;
301  stageNum++;
302  }
303 
304  if (stageNum != 0) { // we have coarsened the matrix
305  hierarchy.addEliminationLevel(matrix, coarseningStages);
306  return true;
307  }
308 
309  return false;
310 }
311 
312 template<class Matrix>
313 count MultiLevelSetup<Matrix>::lowDegreeSweep(const Matrix& matrix, std::vector<bool>& fNode, index stage) const {
314  fNode.resize(matrix.numberOfRows(), true); // first mark all nodes as f nodes
315  count numFNodes = 0;
316  int degreeOffset = stage != 0;
317 
318  for (index i = 0; i < matrix.numberOfRows(); ++i) {
319  if ((int) matrix.nnzInRow(i) - degreeOffset <= (int)SETUP_ELIMINATION_MAX_DEGREE && fNode[i]) { // node i has degree <= 4 and can be eliminated
320  numFNodes++;
321  matrix.forNonZeroElementsInRow(i, [&](index j, edgeweight /*w*/){ // to maintain independence, mark all neighbors as not eliminated
322  if (j != i) { // all neighbors of this f node are c nodes
323  fNode[j] = false;
324  }
325  });
326  } else { // node has high degree, thus it is a c node
327  fNode[i] = false;
328  }
329  }
330 
331  return numFNodes;
332 }
333 
334 template<class Matrix>
335 void MultiLevelSetup<Matrix>::eliminationOperators(const Matrix& matrix, const std::vector<index>& fSet, const std::vector<index>& coarseIndex, Matrix& P, Vector& q) const {
336  std::vector<Triplet> triples;
337  q = Vector(fSet.size());
338  for (index k = 0; k < fSet.size(); ++k) { // Afc
339  matrix.forNonZeroElementsInRow(fSet[k], [&](index j, edgeweight w){
340  if (fSet[k] == j) {
341  q[k] = 1.0 / w;
342  } else {
343  triples.push_back({k, coarseIndex[j], w});
344  }
345  });
346  }
347 
348  for (index i = 0; i < triples.size(); ++i) { // * -Aff^-1
349  triples[i].value *= -q[triples[i].row];
350  }
351 
352  P = Matrix(fSet.size(), coarseIndex.size() - fSet.size(), triples);
353 }
354 
355 
356 
357 template<class Matrix>
358 void MultiLevelSetup<Matrix>::coarseningAggregation(Matrix& matrix, LevelHierarchy<Matrix>& hierarchy, Vector& tv, count numTVVectors) const {
359  Vector B(SETUP_MAX_AGGREGATION_STAGES, std::numeric_limits<double>::max());
360  std::vector<std::vector<index>> S(SETUP_MAX_AGGREGATION_STAGES, std::vector<index>(matrix.numberOfRows(), std::numeric_limits<index>::max()));
361  std::vector<index> status(matrix.numberOfRows(), UNDECIDED);
362  std::vector<count> nc(SETUP_MAX_AGGREGATION_STAGES, matrix.numberOfRows());
363 
364  double alpha = 1.0;
365  double maxCoarseningRatio = SETUP_COARSENING_WORK_GUARD / SETUP_CYCLE_INDEX;
366  count stage = 0;
367  count nC = matrix.numberOfRows();
368 
369 
370  // generate TVs
371  std::vector<Vector> tVs = generateTVs(matrix, tv, numTVVectors);
372 
373  // compute strong adjacency matrix
374  Matrix Wstrong;
375  computeStrongAdjacencyMatrix(matrix, Wstrong);
376 
377  // compute affinityMatrix
378  Matrix affinityMatrix;
379  computeAffinityMatrix(Wstrong, tVs, affinityMatrix);
380 
381  // mark all locally high-degree nodes as seeds
382  addHighDegreeSeedNodes(matrix, status);
383 
384  // aggregate all loose nodes
385  aggregateLooseNodes(Wstrong, status, nC);
386 
387  nc[0] = nC;
388  while (stage < SETUP_MIN_AGGREGATION_STAGES || (alpha >= maxCoarseningRatio && stage < SETUP_MAX_AGGREGATION_STAGES)) {
389  nC = stage > 0? nc[stage - 1] : nc[0];
390 
391  // aggregation stage
392  aggregationStage(matrix, nC, Wstrong, affinityMatrix, tVs, status);
393 
394  alpha = (double) nC / (double) matrix.numberOfRows();
395  alpha <= maxCoarseningRatio? B[stage] = 1.0-alpha : B[stage] = 1.0+alpha;
396 
397  S[stage] = status;
398  nc[stage] = nC;
399  stage++;
400  }
401 
402  double min = B[0];
403  index bestAggregate = 0;
404  for (index i = 1; i < stage; ++i) {
405  if (B[i] < min) {
406  bestAggregate = i;
407  min = B[i];
408  }
409  }
410 
411  for (index i = 0; i < matrix.numberOfRows(); ++i) {
412  if (S[bestAggregate][i] == UNDECIDED) { // undediced nodes become their own seeds
413  S[bestAggregate][i] = i;
414  }
415  }
416 
417  std::vector<index> indexFine(matrix.numberOfRows(), 0);
418  index newIndex = 0;
419  for (index i = 0; i < matrix.numberOfRows(); ++i) {
420  if (S[bestAggregate][i] == i) {
421  indexFine[i] = newIndex++;
422  }
423  }
424 
425  for (index i = 0; i < matrix.numberOfRows(); ++i) {
426  status[i] = indexFine[S[bestAggregate][i]];
427  }
428 
429  assert(newIndex == nc[bestAggregate]);
430 
431  // create interpolation matrix
432  std::vector<Triplet> pTriples(matrix.numberOfRows());
433  std::vector<Triplet> rTriples(matrix.numberOfRows());
434  std::vector<index> PColIndex(matrix.numberOfRows());
435  std::vector<std::vector<index>> PRowIndex(nc[bestAggregate]);
436 
437  for (index i = 0; i < matrix.numberOfRows(); ++i) {
438  pTriples[i] = {i, status[i], 1};
439  rTriples[i] = {status[i], i, 1};
440  PColIndex[i] = status[i];
441  PRowIndex[status[i]].push_back(i);
442  }
443 
444  Matrix P(matrix.numberOfRows(), nc[bestAggregate], pTriples);
445  Matrix R(nc[bestAggregate], matrix.numberOfRows(), rTriples);
446 
447  // create coarsened laplacian
448  galerkinOperator(P, matrix, PColIndex, PRowIndex, matrix);
449 
450  hierarchy.addAggregationLevel(matrix, P, R);
451 }
452 
453 template<class Matrix>
454 std::vector<Vector> MultiLevelSetup<Matrix>::generateTVs(const Matrix& matrix, Vector& tv, count numVectors) const {
455  std::vector<Vector> testVectors(numVectors, Vector(matrix.numberOfColumns()));
456 
457  testVectors[0] = tv;
458 
459  if (numVectors > 1) {
460  Vector b(matrix.numberOfColumns(), 0.0);
461 #pragma omp parallel for
462  for (count i = 1; i < numVectors; ++i) {
463  for (count j = 0; j < matrix.numberOfColumns(); ++j) {
464  testVectors[i][j] = 2 * Aux::Random::probability() - 1;
465  }
466 
467  testVectors[i] = smoother.relax(matrix, b, testVectors[i], SETUP_TV_SWEEPS);
468  }
469  }
470 
471  return testVectors;
472 }
473 
474 template<class Matrix>
475 void MultiLevelSetup<Matrix>::addHighDegreeSeedNodes(const Matrix& matrix, std::vector<index>& status) const {
476  std::vector<count> deg(matrix.numberOfRows());
477 #pragma omp parallel for
478  for (index i = 0; i < matrix.numberOfRows(); ++i) {
479  deg[i] = matrix.nnzInRow(i) - 1;
480  }
481 
482 #pragma omp parallel for
483  for (index i = 0; i < matrix.numberOfRows(); ++i) {
484  double num = 0.0;
485  double denom = 0.0;
486  matrix.forNonZeroElementsInRow(i, [&](index j, double value){
487  if (i != j) {
488  num += std::abs(value) * (double) deg[j];
489  } else {
490  denom = std::abs(value);
491  }
492  });
493 
494 
495  if ((double) deg[i] >= SETUP_AGGREGATION_DEGREE_THRESHOLD * (num / denom)) { // high degree node becomes seed
496  status[i] = i;
497  }
498  }
499 }
500 
501 template<class Matrix>
502 void MultiLevelSetup<Matrix>::aggregateLooseNodes(const Matrix& strongAdjMatrix, std::vector<index>& status, count& nc) const {
503  std::vector<index> looseNodes;
504  for (index i = 0; i < strongAdjMatrix.numberOfRows(); ++i) {
505  double max = std::numeric_limits<double>::min();
506  strongAdjMatrix.forNonZeroElementsInRow(i, [&](index /*j*/, double value) {
507  if (value > max) max = value;
508  });
509 
510  if (std::abs(max) < 1e-9 || max == std::numeric_limits<double>::min()) {
511  looseNodes.push_back(i);
512  }
513  }
514 
515  if (looseNodes.size() > 0) {
516  status[looseNodes[0]] = looseNodes[0]; // mark first as seed
517  for (index k = 1; k < looseNodes.size(); ++k) {
518  status[looseNodes[k]] = looseNodes[0]; // first loose nodes becomes seed
519  }
520 
521  nc -= looseNodes.size() - 1;
522  }
523 }
524 
525 
526 
527 template<class Matrix>
528 void MultiLevelSetup<Matrix>::computeStrongAdjacencyMatrix(const Matrix& matrix, Matrix& strongAdjMatrix) const {
529  std::vector<double> maxNeighbor(matrix.numberOfRows(), std::numeric_limits<double>::min());
530 #pragma omp parallel for
531  for (index i = 0; i < matrix.numberOfRows(); ++i) {
532  matrix.forNonZeroElementsInRow(i, [&](index j, double value) {
533  if (i != j && -value > maxNeighbor[i]) {
534  maxNeighbor[i] = -value;
535  }
536  });
537  }
538 
539  std::vector<index> rowIdx(matrix.numberOfRows()+1, 0);
540 
541  matrix.parallelForNonZeroElementsInRowOrder([&](index i, index j, double value) {
542  if (i != j && std::abs(value) >= 0.1 * std::min(maxNeighbor[i], maxNeighbor[j])) {
543  ++rowIdx[i+1];
544  }
545  });
546 
547  for (index i = 0; i < matrix.numberOfRows(); ++i) {
548  rowIdx[i+1] += rowIdx[i];
549  }
550 
551  count nnz = rowIdx[matrix.numberOfRows()];
552  std::vector<Triplet> triplets(nnz);
553 
554 #pragma omp parallel for
555  for (index i = 0; i < matrix.numberOfRows(); ++i) {
556  index cIdx = rowIdx[i];
557  matrix.forNonZeroElementsInRow(i, [&](index j, double value) {
558  if (i != j && std::abs(value) >= 0.1 * std::min(maxNeighbor[i], maxNeighbor[j])) {
559  triplets[cIdx] = {i,j,-value};
560  ++cIdx;
561  }
562  });
563  }
564 
565  strongAdjMatrix = Matrix(matrix.numberOfRows(), matrix.numberOfColumns(), triplets);
566 }
567 
568 
569 template<class Matrix>
570 void MultiLevelSetup<Matrix>::computeAffinityMatrix(const Matrix& matrix, const std::vector<Vector>& tVs, Matrix& affinityMatrix) const {
571  assert(tVs.size() > 0);
572 
573  std::vector<index> rowIdx(matrix.numberOfRows()+1);
574  std::vector<Triplet> triplets(matrix.nnz());
575 
576 #pragma omp parallel for
577  for (index i = 0; i < matrix.numberOfRows(); ++i) {
578  rowIdx[i+1] = matrix.nnzInRow(i);
579  }
580 
581  for (index i = 0; i < matrix.numberOfRows(); ++i) {
582  rowIdx[i+1] += rowIdx[i];
583  }
584 
585  std::vector<double> normSquared(matrix.numberOfRows(), 0.0);
586 #pragma omp parallel for
587  for (index i = 0; i < matrix.numberOfRows(); ++i) {
588  for (index k = 0; k < tVs.size(); ++k) {
589  normSquared[i] += tVs[k][i] * tVs[k][i];
590  }
591  }
592 
593 #pragma omp parallel for
594  for (index i = 0; i < matrix.numberOfRows(); ++i) {
595  double nir = 1.0 / normSquared[i];
596  index cIdx = rowIdx[i];
597  matrix.forNonZeroElementsInRow(i, [&](index j, double /*val*/) {
598  double ij = 0.0;
599  for (index k = 0; k < tVs.size(); ++k) {
600  ij += tVs[k][i] * tVs[k][j];
601  }
602 
603  double value = (ij * ij) * nir / normSquared[j];
604  triplets[cIdx] = {i,j,value};
605  ++cIdx;
606  });
607  }
608 
609  affinityMatrix = Matrix(matrix.numberOfRows(), matrix.numberOfColumns(), triplets);
610 }
611 
612 template<class Matrix>
613 void MultiLevelSetup<Matrix>::aggregationStage(const Matrix& matrix, count& nc, const Matrix& strongAdjMatrix, const Matrix& affinityMatrix, std::vector<Vector>& tVs, std::vector<index>& status) const {
614  std::vector<std::vector<index>> bins(10);
615  computeStrongNeighbors(affinityMatrix, status, bins);
616 
617  std::vector<double> diag(matrix.numberOfRows(), 0.0);
618 #pragma omp parallel for
619  for (index i = 0 ; i < matrix.numberOfRows(); ++i) {
620  diag[i] = matrix(i,i);
621  }
622 
623  for (index k = bins.size(); k-- > 0;) { // iterate over undecided nodes with strong neighbors in decreasing order of strongest neighbor
624  for (index i : bins[k]) {
625  if (status[i] == UNDECIDED) { // node is still undecided
626  index s = 0;
627  if (findBestSeedEnergyCorrected(strongAdjMatrix, affinityMatrix, diag, tVs, status, i, s)) {
628  status[s] = s; // s becomes seed
629  status[i] = s; // i's seed is s
630  nc--;
631 
632  for (index j = 0; j < tVs.size(); ++j) { // update test vectors
633  tVs[j][i] = tVs[j][s];
634  }
635  }
636  }
637  }
638 
639  if (nc <= matrix.numberOfRows() * SETUP_COARSENING_WORK_GUARD / SETUP_CYCLE_INDEX) {
640  break;
641  }
642  } // iterate over bins
643 }
644 
645 template<class Matrix>
646 void MultiLevelSetup<Matrix>::computeStrongNeighbors(const Matrix& affinityMatrix, const std::vector<index>& status, std::vector<std::vector<index>>& bins) const {
647  std::vector<bool> undecided(affinityMatrix.numberOfRows(), false);
648  std::vector<double> maxNeighbor(affinityMatrix.numberOfRows(), std::numeric_limits<double>::min());
649  double overallMax = 0.0;
650  double overallMin = std::numeric_limits<double>::max();
651 
652  affinityMatrix.parallelForNonZeroElementsInRowOrder([&](index i, index j, double value) { // determine the highest affinity neighbor of each node
653  if (status[i] == UNDECIDED && (status[j] == UNDECIDED || status[j] == j)) { // i is UNDECIDED and its neighbor j is also UNDECIDED or SEED
654  if (value > maxNeighbor[i]) {
655  maxNeighbor[i] = value;
656  }
657  undecided[i] = true;
658  }
659  });
660 
661  for (index i = 0; i < affinityMatrix.numberOfRows(); ++i) {
662  if (maxNeighbor[i] > overallMax) {
663  overallMax = maxNeighbor[i];
664  }
665  if (maxNeighbor[i] < overallMin) {
666  overallMin = maxNeighbor[i];
667  }
668  }
669 
670  double h = fabs(overallMax - overallMin) < 1e-15? 1.0 : (double) bins.size() / (overallMax - overallMin);
671  for (index i = 0; i < affinityMatrix.numberOfRows(); ++i) {
672  if (undecided[i]) { // undecided nodes with strong neighbors
673  index binIndex = (index) std::floor(h * (maxNeighbor[i] - overallMin));
674  if (binIndex == bins.size()) { // last interval is closed on the right
675  binIndex--;
676  }
677 
678  assert(binIndex >= 0 && binIndex < bins.size());
679  bins[binIndex].push_back(i);
680  }
681  }
682 }
683 
684 template<class Matrix>
685 bool MultiLevelSetup<Matrix>::findBestSeedEnergyCorrected(const Matrix& strongAdjMatrix, const Matrix& affinityMatrix, const std::vector<double>& diag, const std::vector<Vector>& tVs, const std::vector<index>& status, const index u, index& s) const {
686  bool foundSeed = false;
687  std::vector<double> r(tVs.size(), 0.0);
688  std::vector<double> q(tVs.size(), 0.0);
689  std::vector<double> E(tVs.size(), 0.0);
690 
691  double d = diag[u];
692  double d2 = 0.5 * diag[u];
693  for (index k = 0; k < tVs.size(); ++k) {
694  double rr = 0.0;
695  double qq = 0.0;
696  strongAdjMatrix.forNonZeroElementsInRow(u, [&](index v, double value) {
697  rr += value * tVs[k][v];
698  qq += value * 0.5 * tVs[k][v] * tVs[k][v];
699  });
700 
701  r[k] = rr;
702  q[k] = qq;
703  double y = rr/d;
704  E[k] = (d2*y - rr)*y + qq;
705  }
706 
707  double maxNeighbor = -1.0;
708  affinityMatrix.forNonZeroElementsInRow(u, [&](index v, double value) {
709  if (status[v] == UNDECIDED || status[v] == v) {
710  double maxMu = -1.0;
711  bool smallRatio = true;
712  for (index k = 0; k < tVs.size(); ++k) {
713  double xv = tVs[k][v];
714  double Ec = (d2*xv - r[k])*xv + q[k];
715  double mu = Ec / (E[k] + 1e-15);
716 
717  if (mu > maxMu) {
718  maxMu = mu;
719  }
720  if (maxMu > 2.5) {
721  smallRatio = false;
722  break;
723  }
724  }
725 
726  if (smallRatio && value > maxNeighbor) {
727  maxNeighbor = value;
728  s = v;
729  foundSeed = true;
730  }
731  }
732  });
733 
734  return foundSeed;
735 }
736 
737 template<class Matrix>
738 bool MultiLevelSetup<Matrix>::isRelaxationFast(const Matrix& A, index lvlIndex, Vector& tv) const {
739  count nu = SETUP_RELAX_ACF_MIN_SWEEPS + 2 * (lvlIndex - 1);
740  count tvNu = SETUP_TV_SWEEPS;
741  count initial = 3;
742 
743  // create testVector in [-1,1]
744  tv = Vector(A.numberOfRows());
745  for (index i = 0; i < tv.getDimension(); ++i) {
746  tv[i] = 2.0 * Aux::Random::probability() - 1.0;
747  }
748 
749  Vector b(A.numberOfRows(), 0.0);
750  Vector x = tv;
751  x = smoother.relax(A, b, x, initial);
752  tv = smoother.relax(A, b, x, tvNu - initial);
753  Vector y = smoother.relax(A, b, tv, nu - tvNu);
754  double relaxAcf = std::pow((y - y.mean()).length() / (x - x.mean()).length(), (double) 1.0 / (double) (nu - initial));
755  return relaxAcf <= SETUP_MAX_COARSE_RELAX_ACF || !canCoarsen(A);
756 }
757 
758 
759 
760 template<class Matrix>
761 void MultiLevelSetup<Matrix>::galerkinOperator(const Matrix& P, const Matrix& A, const std::vector<index>& PColIndex, const std::vector<std::vector<index>>& PRowIndex, Matrix& B) const {
762  std::vector<Triplet> triplets;
763  SparseAccumulator spa(P.numberOfColumns());
764  for (index i = 0; i < P.numberOfColumns(); ++i) {
765  for (index k : PRowIndex[i]) {
766  double Pki = P(k,i);
767  A.forNonZeroElementsInRow(k, [&](index l, double value) {
768  index j = PColIndex[l];
769  spa.scatter(Pki * value * P(l, j), j);
770  });
771  }
772 
773  spa.gather([&](index i, index j, double value) {
774  triplets.push_back({i,j,value});
775  });
776 
777  spa.increaseRow();
778  }
779 
780  B = CSRMatrix(P.numberOfColumns(), P.numberOfColumns(), triplets);
781 }
782 
783 } /* namespace NetworKit */
784 
785 #endif /* MULTILEVELSETUP_H_ */
constexpr count SETUP_MIN_AGGREGATION_STAGES
Definition: LAMGSettings.h:61
constexpr double SETUP_COARSENING_WORK_GUARD
Definition: LAMGSettings.h:59
constexpr count TV_INC
Definition: LAMGSettings.h:18
constexpr double SETUP_ELIMINATION_MIN_ELIM_FRACTION
Definition: LAMGSettings.h:47
constexpr count SETUP_MAX_AGGREGATION_STAGES
Definition: LAMGSettings.h:63
constexpr count SETUP_RELAX_ACF_MIN_SWEEPS
Definition: LAMGSettings.h:32
#define UNDECIDED
Definition: MultiLevelSetup.h:26
uint64_t index
Typedefs.
Definition: Globals.h:20
#define DEBUG(...)
Definition: Log.h:82
constexpr count SETUP_TV_SWEEPS
Definition: LAMGSettings.h:22
constexpr double SETUP_MAX_COARSE_RELAX_ACF
Definition: LAMGSettings.h:34
constexpr count TV_MAX
Definition: LAMGSettings.h:20
constexpr count SETUP_ELIMINATION_MAX_STAGES
Definition: LAMGSettings.h:45
uint64_t count
Definition: Globals.h:21
std::vector< std::vector< count > > Matrix
Definition: DynamicNMIDistance.h:16
constexpr count TV_NUM
Definition: LAMGSettings.h:16
Implements the setup phase of LAMG (Lean Algebraic Multigrid by Livne et al.).
Definition: MultiLevelSetup.h:24
constexpr count SETUP_ELIMINATION_MAX_DEGREE
Definition: LAMGSettings.h:43
constexpr count MAX_DIRECT_SOLVE_SIZE
Definition: LAMGSettings.h:24
constexpr count SETUP_AGGREGATION_DEGREE_THRESHOLD
Definition: LAMGSettings.h:55
constexpr double SETUP_CYCLE_INDEX
Definition: LAMGSettings.h:30
void setup(const CSRMatrix &matrix, LevelHierarchy< CSRMatrix > &hierarchy) const
Definition: MultiLevelSetup.cpp:13
double probability()
Definition: Random.cpp:90
double edgeweight
Definition: Globals.h:24