All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SolverLamg.h
Go to the documentation of this file.
1 /*
2  * SolverLamg.h
3  *
4  * Created on: 12.01.2015
5  * Author: Michael
6  */
7 
8 #ifndef SOLVERLAMG_H_
9 #define SOLVERLAMG_H_
10 
11 #include "LevelHierarchy.h"
12 #include "../Smoother.h"
13 #include "../../algebraic/DenseMatrix.h"
14 
15 namespace NetworKit {
16 
21  // in
22  count maxIters = std::numeric_limits<count>::max(); // maximum number of iterations
23  count maxConvergenceTime = std::numeric_limits<count>::max(); // maximum time in milliseconds spent to solve the system
24  double desiredResidualReduction = 1e-8; // desired reduction of the initial residual (finalResidual <= desiredResReduction * initialResidual)
25  count numPreSmoothIters = 1; // number of pre smoothing iterations
26  count numPostSmoothIters = 2; // number of post smoothing iterations
27 
28  // out
29  count numIters; // number of iterations needed during solve phase
30  double residual; // absolute final residual
31  bool converged; // flag of conversion status
32  std::vector<double> residualHistory; // history of absolute residuals
33 };
34 
39 template<class Matrix>
40 class SolverLamg {
41 private:
42  LevelHierarchy<Matrix> &hierarchy;
43  const Smoother<Matrix> &smoother;
44 
45  // data structures for iterate recombination
46  std::vector<std::vector<Vector>> history;
47  std::vector<std::vector<Vector>> rHistory;
48  std::vector<index> latestIterate;
49  std::vector<count> numActiveIterates;
50 
51  // bStages for Elimination Levels
52  std::vector<std::vector<Vector>> bStages;
53 
54  void solveCycle(Vector& x, const Vector& b, int finest, LAMGSolverStatus& status);
55  void cycle(Vector& x, const Vector& b, int finest, int coarsest, std::vector<count>& numVisits, std::vector<Vector>& X, std::vector<Vector>& B, const LAMGSolverStatus& status);
56  void multigridCycle(index level, Vector& xf, const Vector& bf);
57  void saveIterate(index level, const Vector& x, const Vector& r);
58  void clearHistory(index level);
59  void minRes(index level, Vector& x, const Vector& r);
60 
61 public:
68  SolverLamg(LevelHierarchy<Matrix>& hierarchy, const Smoother<Matrix>& smoother) : hierarchy(hierarchy), smoother(smoother), bStages(hierarchy.size(), std::vector<Vector>()) {}
69 
70  SolverLamg (const SolverLamg<Matrix>& other) = default;
71 
72  SolverLamg (SolverLamg<Matrix>&& other) = default;
73 
74  virtual ~SolverLamg() = default;
75 
76  SolverLamg& operator=(SolverLamg<Matrix>&& other) = default;
77 
78  SolverLamg& operator=(const SolverLamg<Matrix>& other) = default;
79 
88  void solve(Vector& x, const Vector& b, LAMGSolverStatus& status);
89 };
90 
91 template<class Matrix>
93  bStages = std::vector<std::vector<Vector>>(hierarchy.size(), std::vector<Vector>());
94  if (hierarchy.size() >= 2) {
95  Vector bc = b;
96  Vector xc = x;
97  int finest = 0;
98 
99  if (hierarchy.getType(1) == ELIMINATION) {
100  hierarchy.at(1).restrict(b, bc, bStages[1]);
101  if (hierarchy.at(1).getLaplacian().numberOfRows() == 1) {
102  x = 0.0;
103  return;
104  } else {
105  hierarchy.at(1).coarseType(x, xc);
106  finest = 1;
107  }
108  }
109  solveCycle(xc, bc, finest, status);
110 
111  if (finest == 1) { // interpolate from finest == ELIMINATION level back to actual finest level
112  hierarchy.at(1).interpolate(xc, x, bStages[1]);
113  } else {
114  x = xc;
115  }
116  } else {
117  solveCycle(x, b, 0, status);
118  }
119 
120  double residual = (b - hierarchy.at(0).getLaplacian() * x).length();
121  status.residual = residual;
122 }
123 
124 template<class Matrix>
125 void SolverLamg<Matrix>::solveCycle(Vector& x, const Vector& b, int finest, LAMGSolverStatus& status) {
126  Aux::Timer timer;
127  timer.start();
128 
129  // data structures for iterate recombination
130  history = std::vector<std::vector<Vector>>(hierarchy.size());
131  rHistory = std::vector<std::vector<Vector>>(hierarchy.size());
132  latestIterate = std::vector<index>(hierarchy.size(), 0);
133  numActiveIterates = std::vector<count>(hierarchy.size(), 0);
134  int coarsest = hierarchy.size() - 1;
135  std::vector<count> numVisits(coarsest);
136  std::vector<Vector> X(hierarchy.size());
137  std::vector<Vector> B(hierarchy.size());
138 
139  for (index i = 0; i < hierarchy.size(); ++i) {
140  history[i] = std::vector<Vector>(MAX_COMBINED_ITERATES, Vector(hierarchy.at(i).getNumberOfNodes()));
141  rHistory[i] = std::vector<Vector>(MAX_COMBINED_ITERATES, Vector(hierarchy.at(i).getNumberOfNodes()));
142  }
143 
144  Vector r = b - hierarchy.at(finest).getLaplacian() * x;
145  double residual = r.length();
146  double finalResidual = residual * status.desiredResidualReduction;
147  double bestResidual = std::numeric_limits<double>::max();
148 
149  count iterations = 0;
150  status.residualHistory.emplace_back(residual);
151  count noResReduction = 0;
152  while (residual > finalResidual && noResReduction < 5 && iterations < status.maxIters && timer.elapsedMilliseconds() <= status.maxConvergenceTime ) {
153  cycle(x, b, finest, coarsest, numVisits, X, B, status);
154  r = b - hierarchy.at(finest).getLaplacian() * x;
155  residual = r.length();
156  status.residualHistory.emplace_back(residual);
157  if (residual < bestResidual) {
158  noResReduction = 0;
159  bestResidual = residual;
160  } else {
161  ++noResReduction;
162  }
163  iterations++;
164  }
165 
166  timer.stop();
167 
168  status.numIters = iterations;
169  status.residual = r.length();
170  status.converged = r.length() <= finalResidual;
171 }
172 
173 template<class Matrix>
174 void SolverLamg<Matrix>::cycle(Vector& x, const Vector& b, int finest, int coarsest, std::vector<count>& numVisits, std::vector<Vector>& X, std::vector<Vector>& B, const LAMGSolverStatus& status) {
175  std::fill(numVisits.begin(), numVisits.end(), 0);
176  X[finest] = x;
177  B[finest] = b;
178 
179  int currLvl = finest;
180  int nextLvl = finest;
181  double maxVisits = 0.0;
182 
183  saveIterate(currLvl, X[currLvl], B[currLvl] - hierarchy.at(currLvl).getLaplacian() * X[currLvl]);
184  while (true) {
185  if (currLvl == coarsest) {
186  nextLvl = currLvl - 1;
187  if (currLvl == finest) { // finest level
188  X[currLvl] = smoother.relax(hierarchy.at(currLvl).getLaplacian(), B[currLvl], X[currLvl], status.numPreSmoothIters);
189  } else {
190  Vector bCoarse(B[currLvl].getDimension()+1, 0.0);
191  for (index i = 0; i < B[currLvl].getDimension(); ++i) {
192  bCoarse[i] = B[currLvl][i];
193  }
194 
195  Vector xCoarse = DenseMatrix::LUSolve(hierarchy.getCoarseMatrix(), bCoarse);
196  for (index i = 0; i < X[currLvl].getDimension(); ++i) {
197  X[currLvl][i] = xCoarse[i];
198  }
199  }
200  } else {
201  if (currLvl == finest) {
202  maxVisits = 1.0;
203  } else {
204  maxVisits = hierarchy.cycleIndex(currLvl) * numVisits[currLvl-1];
205  }
206 
207  if (numVisits[currLvl] < maxVisits) {
208  nextLvl = currLvl + 1;
209  } else {
210  nextLvl = currLvl - 1;
211  }
212  }
213 
214  if (nextLvl < finest) break;
215 
216  if (nextLvl > currLvl) { // preProcess
217  numVisits[currLvl]++;
218 
219  if (hierarchy.getType(nextLvl) != ELIMINATION) {
220  X[currLvl] = smoother.relax(hierarchy.at(currLvl).getLaplacian(), B[currLvl], X[currLvl], status.numPreSmoothIters);
221  }
222 
223  if (hierarchy.getType(nextLvl) == ELIMINATION) {
224  hierarchy.at(nextLvl).restrict(B[currLvl], B[nextLvl], bStages[nextLvl]);
225  } else {
226  hierarchy.at(nextLvl).restrict(B[currLvl] - hierarchy.at(currLvl).getLaplacian() * X[currLvl], B[nextLvl]);
227  }
228 
229  hierarchy.at(nextLvl).coarseType(X[currLvl], X[nextLvl]);
230 
231  clearHistory(nextLvl);
232  } else { // postProcess
233  if (currLvl == coarsest || hierarchy.getType(currLvl+1) != ELIMINATION) {
234  minRes(currLvl, X[currLvl], B[currLvl] - hierarchy.at(currLvl).getLaplacian() * X[currLvl]);
235  }
236 
237  if (nextLvl > finest) {
238  saveIterate(nextLvl, X[nextLvl], B[nextLvl] - hierarchy.at(nextLvl).getLaplacian() * X[nextLvl]);
239  }
240 
241  if (hierarchy.getType(currLvl) == ELIMINATION) {
242  hierarchy.at(currLvl).interpolate(X[currLvl], X[nextLvl], bStages[currLvl]);
243  } else {
244  Vector xf = X[nextLvl];
245  hierarchy.at(currLvl).interpolate(X[currLvl], xf);
246  X[nextLvl] += xf;
247  }
248 
249  if (hierarchy.getType(currLvl) != ELIMINATION) {
250  X[nextLvl] = smoother.relax(hierarchy.at(nextLvl).getLaplacian(), B[nextLvl], X[nextLvl], status.numPostSmoothIters);
251  }
252 
253  }
254 
255  currLvl = nextLvl;
256  } // while
257 
258  // post-cycle finest
259  if ((int64_t) hierarchy.size() > finest + 1 && hierarchy.getType(finest+1) != ELIMINATION) { // do an iterate recombination on calculated solutions
260  minRes(finest, X[finest], B[finest] - hierarchy.at(finest).getLaplacian() * X[finest]);
261  }
262 
263 
264  X[finest] -= X[finest].mean();
265  x = X[finest];
266 }
267 
268 template<class Matrix>
269 void SolverLamg<Matrix>::saveIterate(index level, const Vector& x, const Vector& r) {
270  // update latest pointer
271  index i = latestIterate[level];
272  latestIterate[level] = (i+1) % MAX_COMBINED_ITERATES;
273 
274 
275  // update numIterates
276  if (numActiveIterates[level] < MAX_COMBINED_ITERATES) {
277  numActiveIterates[level]++;
278  }
279 
280  // update history array
281  history[level][i] = x;
282  rHistory[level][i] = r;
283 }
284 
285 template<class Matrix>
286 void SolverLamg<Matrix>::clearHistory(index level) {
287  latestIterate[level] = 0;
288  numActiveIterates[level] = 0;
289 }
290 
291 template<class Matrix>
292 void SolverLamg<Matrix>::minRes(index level, Vector& x, const Vector& r) {
293  if (numActiveIterates[level] > 0) {
294  count n = numActiveIterates[level];
295 
296  std::vector<index> ARowIdx(r.getDimension()+1);
297  std::vector<index> ERowIdx(r.getDimension()+1);
298 
299 #pragma omp parallel for
300  for (index i = 0; i < r.getDimension(); ++i) {
301  for (index k = 0; k < n; ++k) {
302  double AEvalue = r[i] - rHistory[level][k][i];
303  if (std::fabs(AEvalue) > 1e-25) {
304  ++ARowIdx[i+1];
305  }
306 
307  double Eval = history[level][k][i] - x[i];
308  if (std::fabs(Eval) > 1e-25) {
309  ++ERowIdx[i+1];
310  }
311  }
312  }
313 
314  for (index i = 0; i < r.getDimension(); ++i) {
315  ARowIdx[i+1] += ARowIdx[i];
316  ERowIdx[i+1] += ERowIdx[i];
317  }
318 
319  std::vector<index> AColumnIdx(ARowIdx[r.getDimension()]);
320  std::vector<double> ANonZeros(ARowIdx[r.getDimension()]);
321 
322  std::vector<index> EColumnIdx(ERowIdx[r.getDimension()]);
323  std::vector<double> ENonZeros(ERowIdx[r.getDimension()]);
324 
325 #pragma omp parallel for
326  for (index i = 0; i < r.getDimension(); ++i) {
327  for (index k = 0, aIdx = ARowIdx[i], eIdx = ERowIdx[i]; k < n; ++k) {
328  double AEvalue = r[i] - rHistory[level][k][i];
329  if (std::fabs(AEvalue) > 1e-25) {
330  AColumnIdx[aIdx] = k;
331  ANonZeros[aIdx] = AEvalue;
332  ++aIdx;
333  }
334 
335  double Eval = history[level][k][i] - x[i];
336  if (std::fabs(Eval) > 1e-25) {
337  EColumnIdx[eIdx] = k;
338  ENonZeros[eIdx] = Eval;
339  ++eIdx;
340  }
341  }
342  }
343 
344  CSRMatrix AE(r.getDimension(), n, ARowIdx, AColumnIdx, ANonZeros, 0.0, true);
345  CSRMatrix E(r.getDimension(), n, ERowIdx, EColumnIdx, ENonZeros, 0.0, true);
346 
347  Vector alpha = smoother.relax(CSRMatrix::mTmMultiply(AE, AE), CSRMatrix::mTvMultiply(AE, r), Vector(n, 0.0), 10);
348  x += E * alpha;
349  }
350 
351 }
352 
353 } /* namespace NetworKit */
354 
355 #endif /* SOLVERLAMG_H_ */
SolverLamg(LevelHierarchy< Matrix > &hierarchy, const Smoother< Matrix > &smoother)
Constructs a new solver instance for the specified hierarchy.
Definition: SolverLamg.h:68
static Vector mTvMultiply(const CSRMatrix &matrix, const Vector &vector)
Computes matrix^T * vector.
Definition: CSRMatrix.cpp:492
bool converged
Definition: SolverLamg.h:31
count numPreSmoothIters
Definition: SolverLamg.h:25
std::vector< double > residualHistory
Definition: SolverLamg.h:32
double & at(const index idx)
Returns a reference to the element at index idx.
Definition: Vector.h:128
double desiredResidualReduction
Definition: SolverLamg.h:24
A timer for running time measurements.
Definition: Timer.h:25
count maxIters
Definition: SolverLamg.h:22
static Vector LUSolve(const DenseMatrix &LU, const Vector &b)
Computes the solution vector x to the system LU * x = b where LU is a matrix decomposed into L and U...
Definition: DenseMatrix.cpp:213
uint64_t index
Typedefs.
Definition: Globals.h:20
count numPostSmoothIters
Definition: SolverLamg.h:26
double residual
Definition: SolverLamg.h:30
Definition: Level.h:16
count maxConvergenceTime
Definition: SolverLamg.h:23
constexpr count MAX_COMBINED_ITERATES
Definition: LAMGSettings.h:36
count numIters
Definition: SolverLamg.h:29
uint64_t count
Definition: Globals.h:21
Implements the solve phase of LAMG (Lean Algebraic Multigrid by Livne et al.).
Definition: SolverLamg.h:40
virtual ~SolverLamg()=default
SolverLamg & operator=(SolverLamg< Matrix > &&other)=default
double length() const
Calculates and returns the Euclidean length of this vector.
Definition: Vector.cpp:34
The Vector class represents a basic vector with double coefficients.
Definition: Vector.h:25
virtual uint64_t elapsedMilliseconds() const
The number of milliseconds since the current time that the Timer object was created.
Definition: Timer.cpp:43
virtual my_steady_clock::time_point stop()
Stops the clock permanently for the instance of the Timer.
Definition: Timer.cpp:21
void solve(Vector &x, const Vector &b, LAMGSolverStatus &status)
Solves the system A*x = b for the given initial x and right-hand side b.
Definition: SolverLamg.h:92
Status parameters of the solver.
Definition: SolverLamg.h:20
Abstract base class of a smoother.
Definition: Smoother.h:24
static CSRMatrix mTmMultiply(const CSRMatrix &A, const CSRMatrix &B)
Computes A^T * B.
Definition: CSRMatrix.cpp:431
virtual my_steady_clock::time_point start()
Start the clock.
Definition: Timer.cpp:15
Definition: LevelHierarchy.h:24