以下是科研实习第四周的进度
代码修正
在上周的演示中，主要使用了Jacobi迭代法对拉普拉斯方程进行了求解，但是经过测试，在较为复杂的边界条件上，Jacobi法迭代速度较慢甚至过慢，需要很长的时间才能得到解。
对此进行分析，对于这样的一个PDE，使用有限差分法我们很容易将其写成一个线性方程组，对于一个100 * 100 的拉普拉斯方程，其对应的线性方程组具有100*100*100*100的矩阵，有限差分法需要从边界向中心迭代，在开启并行的情况下，中心部分可能需要很多个Epoch才能发生变化。
因此对于大型矩阵，我们一方面可以采用迭代速度更快的G-S迭代法以及松弛系数在1-2的松弛迭代法，另一方面，考虑到拉普拉斯方程对应的系数矩阵是一个非常稀疏的矩阵，我们可以采用稀疏矩阵的求解法对其进行求解。
下面使用基于Givens旋转的QR分解对拉普拉斯方程进行求解（非稀疏矩阵特化，稀疏矩阵特化与另外一个同学讨论并进行实验，目前仍旧存在一些问题，因此求解较慢，使用规模较小的网格)
Unit Test
我们使用 Boost.ut 进行单元测试，具体测试代码如下所示
vector<double> x = solve_linear_system(A,b);
"Serial"_test = [&]{
for (int i = 0; i < 400; ++i) {
expect(fabs(x[i]-1<1e-6));
}
};
x = solve_linear_system_omp(A,b);
"OpenMP"_test = [&]{
for (int i = 0; i < 400; ++i) {
expect(fabs(x[i]-1<1e-6));
}
};
x = solve_linear_system_cu(A,b);
"CUDA"_test = [&]{
for (int i = 0; i < 400; ++i) {
expect(fabs(x[i]-1<1e-6));
}
};
输出结果为：
Suite 'global': all tests passed (1200 asserts in 3 tests)
注：由于并行度不高，CUDA的优化大多不是来源于并行，而是显卡本身对数组的操作速率高于CPU
注：CUDA运行的结果不稳定，原因可能是具有访址冲突，但是Debug过程中存在一定的修复问题，因此还未完全修复
Benchmark
我们使用NanoBench对代码进行性能测试，我们采用的方式是对一个代码连续运行20次取其平均效果，代码如下：
Bench().run("Serial",[&]{
vector<vector<double>> A_temp(A);
vector<double> b_temp(b);
solve_linear_system(A_temp,b_temp);
});
Bench().run("OpenMP",[&]{
vector<vector<double>> A_temp(A);
vector<double> b_temp(b);
solve_linear_system_omp(A_temp,b_temp);
});
Bench().run("CUDA",[&]{
vector<vector<double>> A_temp(A);
vector<double> b_temp(b);
solve_linear_system_cu(A_temp,b_temp);
});
其输出结果如下：
| ns/op | op/s | err% | ins/op | cyc/op | IPC | bra/op | miss% | total | benchmark
|--------------------:|--------------------:|--------:|----------------:|----------------:|-------:|---------------:|--------:|----------:|:----------
| 2,454,934,332.00 | 0.41 | 0.1% |9,020,857,795.00 |5,620,531,896.00 | 1.605 |1,152,369,662.00 | 0.0% | 27.12 | `Serial`
| 1,411,725,552.00 | 0.71 | 0.7% |2,779,103,816.00 |3,232,946,976.00 | 0.860 | 500,965,126.00 | 3.5% | 15.48 | `OpenMP`
| 51,278,734.00 | 19.50 | 0.2% | 149,190,381.00 | 117,345,448.00 | 1.271 | 27,837,262.00 | 6.0% | 0.74 | `CUDA`