Go to the documentation of this file.
1 /*
3  *
4  * Created on: 15.06.2014
5  * Author: Daniel Hoske and Michael Wegner
6  */
7
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
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...
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.
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.
count numIters
Definition: LinearSolver.h:20
void setupConnected(const Matrix &matrix)
Sets the solver up for the specified matrix where the underlying graph has to be connected.
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