并行求解线性方程组：使用 Conjugate Gradient 算法的 C++ 示例

本文提供了一个使用 Conjugate Gradient (CG) 算法并行求解线性方程组的 C++ 示例代码：

#include <iostream>
#include <vector>
#include <cmath>
#include <chrono>
#include <omp.h>

using namespace std;
using namespace chrono;

const int N = 10000; // 矩阵大小
const int MAX_ITER = 5000; // 最大迭代次数
const double EPSILON = 1e-10; // 精度

vector<vector<double>> A(N, vector<double>(N)); // 系数矩阵
vector<double> b(N); // 常数矩阵
vector<double> x(N); // 未知数矩阵
vector<double> r(N); // 残差矩阵
vector<double> p(N); // 方向矩阵
vector<double> Ap(N); // AP矩阵

// 计算向量内积
double dot(const vector<double>& v1, const vector<double>& v2) {
    double res = 0;
    for (int i = 0; i < N; i++) {
        res += v1[i] * v2[i];
    }
    return res;
}

// 计算向量1加上向量2乘以常数c的结果
vector<double> add(const vector<double>& v1, const vector<double>& v2, double c) {
    vector<double> res(N);
    for (int i = 0; i < N; i++) {
        res[i] = v1[i] + v2[i] * c;
    }
    return res;
}

// 初始化矩阵
void init() {
    #pragma omp parallel for
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (i == j) {
                A[i][j] = 2.0;
            } else if (abs(i - j) == 1) {
                A[i][j] = -1.0;
            }
        }
        b[i] = 1.0 / (i + 1);
        x[i] = 0.0;
    }
}

// CG算法主体
void cg() {
    double alpha, beta, r_norm, r_norm_old;
    int iter = 0;

    r = add(b, A[0], -dot(x, A[0])); // 计算初始残差
    p = r; // 初始方向
    r_norm = dot(r, r); // 初始残差范数

    while (iter < MAX_ITER && sqrt(r_norm) > EPSILON) {
        Ap = vector<double>(N, 0.0);
        #pragma omp parallel for
        for (int i = 0; i < N; i++) {
            for (int j = 0; j < N; j++) {
                Ap[i] += A[i][j] * p[j];
            }
        }
        alpha = r_norm / dot(p, Ap); // 计算步长
        x = add(x, p, alpha); // 更新解向量
        r_norm_old = r_norm;
        r = add(r, Ap, -alpha); // 更新残差
        r_norm = dot(r, r); // 更新残差范数
        beta = r_norm / r_norm_old; // 计算方向系数
        p = add(r, p, beta); // 更新方向
        iter++;
    }

    cout << "CG迭代次数: " << iter << endl;
}

int main() {
    init();
    high_resolution_clock::time_point t1 = high_resolution_clock::now();
    cg();
    high_resolution_clock::time_point t2 = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(t2 - t1).count();
    cout << "总用时: " << duration << "us" << endl;
    return 0;
}

在这个示例中，我们首先定义了一个大小为N的系数矩阵A、常数矩阵b、未知数矩阵x、残差矩阵r、方向矩阵p和AP矩阵Ap。然后实现了dot和add函数用于计算向量内积和向量加法。在初始化函数init中，我们按照题目要求生成了系数矩阵和常数矩阵。在CG算法主体cg中，我们按照算法步骤依次计算了初值和各个迭代步骤中需要的各个向量，并更新了解向量、残差和方向。在主函数中，我们调用init和cg函数，并通过chrono库计算了程序总用时。

需要注意的是，在计算AP矩阵时，我们使用了OpenMP的并行for循环来加速计算过程。如果你的电脑上没有OpenMP库，可以通过以下命令安装：

sudo apt-get install libomp-dev

最后，将代码保存为cg.cpp并使用g++编译器编译并运行：

g++ -fopenmp cg.cpp -o cg
./cg