Editorial StaffThe Conjugate Gradient (CG) method is an iterative algorithm used for solving linear systems of equations. It is particularly effective for symmetric positive definite matrices.
Here’s an example implementation of the Conjugate Gradient method in C++:
The conjugateGradient function implements the Conjugate Gradient method to solve a system of linear equations Ax = b, where A is a symmetric positive definite matrix and b is the right-hand side vector. The function returns the solution vector x.
The dotProduct function calculates the dot product of two vectors and the matrixVectorMultiply function performs matrix-vector multiplication.
In the main function, an example usage is shown, where a 2×2 matrix A and a 2-dimensional vector b are defined. The conjugateGradient the function is called to solve the system, and the resulting solution vector x is printed.
Please note that this is a basic implementation and may not handle all possible edge cases or provide optimal performance.
#include <iostream>
#include <vector>
// Function to calculate the dot product of two vectors
double dotProduct(const std::vector<double>& a, const std::vector<double>& b) {
double result = 0.0;
for (size_t i = 0; i < a.size(); ++i) {
result += a[i] * b[i];
}
return result;
}
// Function to perform matrix-vector multiplication
std::vector<double> matrixVectorMultiply(const std::vector<std::vector<double>>& A, const std::vector<double>& x) {
size_t n = A.size();
std::vector<double> result(n, 0.0);
for (size_t i = 0; i < n; ++i) {
for (size_t j = 0; j < n; ++j) {
result[i] += A[i][j] * x[j];
}
}
return result;
}
// Conjugate Gradient method
std::vector<double> conjugateGradient(const std::vector<std::vector<double>>& A, const std::vector<double>& b) {
size_t n = A.size();
std::vector<double> x(n, 0.0); // Initial guess for the solution
std::vector<double> r = b; // Residual vector
std::vector<double> p = r; // Search direction vector
for (size_t k = 0; k < n; ++k) {
std::vector<double> Ap = matrixVectorMultiply(A, p);
double alpha = dotProduct(r, r) / dotProduct(p, Ap);
for (size_t i = 0; i < n; ++i) {
x[i] += alpha * p[i];
r[i] -= alpha * Ap[i];
}
double beta = dotProduct(r, r) / dotProduct(p, p);
for (size_t i = 0; i < n; ++i) {
p[i] = r[i] + beta * p[i];
}
}
return x;
}
int main() {
// Example usage
std::vector<std::vector<double>> A = { {4, 1}, {1, 3} }; // Symmetric positive definite matrix
std::vector<double> b = { 1, 2 }; // Right-hand side vector
std::vector<double> x = conjugateGradient(A, b);
std::cout << "Solution: ";
for (const auto& value : x) {
std::cout << value << " ";
}
std::cout << std::endl;
return 0;
}
