开源并行平台上共轭梯度法求解矩阵的实施步骤：步骤记录和关于并行矩阵求解的疑问
一、求解方程
以二扩散方程为例，采用有限体积法求解
- \mu \bigtriangledown ^{2} u = 0
- \int\limits_{\bigtriangleup V} \mu (\frac{\partial^{2}u }{\partial x^2} + \frac{\partial^{2}u }{\partial y^2}) dV= 0
（这里因为推公式的时候用的NS方程，将扩散项从右侧移到了左侧带上了负号，不影响单独计算扩散方程时的结果）
离散得到
- \mu_{e}A_{e}(\frac{\partial u}{\partial x})_{e}
+ \mu_{w}A_{w}(\frac{\partial u}{\partial x})_{w}
- \mu_{n}A_{n}(\frac{\partial u}{\partial x})_{n}
+ \mu_{s}A_{s}(\frac{\partial u}{\partial x})_{s} =0
由于采用的是同位网格，式中偏导项采用差分格式很容易得到离散形式
二、网格分区和矩阵并行的关联 及 并行Jacobi法求解的介绍
为了具体的写出矩阵，给出网格、边界条件和扩散系数的赋值：在4×4笛卡尔网格下，左右边界取定值（左边界 u = 100，右边界为 u = 500），上下边界为纽曼边界条件（边界与靠近边界格子的值相等），$\mu$ 取 100。
程序采用单核和两核运行对比：
图1 单进程下一个进程负责所有网格
图2 两进程下各进程负责的网格
从左到右，从上到下，对格子依次编号 1～16
在单进程下，组成的矩阵为：
{\tiny \begin{bmatrix}
400& -100& & & -100& & & & & & & & & & & \\
-100& 400& -100& & & -100& & & & & & & & & & \\
& -100& 400& -100& & & -100& & & & & & & & & \\
& & -100& 400& & & & -100& & & & & & & & \\
-100& & & & 400& -100& & & -100& & & & & & & \\
& -100& & & -100& 400& -100& & & -100& & & & & & \\
& & -100& & & -100& 400& -100& & & -100& & & & & \\
& & & -100& & & -100& 400& & & & -100& & & & \\
& & & & -100& & & & 400& -100& & & -100& & & \\
& & & & & -100& & & -100& 400& -100& & & -100& & \\
& & & & & & -100& & & -100& 400& -100& & & -100& \\
& & & & & & & -100& & & -100& 400& & & & -100\\
& & & & & & & & -100& & & & 400& -100& & \\
& & & & & & & & & -100& & & -100& 400& -100& \\
& & & & & & & & & & -100& & & -100& 400& -100\\
& & & & & & & & & & & -100& & & -100& 400
\end{bmatrix}
\begin{bmatrix}
u_1\\
u_2\\
u_3\\
u_4\\
u_5\\
u_6\\
u_7\\
u_8\\
u_9\\
u_{10}\\
u_{11}\\
u_{12}\\
u_{13}\\
u_{14}\\
u_{15}\\
u_{16}
\end{bmatrix}
=
\begin{bmatrix}
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000
\end{bmatrix}}
但在两进程下，由于网格分别在两进程，存在网格的“内边界”，只能访问内边界在上一迭代步的值：也就是说，黄色区域的格子在计算第 n 步时，将网格第三行紫色区域的格子的在第 n-1 迭代步计算得到的 u 作为边界，进行第 n 步的计算。
反映到矩阵上，就是如下这个变换，红色系数代表边界格子对该单元的影响（比如图2 网格第三行紫色格子作为黄色区域的边界，网格第二行黄色格子作为紫色区域的边界），因为利用了上一步的值，是已知量，移项到右侧，将矩阵根据网格分区划分为两个矩阵：
{\tiny \begin{bmatrix}
400& -100& & & -100& & & & & & & & & & & \\
-100& 400& -100& & & -100& & & & & & & & & & \\
& -100& 400& -100& & & -100& & & & & & & & & \\
& & -100& 400& & & & -100& & & & & & & & \\
-100& & & & 400& -100& & & {\color{Red} -100} & & & & & & & \\
& -100& & & -100& 400& -100& & & {\color{Red}-100}& & & & & & \\
& & -100& & & -100& 400& -100& & & {\color{Red}-100}& & & & & \\
& & & -100& & & -100& 400& & & & {\color{Red}-100}& & & & \\
& & & & {\color{Violet}-100}& & & & 400& -100& & & -100& & & \\
& & & & & {\color{Violet}-100}& & & -100& 400& -100& & & -100& & \\
& & & & & & {\color{Violet}-100}& & & -100& 400& -100& & & -100& \\
& & & & & & & {\color{Violet}-100}& & & -100& 400& & & & -100\\
& & & & & & & & -100& & & & 400& -100& & \\
& & & & & & & & & -100& & & -100& 400& -100& \\
& & & & & & & & & & -100& & & -100& 400& -100\\
& & & & & & & & & & & -100& & & -100& 400
\end{bmatrix}
\begin{bmatrix}
u_1\\
u_2\\
u_3\\
u_4\\
{\color{Violet} u_5} \\
{\color{Violet} u_6}\\
{\color{Violet} u_7}\\
{\color{Violet} u_8}\\
{\color{Red} u_9}\\
{\color{Red} u_{10}}\\
{\color{Red} u_{11}}\\
{\color{Red} u_{12}}\\
u_{13}\\
u_{14}\\
u_{15}\\
u_{16}
\end{bmatrix}
=
\begin{bmatrix}
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000
\end{bmatrix}}
其中， u_4 至 u_8 是网格黄色区域对紫色区域的边界， u_{9} 至 u_{12} 是网格第三行紫色区域的边界，分别相对两个计算域他们都是取上一步迭代出的结果，所以 u_4 至 u_{12} 都是已知量。
接下来将这些已知量移项至右侧：
{\tiny \begin{bmatrix}
400& -100& & & -100& & & & & & & & & & & \\
-100& 400& -100& & & -100& & & & & & & & & & \\
& -100& 400& -100& & & -100& & & & & & & & & \\
& & -100& 400& & & & -100& & & & & & & & \\
-100& & & & 400& -100& & & 0 & & & & & & & \\
& -100& & & -100& 400& -100& & & 0 & & & & & & \\
& & -100& & & -100& 400& -100& & & 0 & & & & & \\
& & & -100& & & -100& 400& & & & 0& & & & \\
& & & &0 & & & & 400& -100& & & -100& & & \\
& & & & & 0 & & & -100& 400& -100& & & -100& & \\
& & & & & & 0 & & & -100& 400& -100& & & -100& \\
& & & & & & & 0 & & & -100& 400& & & & -100\\
& & & & & & & & -100& & & & 400& -100& & \\
& & & & & & & & & -100& & & -100& 400& -100& \\
& & & & & & & & & & -100& & & -100& 400& -100\\
& & & & & & & & & & & -100& & & -100& 400
\end{bmatrix}
\begin{bmatrix}
u_1\\
u_2\\
u_3\\
u_4\\
{\color{Violet} u_5} \\
{\color{Violet} u_6}\\
{\color{Violet} u_7}\\
{\color{Violet} u_8}\\
{\color{Red} u_9}\\
{\color{Red} u_{10}}\\
{\color{Red} u_{11}}\\
{\color{Red} u_{12}}\\
u_{13}\\
u_{14}\\
u_{15}\\
u_{16}
\end{bmatrix}
=
\begin{bmatrix}
10000\\
0\\
0\\
50000\\
10000 {\color{red}+ 100 u_9}\\
0 {\color{red}+ 100 u_{10}}\\
0{\color{red}+ 100 u_{11}}\\
50000{\color{red}+ 100 u_{12}}\\
10000{\color{Violet}+ 100 u_{5}}\\
0{\color{Violet}+ 100 u_{6}}\\
0{\color{Violet}+ 100 u_{7}}\\
50000{\color{Violet}+ 100 u_{8}}\\
10000\\
0\\
0\\
50000
\end{bmatrix}}
{\tiny \begin{bmatrix}
400& -100& & & -100& & & & &&& & & & & & & \\
-100& 400& -100& & & -100& & & && & & & & & & & \\
& -100& 400& -100& & & -100& & && & & & & & & & \\
& & -100& 400& & & & -100& && & & & & & & & \\
-100& & & & 400& -100& & & &&0 & & & & & & & \\
& -100& & & -100& 400& -100& & && & 0 & & & & & & \\
& & -100& & & -100& 400& -100& && & & 0 & & & & & \\
& & & -100& & & -100& 400& &&& & & 0& & & & \\
\\
\\
& & & &0 & & & &&& 400& -100& & & -100& & & \\
& & & & & 0 & & & && -100& 400& -100& & & -100& & \\
& & & & & & 0 & & & && -100& 400& -100& & & -100& \\
& & & & & & & 0 & & & && -100& 400& & & & -100\\
& & & & & & & &&& -100& & & & 400& -100& & \\
& & & & & & & & && & -100& & & -100& 400& -100& \\
& & & & & & & & && & & -100& & & -100& 400& -100\\
& & & & & & & & && & & & -100& & & -100& 400
\end{bmatrix}
\begin{bmatrix}
u_1\\
u_2\\
u_3\\
u_4\\
{\color{Violet} u_5} \\
{\color{Violet} u_6}\\
{\color{Violet} u_7}\\
{\color{Violet} u_8}\\
\\
\\
{\color{Red} u_9}\\
{\color{Red} u_{10}}\\
{\color{Red} u_{11}}\\
{\color{Red} u_{12}}\\
u_{13}\\
u_{14}\\
u_{15}\\
u_{16}
\end{bmatrix}
=
\begin{bmatrix}
10000\\
0\\
0\\
50000\\
10000 {\color{red}+ 100 u_9}\\
0 {\color{red}+ 100 u_{10}}\\
0{\color{red}+ 100 u_{11}}\\
50000{\color{red}+ 100 u_{12}}\\
\\
\\
10000{\color{Violet}+ 100 u_{5}}\\
0{\color{Violet}+ 100 u_{6}}\\
0{\color{Violet}+ 100 u_{7}}\\
50000{\color{Violet}+ 100 u_{8}}\\
10000\\
0\\
0\\
50000
\end{bmatrix}}
容易看出，矩阵被划分为了两个部分，分别存在在两个进程中
进程0中的矩阵：
{\tiny \begin{bmatrix}
400& -100& & & -100& & & & \\
-100& 400& -100& & & -100& & \\
& -100& 400& -100& & & -100& \\
& & -100& 400& & & & -100& \\
-100& & & & 400& -100& & & \\
& -100& & & -100& 400& -100 \\
& & -100& & & -100& 400& -100 \\
& & & -100& & & -100& 400
\end{bmatrix}
\begin{bmatrix}
u_1\\
u_2\\
u_3\\
u_4\\
{ u_5} \\
{ u_6}\\
{ u_7}\\
{ u_8}\\
\end{bmatrix}
=
\begin{bmatrix}
10000\\
0\\
0\\
50000\\
10000 {\color{red}+ 100 u_9}\\
0 {\color{red}+ 100 u_{10}}\\
0{\color{red}+ 100 u_{11}}\\
50000{\color{red}+ 100 u_{12}}\\
\end{bmatrix}}
进程1中的矩阵：
{\tiny \begin{bmatrix}
400& -100& & & -100& & & \\
-100& 400& -100& & & -100& & \\
& -100& 400& -100& & & -100& \\
&& -100& 400& & & & -100\\
-100& & & & 400& -100& & \\
& -100& & & -100& 400& -100& \\
& & -100& & & -100& 400& -100\\
& & & -100& & & -100& 400
\end{bmatrix}
\begin{bmatrix}
{ u_9}\\
{ u_{10}}\\
{ u_{11}}\\
{ u_{12}}\\
u_{13}\\
u_{14}\\
u_{15}\\
u_{16}
\end{bmatrix}
=
\begin{bmatrix}
10000{\color{Violet}+ 100 u_{5}}\\
0{\color{Violet}+ 100 u_{6}}\\
0{\color{Violet}+ 100 u_{7}}\\
50000{\color{Violet}+ 100 u_{8}}\\
10000\\
0\\
0\\
50000
\end{bmatrix}}
两个进程的矩阵分在各自进程中进行Jacobi迭代（别的迭代方法也可以），可以发现，所有的格子都会获得新的值，包括作为内边界的 u_5 ～ u_{12} ，因此进程间需要通信，更新边界值，用计算出的新 u_5 ～ u_{12} 作为下一步迭代计算的值。
这是并行Jacobi方法，主要引自王承尧等的计算流体力学及其并行算法。
三、共轭梯度法的并行和我的疑问
在我学习并行Jacobi法求解矩阵之后，在思考共轭梯度法能否仿照并行Jacobi的求解。
但是最终结论是：我虽然成功在多进程下使用了CG法求解矩阵，但并不是按照上述Jacobi那样按网格分区分解矩阵，交换边界，再各自进程迭代的模式求解，碰到了一个还不太明白为什么的坑。
原始的CG法，也就是单进程下的CG法简述如下：
step 1: d_0 = r_0 = b - Ax_0
step 2: \alpha_i = \frac{r_i^Tr_i}{d_i^TAd_i}
step 3: x_{i+1} = x_i + \alpha_i A d_i
*step 4: r_{i+1} = r_i -\alpha_i A d_i
step 5: \beta_{i+1} = \frac{r_{i+1}^Tr_{i+1}}{r_{i}^Tr_{i}}
step 6: d_{i+1} = r_{i+1} + \beta_{i+1}d_i
step 7: 判断残差，不然返回 step2
那么仿照Jacobi在分区网格的并行，CG法在分区网格的并行就变成了：
step 1: d_0 = r_0 =b-Ax_0 （这里的 A 是该进程的 A ，也是如上边 Jacobi 并行迭代一样，经过了边界单元已知量的移项，将总矩阵拆分后的本地矩阵，后边会将各自进程分别进行的操作称为本地进程操作）
step 2: \alpha_i = \frac{r_i^Tr_i}{d_i^TAd_i} （这里分子上的残差内积实际上是本地进程里单元的残差内积）
step 3: x_{i+1} = x_i + \alpha_i A d_i
*step 4: r_{i+1} = b-Ax 这里计算新的残差计算需要引用到交换边界带来的新的 b（和原本方法的 step 4 不一样）
step 5: \beta_{i+1} = \frac{r_{i+1}^Tr_{i+1}}{r_{i}^Tr_{i}} （这里分子分母使用的残差内积都是本地单元的残差）
step 6: d_{i+1} = r_{i+1} + \beta_{i+1}d_i
step 7: 进程间交换边界
step 8: 各进程计算新的 b
step 9: 判断残差，不然返回 step 2
但是结果不尽如人意，单进程计算是成功的，但在两进程下如此生硬的套上“仿Jacobi法并行”的过程是无法求解的，计算的样子是在残差还很大的时候， \beta 就接近于 1 ，无法收敛了。
在经历的很久的琢磨之后，放弃了矩阵的分解形式的算法，回到矩阵移项之前的样子，也就是：
{\tiny \begin{bmatrix}
400& -100& & & -100& & & & & & & & & & & \\
-100& 400& -100& & & -100& & & & & & & & & & \\
& -100& 400& -100& & & -100& & & & & & & & & \\
& & -100& 400& & & & -100& & & & & & & & \\
-100& & & & 400& -100& & & {\color{Red} -100} & & & & & & & \\
& -100& & & -100& 400& -100& & & {\color{Red}-100}& & & & & & \\
& & -100& & & -100& 400& -100& & & {\color{Red}-100}& & & & & \\
& & & -100& & & -100& 400& & & & {\color{Red}-100}& & & & \\
& & & & {\color{Violet}-100}& & & & 400& -100& & & -100& & & \\
& & & & & {\color{Violet}-100}& & & -100& 400& -100& & & -100& & \\
& & & & & & {\color{Violet}-100}& & & -100& 400& -100& & & -100& \\
& & & & & & & {\color{Violet}-100}& & & -100& 400& & & & -100\\
& & & & & & & & -100& & & & 400& -100& & \\
& & & & & & & & & -100& & & -100& 400& -100& \\
& & & & & & & & & & -100& & & -100& 400& -100\\
& & & & & & & & & & & -100& & & -100& 400
\end{bmatrix}
\begin{bmatrix}
u_1\\
u_2\\
u_3\\
u_4\\
{\color{Violet} u_5} \\
{\color{Violet} u_6}\\
{\color{Violet} u_7}\\
{\color{Violet} u_8}\\
{\color{Red} u_9}\\
{\color{Red} u_{10}}\\
{\color{Red} u_{11}}\\
{\color{Red} u_{12}}\\
u_{13}\\
u_{14}\\
u_{15}\\
u_{16}
\end{bmatrix}
=
\begin{bmatrix}
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000\\
10000\\
0\\
0\\
50000
\end{bmatrix}}
求解步骤也就相应变成了：
step 1: d_0 = r_0 = b - Ax_0 （矩阵与向量的乘是向量，这里向量的每个值都储存在每个网格单元的数据结构里）
step 2: \alpha_i = \frac{r_i^Tr_i}{d_i^TAd_i} （这里残差的内积需要进行进程间通信进行累加，变成全局的残差内积，分母 d_i^TAd_i 同理）
step 3: x_{i+1} = x_i + \alpha_i A d_i
step 4: r_{i+1} = r_i -\alpha_i A d_i
step 5: \beta_{i+1} = \frac{r_{i+1}^Tr_{i+1}}{r_{i}^Tr_{i}} （同样的，这里的残差内积也是全局的残差内积）
step 6: d_{i+1} = r_{i+1} + \beta_{i+1}d_i
step 7: 判断残差，不然返回 step2
这样才能够成功求解，按这种方法，实际上并行的CG法每一步得到的每一个计算结果，和单进程的CG法全部都一样。
所以我的疑问是：
1、为什么并行CG法不能套用并行Jacobi的方法？（我自己的感觉是可能作为共轭梯度解的下降方向是需要全局的残差进行指向的）
2、经过这次挖坑和艰难埋坑，深感自己在并行求解矩阵方面知识的浅薄，希望能有大佬介绍几本并行矩阵求解的从基础到进阶的书籍，作为一个引路，先在这里感谢大佬了！
以上就是我的收获，总结和疑惑了，感谢大家的阅读！
更新 1
对问题 1 引用群里赞哥的解答：
- 不能这样做的原因需要知道一些Krylov subspace的数学理论知识 简单概括下就是 你这样把A进行分区并行的话 并行时候每一步迭代构建的Krylov subspace不是原来整个A的Krylov subspace 所以CG的收敛性的理论在你这都不成立 就不能保证收敛
对问题 2 ，也是赞哥推荐的参考资料
致谢
特别感谢赞哥和连续介质模拟群里各位老哥的鼎力相助。
