C++ 多线程共轭梯度法求解线性方程组

共轭梯度法是一种高效的优化算法，可以用于求解线性方程组或最小化二次函数的问题。多线程可以加速算法的运行速度，特别是对于大规模的问题。以下是 C++ 实现共轭梯度法多线程求解的示例代码：

#include <iostream>
#include <cmath>
#include <thread>
#include <vector>

using namespace std;

const int MAX_ITER = 1000; // 最大迭代次数
const double EPS = 1e-6; // 精度
const int THREAD_NUM = 4; // 线程数量

// 定义向量类
class Vector {
private:
    int n;
    double* data;

public:
    Vector(int n) : n(n) {
        data = new double[n];
        for (int i = 0; i < n; i++) {
            data[i] = 0;
        }
    }

    Vector(const Vector& other) : n(other.n) {
        data = new double[n];
        for (int i = 0; i < n; i++) {
            data[i] = other.data[i];
        }
    }

    ~Vector() {
        delete[] data;
    }

    int size() const {
        return n;
    }

    double& operator[](int i) {
        return data[i];
    }

    const double& operator[](int i) const {
        return data[i];
    }

    double dot(const Vector& other) const {
        double result = 0;
        for (int i = 0; i < n; i++) {
            result += data[i] * other.data[i];
        }
        return result;
    }

    Vector operator+(const Vector& other) const {
        Vector result(n);
        for (int i = 0; i < n; i++) {
            result.data[i] = data[i] + other.data[i];
        }
        return result;
    }

    Vector operator-(const Vector& other) const {
        Vector result(n);
        for (int i = 0; i < n; i++) {
            result.data[i] = data[i] - other.data[i];
        }
        return result;
    }

    Vector operator*(double scalar) const {
        Vector result(n);
        for (int i = 0; i < n; i++) {
            result.data[i] = data[i] * scalar;
        }
        return result;
    }

    Vector& operator+=(const Vector& other) {
        for (int i = 0; i < n; i++) {
            data[i] += other.data[i];
        }
        return *this;
    }

    Vector& operator-=(const Vector& other) {
        for (int i = 0; i < n; i++) {
            data[i] -= other.data[i];
        }
        return *this;
    }

    Vector& operator*=(double scalar) {
        for (int i = 0; i < n; i++) {
            data[i] *= scalar;
        }
        return *this;
    }
};

// 定义矩阵类
class Matrix {
private:
    int n;
    double** data;

public:
    Matrix(int n) : n(n) {
        data = new double*[n];
        for (int i = 0; i < n; i++) {
            data[i] = new double[n];
            for (int j = 0; j < n; j++) {
                data[i][j] = 0;
            }
        }
    }

    Matrix(const Matrix& other) : n(other.n) {
        data = new double*[n];
        for (int i = 0; i < n; i++) {
            data[i] = new double[n];
            for (int j = 0; j < n; j++) {
                data[i][j] = other.data[i][j];
            }
        }
    }

    ~Matrix() {
        for (int i = 0; i < n; i++) {
            delete[] data[i];
        }
        delete[] data;
    }

    int size() const {
        return n;
    }

    double* operator[](int i) {
        return data[i];
    }

    const double* operator[](int i) const {
        return data[i];
    }

    Vector operator*(const Vector& x) const {
        Vector result(n);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                result[i] += data[i][j] * x[j];
            }
        }
        return result;
    }
};

// 计算向量范数的平方
double norm_squared(const Vector& x) {
    return x.dot(x);
}

// 计算向量的二范数
double norm(const Vector& x) {
    return sqrt(norm_squared(x));
}

// 计算矩阵的转置
Matrix transpose(const Matrix& A) {
    Matrix result(A.size());
    for (int i = 0; i < A.size(); i++) {
        for (int j = 0; j < A.size(); j++) {
            result[j][i] = A[i][j];
        }
    }
    return result;
}

// 计算共轭梯度法的迭代步骤
void cg_step(int start, int end, const Matrix& A, const Vector& b, Vector& x, Vector& r, Vector& p) {
    for (int k = start; k < end; k++) {
        double alpha = r.dot(r) / p.dot(A * p);
        x += p * alpha;
        Vector r_new = r - A * p * alpha;
        if (norm(r_new) < EPS) {
            break;
        }
        double beta = r_new.dot(r_new) / r.dot(r);
        p = r_new + p * beta;
        r = r_new;
    }
}

// 求解线性方程组 Ax=b
Vector solve(const Matrix& A, const Vector& b) {
    int n = A.size();
    Vector x(n);
    Vector r = b - A * x;
    Vector p = r;
    for (int i = 0; i < MAX_ITER; i++) {
        vector<thread> threads;
        for (int j = 0; j < THREAD_NUM; j++) {
            int start = j * (n / THREAD_NUM);
            int end = (j + 1) * (n / THREAD_NUM);
            threads.emplace_back(cg_step, start, end, ref(A), ref(b), ref(x), ref(r), ref(p));
        }
        for (auto& thread : threads) {
            thread.join();
        }
        if (norm(r) < EPS) {
            break;
        }
    }
    return x;
}

// 测试
int main() {
    Matrix A(3);
    A[0][0] = 4; A[0][1] = 1; A[0][2] = 1;
    A[1][0] = 1; A[1][1] = 3; A[1][2] = 1;
    A[2][0] = 1; A[2][1] = 1; A[2][2] = 2;

    Vector b(3);
    b[0] = 1; b[1] = 2; b[2] = 3;

    Vector x = solve(A, b);

    cout << 'Solution: ';
    for (int i = 0; i < x.size(); i++) {
        cout << x[i] << ' ';
    }
    cout << endl;

    return 0;
}

在以上代码中，向量和矩阵分别用 Vector 和 Matrix 类表示，其中 Vector 类包含了向量的基本运算，Matrix 类包含了矩阵的转置和向量乘法运算。共轭梯度法的求解过程被分成了多个线程，在每个线程中执行一部分迭代步骤，然后等待所有线程执行完毕。在求解线性方程组时，调用 solve 函数即可。在测试中，使用一个 3x3 的矩阵和一个 3 维的向量作为例子进行求解。