Lamg.h
Go to the documentation of this file.
1 /*
2  * Lamg.h
3  *
4  * Created on: Oct 20, 2015
5  * Author: Michael Wegner (michael.wegner@student.kit.edu)
6  */
7
8 #ifndef NETWORKIT_CPP_NUMERICS_LAMG_LAMG_H_
9 #define NETWORKIT_CPP_NUMERICS_LAMG_LAMG_H_
10
11 #include <vector>
12
13 #include "../LinearSolver.h"
14 #include "MultiLevelSetup.h"
15 #include "SolverLamg.h"
16 #include "../GaussSeidelRelaxation.h"
17 #include "../../algebraic/MatrixTools.h"
18 #include "../../components/ParallelConnectedComponents.h"
19 #include "omp.h"
20
21 namespace NetworKit {
22
29 template<class Matrix>
30 class Lamg : public LinearSolver<Matrix> {
31 private:
32  bool validSetup;
34  MultiLevelSetup<Matrix> lamgSetup;
35  Matrix laplacianMatrix;
36  std::vector<LevelHierarchy<Matrix>> compHierarchies;
37  std::vector<SolverLamg<Matrix>> compSolvers;
38  std::vector<LAMGSolverStatus> compStati;
39
40  std::vector<Vector> initialVectors;
41  std::vector<Vector> rhsVectors;
42
43  count numComponents;
44  std::vector<std::vector<index>> components;
45  std::vector<index> graph2Components;
46
47  void initializeForOneComponent();
48
49 public:
55  Lamg(const double tolerance = 1e-6) : LinearSolver<Matrix>(tolerance), validSetup(false), lamgSetup(smoother), numComponents(0) {}
57  ~Lamg() = default;
58
65  void setup(const Matrix& laplacianMatrix);
66
72  void setupConnected(const Matrix& laplacianMatrix);
73
84  SolverStatus solve(const Vector& rhs, Vector& result, count maxConvergenceTime = 5 * 60 * 1000, count maxIterations = std::numeric_limits<count>::max());
85
95  void parallelSolve(const std::vector<Vector>& rhs, std::vector<Vector>& results, count maxConvergenceTime = 5 * 60 * 1000, count maxIterations = std::numeric_limits<count>::max());
96
97 };
98
99 template<class Matrix>
100 void Lamg<Matrix>::initializeForOneComponent() {
101  compHierarchies = std::vector<LevelHierarchy<Matrix>>(1);
102  lamgSetup.setup(laplacianMatrix, compHierarchies[0]);
103  compSolvers.clear();
104  compSolvers.push_back(SolverLamg<Matrix>(compHierarchies[0], smoother));
105  validSetup = true;
106 }
107
108 template<class Matrix>
109 void Lamg<Matrix>::setupConnected(const Matrix& laplacianMatrix) {
110  this->laplacianMatrix = laplacianMatrix;
111  initializeForOneComponent();
112  numComponents = 1;
113 }
114
115 template<class Matrix>
116 void Lamg<Matrix>::setup(const Matrix& laplacianMatrix) {
117  this->laplacianMatrix = laplacianMatrix;
118  Graph G = MatrixTools::matrixToGraph(laplacianMatrix);
119  ParallelConnectedComponents con(G, false);
120  con.run();
121  numComponents = con.numberOfComponents();
122  if (numComponents == 1) {
123  initializeForOneComponent();
124  } else {
125  graph2Components = std::vector<index>(G.numberOfNodes());
126
127  initialVectors = std::vector<Vector>(numComponents);
128  rhsVectors = std::vector<Vector>(numComponents);
129
130  components = std::vector<std::vector<index>>(numComponents);
131  compHierarchies = std::vector<LevelHierarchy<Matrix>>(numComponents);
132  compSolvers.clear();
133  compStati = std::vector<LAMGSolverStatus>(numComponents);
134
135  // create solver for every component
136  index compIdx = 0;
137  for (auto component : con.getPartition().getSubsets()) {
138  components[compIdx] = std::vector<index>(component.begin(), component.end());
139
140  std::vector<Triplet> triplets;
141  index idx = 0;
142  for (node u : components[compIdx]) {
143  graph2Components[u] = idx;
144  idx++;
145  }
146
147  for (node u : components[compIdx]) {
148  G.forNeighborsOf(u, [&](node v, edgeweight w) {
149  triplets.push_back({graph2Components[u], graph2Components[v], w});
150  });
151  }
152
153  Matrix compMatrix(component.size(), component.size(), triplets);
154  initialVectors[compIdx] = Vector(component.size());
155  rhsVectors[compIdx] = Vector(component.size());
156  lamgSetup.setup(compMatrix, compHierarchies[compIdx]);
157  compSolvers.push_back(SolverLamg<Matrix>(compHierarchies[compIdx], smoother));
158  LAMGSolverStatus status;
159  status.desiredResidualReduction = this->tolerance * component.size() / G.numberOfNodes();
160  compStati[compIdx] = status;
161
162  compIdx++;
163  }
164
165  validSetup = true;
166  }
167 }
168
169 template<class Matrix>
170 SolverStatus Lamg<Matrix>::solve(const Vector& rhs, Vector& result, count maxConvergenceTime, count maxIterations) {
171  if (!validSetup || result.getDimension() != laplacianMatrix.numberOfColumns()
172  || rhs.getDimension() != laplacianMatrix.numberOfRows()) {
173  throw std::runtime_error("No or wrong matrix is setup for given vectors.");
174  }
175
176  SolverStatus status;
177
178  if (numComponents == 1) {
179  LAMGSolverStatus stat;
180  stat.desiredResidualReduction = this->tolerance * rhs.length() / (laplacianMatrix * result - rhs).length();
181  stat.maxIters = maxIterations;
182  stat.maxConvergenceTime = maxConvergenceTime;
183  compSolvers[0].solve(result, rhs, stat);
184
185  status.residual = stat.residual;
186  status.numIters = stat.numIters;
187  status.converged = stat.converged;
188  } else {
189  // solve on every component
190  count maxIters = 0;
191  for (index i = 0; i < components.size(); ++i) {
192  for (auto element : components[i]) {
193  initialVectors[i][graph2Components[element]] = result[element];
194  rhsVectors[i][graph2Components[element]] = rhs[element];
195  }
196
197  double resReduction = this->tolerance * rhsVectors[i].length() / (compHierarchies[i].at(0).getLaplacian() * initialVectors[i] - rhsVectors[i]).length();
198  compStati[i].desiredResidualReduction = resReduction * components[i].size() / laplacianMatrix.numberOfRows();
199  compStati[i].maxIters = maxIterations;
200  compStati[i].maxConvergenceTime = maxConvergenceTime;
201  compSolvers[i].solve(initialVectors[i], rhsVectors[i], compStati[i]);
202
203  for (auto element : components[i]) { // write solution back to result
204  result[element] = initialVectors[i][graph2Components[element]];
205  }
206
207  maxIters = std::max(maxIters, compStati[i].numIters);
208  }
209
210  status.residual = (rhs - laplacianMatrix * result).length();
211  status.converged = status.residual <= this->tolerance;
212  status.numIters = maxIters;
213  }
214
215  return status;
216 }
217
218 template<class Matrix>
219 void Lamg<Matrix>::parallelSolve(const std::vector<Vector>& rhs, std::vector<Vector>& results, count maxConvergenceTime, count maxIterations) {
220  if (numComponents == 1) {
221  assert(rhs.size() == results.size());
223  if (compSolvers.size() != numThreads) {
224  compSolvers.clear();
225
226  for (index i = 0; i < (index) numThreads; ++i) {
227  compSolvers.push_back(SolverLamg<Matrix>(compHierarchies[0], smoother));
228  }
229  }
230
231  bool nested = omp_get_nested();
232  if (nested) omp_set_nested(false);
233
234 #pragma omp parallel for
235  for (index i = 0; i < rhs.size(); ++i) {
237  LAMGSolverStatus stat;
238  stat.desiredResidualReduction = this->tolerance * rhs[i].length() / (laplacianMatrix * results[i] - rhs[i]).length();
239  stat.maxIters = maxIterations;
240  stat.maxConvergenceTime = maxConvergenceTime;
242  }
243
244  if (nested) omp_set_nested(true);
245  }
246 }
247
248 } /* namespace NetworKit */
249
250 #endif /* NETWORKIT_CPP_NUMERICS_LAMG_LAMG_H_ */
Lamg(const double tolerance=1e-6)
Construct a solver with the given tolerance.
Definition: Lamg.h:55
bool converged
Definition: SolverLamg.h:31
void parallelSolve(const std::vector< Vector > &rhs, std::vector< Vector > &results, count maxConvergenceTime=5 *60 *1000, count maxIterations=std::numeric_limits< count >::max())
Compute the results for the matrix currently setup and the right-hand sides rhs.
Definition: Lamg.h:219
SolverStatus solve(const Vector &rhs, Vector &result, count maxConvergenceTime=5 *60 *1000, count maxIterations=std::numeric_limits< count >::max())
Computes the result for the matrix currently setup and the right-hand side rhs.
Definition: Lamg.h:170
count getDimension() const
Definition: Vector.h:74
double desiredResidualReduction
Definition: SolverLamg.h:24
bool converged
Definition: LinearSolver.h:22
count maxIters
Definition: SolverLamg.h:22
Determines the connected components of an undirected graph.
Definition: ParallelConnectedComponents.h:22
uint64_t index
Typedefs.
Definition: Globals.h:20
void setup(const Matrix &laplacianMatrix)
Compute the multigrid hierarchy for the given Laplacian matrix laplacianMatrix.
Definition: Lamg.h:116
double residual
Definition: SolverLamg.h:30
count numIters
Definition: LinearSolver.h:20
void run()
This method determines the connected components for the graph g.
Definition: ParallelConnectedComponents.cpp:20
void setupConnected(const Matrix &laplacianMatrix)
Compute the multigrid hierarchy for te given Laplacian matrix laplacianMatrix.
Definition: Lamg.h:109
count numberOfComponents()
This method returns the number of connected components.
Definition: ParallelConnectedComponents.cpp:174
Describes the status of a LinearSolver after the solver finished.
Definition: LinearSolver.h:19
count maxConvergenceTime
Definition: SolverLamg.h:23
count numIters
Definition: SolverLamg.h:29
uint64_t count
Definition: Globals.h:21
std::vector< std::vector< count > > Matrix
Definition: DynamicNMIDistance.h:16
Abstract base class for solvers that solve linear systems.
Definition: LinearSolver.h:29
index node
Definition: Globals.h:23
Implements the setup phase of LAMG (Lean Algebraic Multigrid by Livne et al.).
Definition: MultiLevelSetup.h:24
Implements the solve phase of LAMG (Lean Algebraic Multigrid by Livne et al.).
Definition: SolverLamg.h:40
double tolerance
Definition: LinearSolver.h:31
double length() const
Calculates and returns the Euclidean length of this vector.
Definition: Vector.cpp:34
Partition getPartition()
Return a Partition that represents the components.
Definition: ParallelConnectedComponents.cpp:170
The Vector class represents a basic vector with double coefficients.
Definition: Vector.h:25
std::set< std::set< index > > getSubsets() const
Definition: Partition.cpp:156
NetworKit::Graph matrixToGraph(const Matrix &matrix)
Interprets the matrix as adjacency matrix of a graph.
Definition: MatrixTools.h:107
double residual
Definition: LinearSolver.h:21
count numberOfNodes() const
Return the number of nodes in the graph.
Definition: Graph.h:706
Represents the interface to the Lean Algebraic Multigrid (LAMG) graph Laplacian linear solver by Oren...
Definition: Lamg.h:30
Status parameters of the solver.
Definition: SolverLamg.h:20
void forNeighborsOf(node u, L handle) const
Iterate over all neighbors of a node and call handle (lamdba closure).
Definition: Graph.h:1378
Implementation of the Gauss-Seidel smoother.
Definition: GaussSeidelRelaxation.h:20
A graph (with optional weights) and parallel iterator methods.
Definition: Graph.h:79
~Lamg()=default
Default destructor.
double edgeweight
Definition: Globals.h:24