跳转至

共轭梯度法

背景介绍

在各类数值模拟计算问题中,求解线性方程是一个非常常见且重要的步骤。比如在流体数值计算(CFD)中,我们常常需要求解一个关于流体压力的泊松类型问题:

\frac{\delta t}{\rho}\nabla p = \nabla \cdot \vec {u}

其中方程的左边是一个压力场 p 的二阶导数,而右边是速度场(矢量)的散度。通常我们想要在已知速度场的情况下,求解其对应的压力场 p ,在特定的网格设定下,我们可以使用有限差分的方法我们可以将上面的微分方程转化为线性方程组。比如,下图演示了一个流体数值模拟中常用的,将压力存储在网格中心,将速度存储在网格界面上的离散形式:

其对应的离散化方程式是:

\frac{\delta t}{\rho}\left[\frac{1}{dx}\cdot \left(\frac{p_{i+1,j}-p_{i,j}}{dx} - \frac{p_{i,j}-p_{i-1,j}}{dx} \right) + \frac{1}{dy} \cdot \left( \frac{p_{i,j+1}-p_{i,j}}{dy} - \frac{p_{i,j}-p_{i, j-1}}{dy} \right) \right] = \frac{u_{i,j}-u_{i-1,j}}{dx} + \frac{v_{i,j}-v_{i,j-1}}{dy}

其中 ij 都是 [1, n] 的整数,我们可以对整个计算区域的所有格子列出类似的方程共 n×n 个。

对上式稍作整理可以发现当格子大小 dx=dy 时,方程左边可以简化成

\frac{\delta t}{\rho dx^2}(4p_{i,j}-p_{i+1,j}-p_{i-1,j}-p_{i,j-1}-p_{i,j+1})

其中, \frac{\delta t}{\rho dx^2} 是一个常数系数,为了简化问题我们可以直接将其省略。于是,我们可以将 p 的系数存入一个矩阵 A ,而将 p 保存在一个向量中,上面的表达式就变成了如下的一个矩阵和向量的乘:

于是这个问题就被转化为了求解线性方程组的问题。在各类数值仿真中,这样的例子还有很多很多。

共轭梯度法(Conjugate-gradient,以下简称为 CG 方法)是一种求解系数矩阵为对称正定矩阵的线性方程组的迭代解法。假设我们要求解如下的线性系统:

Ax=b

其中 A 是一个对称且正定,大小为 n×n 的实数矩阵, x 为待求解的大小为 n×1 的未知数向量。

整个共轭梯度法的过程可以用伪代码简单描述如下:

其中 x_k​ 是迭代求解过程中第 k 个循环时的解, r_k​ 是求解过程中在第 k 个迭代中的残差:

r_k ​= b − Ax_k

在循环迭代到 r_k​ 的模小于规定值时即认为求解完成。

输入输出

输入格式

本题中,约定系数矩阵 A 是一个对角线元素为 4.0,其余邻居对应元素为 -1.0 的五对角矩阵(five-diagonal matrices,除对角线和相邻元素外其他元素为 0)。

b 对应的是一个 n^2\times 1 的向量,共 20 组测试点,是二进制文件,小端法,储存着 n\times n 个 float32(4Byte)。命名规则为 b_case_n_k.bin ,如 b_1_256_1.bin ,输入文件放置在程序运行的相同文件夹下。

其中,case 是后文评分标准中的测试点编号,k 表示当前测试点的第 k 个数据点。程序运行时,需读入所有数据点,进行处理后输出。

输出格式

输出结果到 20 个二进制文件,小端法。命名规则为 ans_case_n_k.bin ,数据类型为 4byte 浮点数。

提交说明

支持pythonCC++CUDA 实现。提交一个tar 包,其中有start.sh,负责编译或者环境配置。

对于python,评测环境为module load pyenv。评测节点有网络连接,如需要更多的环境包,请在start.sh 中自行安装。

对于CC++的依赖,评测机与算力平台环境一致,可以使用sudo apt install 进行安装(脚本并非跑在 root 下)

对于CUDA,评测机环境为 cuda 11.8,在SCOW 平台module load nvtools 即可使用,或使用我们的GPU 类型算力。

对于其他任何语言,只要能编译、运行即可提交,可以详细咨询组委会

评测环境

8 Core 32G, A100 * 1(与选手可申请的环境相同,算力平台中的GPU 类型)

评测说明

本题暂时不支持实时评测,由组委会进行评测并审核。视提交时间,可能会有 10 分钟到 4 小时不等的延迟。

评分标准

将在规定的硬件平台上测算选手的程序,分别对于如下大小规模的问题进行 5 次运行,取平均耗时作为最终结果:

我们将把选手时间减去IO 时间和CUDA 初始化时间,作为耗时。

测试点 未知数数量 重复次数
1 256^2 5
2 512^2 5
3 1024^2 5
4 2048^2 5

我们提供了一个Taichi 实现作为例子,选手可以参考。

$ python3 main.py

附件

python 样例

下载cg.zip