All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConjugateGradient.h
Go to the documentation of this file.
1 /*
2  * ConjugateGradient.h
3  *
4  * Created on: 15.06.2014
5  * Author: Daniel Hoske and Michael Wegner
6  */
7 
8 #ifndef CONJUGATE_GRADIENT_H_
9 #define CONJUGATE_GRADIENT_H_
10 
11 #include <cstdint>
12 #include <utility>
13 
14 #include "LinearSolver.h"
15 #include "../algebraic/Vector.h"
16 #include "../algebraic/CSRMatrix.h"
17 
18 namespace NetworKit {
19 
24 template<class Matrix, class Preconditioner>
25 class ConjugateGradient : public LinearSolver<Matrix> {
26 public:
28 
29  void setup(const Matrix& matrix) {
30  this->matrix = matrix;
31  precond = Preconditioner(matrix);
32  }
33 
34  void setupConnected(const Matrix& matrix) {
35  this->matrix = matrix;
36  precond = Preconditioner(matrix);
37  }
38 
50  SolverStatus solve(const Vector& rhs, Vector& result, count maxConvergenceTime = 5 * 60 * 1000, count maxIterations = std::numeric_limits<count>::max());
51 
59  void parallelSolve(const std::vector<Vector>& rhs, std::vector<Vector>& results, count maxConvergenceTime = 5 * 60 * 1000, count maxIterations = std::numeric_limits<count>::max());
60 
61 private:
62  Matrix matrix;
63  Preconditioner precond;
64 
65 };
66 
67 template<class Matrix, class Preconditioner>
68 SolverStatus ConjugateGradient<Matrix, Preconditioner>::solve(const Vector& rhs, Vector& result, count maxConvergenceTime, count maxIterations) {
69  assert(matrix.numberOfRows() == rhs.getDimension());
70 
71  // Absolute residual to achieve
72  double sqr_desired_residual = this->tolerance * this->tolerance * (rhs.length() * rhs.length());
73 
74  // Main loop. See: http://en.wikipedia.org/wiki/Conjugate_gradient_method#The_resulting_algorithm
75  Vector residual_dir = rhs - matrix*result;
76  Vector conjugate_dir = precond.rhs(residual_dir);
77  double sqr_residual = Vector::innerProduct(residual_dir, residual_dir);
78  double sqr_residual_precond = Vector::innerProduct(residual_dir, conjugate_dir);
79 
80  count niters = 0;
81  Vector tmp, residual_precond;
82  while (sqr_residual > sqr_desired_residual) {
83  niters++;
84  if (niters > maxIterations) {
85  break;
86  }
87 
88  tmp = matrix * conjugate_dir;
89  double step = sqr_residual_precond / Vector::innerProduct(conjugate_dir, tmp);
90  result += step * conjugate_dir;
91  residual_dir -= step * tmp;
92  sqr_residual = Vector::innerProduct(residual_dir, residual_dir);
93 
94  residual_precond = precond.rhs(residual_dir);
95  double new_sqr_residual_precond = Vector::innerProduct(residual_dir, residual_precond);
96  conjugate_dir = (new_sqr_residual_precond / sqr_residual_precond) * conjugate_dir + residual_precond;
97  sqr_residual_precond = new_sqr_residual_precond;
98  }
99 
100  SolverStatus status;
101  status.numIters = niters;
102  status.residual = (rhs - matrix*result).length();
103  status.converged = status.residual / rhs.length() <= this->tolerance;
104 
105  return status;
106 }
107 
108 template<class Matrix, class Preconditioner>
109 void ConjugateGradient<Matrix, Preconditioner>::parallelSolve(const std::vector<Vector>& rhs, std::vector<Vector>& results, count maxConvergenceTime, count maxIterations) {
110 #pragma omp parallel for
111  for (index i = 0; i < rhs.size(); ++i) {
112  this->solve(rhs[i], results[i], maxConvergenceTime, maxIterations);
113  }
114 }
115 
116 
117 
118 
119 
120 
121 } /* namespace NetworKit */
122 
123 #endif /* CONJUGATE_GRADIENT_H_ */
SolverStatus solve(const Vector &rhs, Vector &result, count maxConvergenceTime=5 *60 *1000, count maxIterations=std::numeric_limits< count >::max())
Solves the linear system using the conjugate gradient method with a given preconditioner and with in...
Definition: ConjugateGradient.h:68
ConjugateGradient(double tolerance=1e-5)
Definition: ConjugateGradient.h:27
static double innerProduct(const Vector &v1, const Vector &v2)
Computes the inner product (dot product) of the vectors v1 and v2.
Definition: Vector.cpp:61
count getDimension() const
Definition: Vector.h:74
void parallelSolve(const std::vector< Vector > &rhs, std::vector< Vector > &results, count maxConvergenceTime=5 *60 *1000, count maxIterations=std::numeric_limits< count >::max())
Solves the linear systems in parallel.
Definition: ConjugateGradient.h:109
bool converged
Definition: LinearSolver.h:22
uint64_t index
Typedefs.
Definition: Globals.h:20
void setup(const Matrix &matrix)
Sets the solver up for the specified matrix.
Definition: ConjugateGradient.h:29
count numIters
Definition: LinearSolver.h:20
Implementation of Conjugate Gradient.
Definition: ConjugateGradient.h:25
void setupConnected(const Matrix &matrix)
Sets the solver up for the specified matrix where the underlying graph has to be connected.
Definition: ConjugateGradient.h:34
Describes the status of a LinearSolver after the solver finished.
Definition: LinearSolver.h:19
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
double tolerance
Definition: LinearSolver.h:31
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
double residual
Definition: LinearSolver.h:21