GPU Gems 2 · Chapter 44. A GPU Framework for Solving Systems of Linear Equations

把偏微分方程离散成网格之后，问题不再是“画出一张纹理”，而是不断对大量数值做加法、乘法、点积和求解。CPU 可以直接使用数组和指针；GPU 则更擅长纹理取样和并行写入，因此需要先约定数值对象在纹理里的形状。

本章的框架把底层存储和 fragment program 封装成向量、矩阵、归约等可组合的操作。上层算法只关心 A * x、x + αp 以及 residual 是否变小；底层负责把它们变成纹理分配、坐标、渲染目标和临时内存。

1. 统一表示：标量和向量也要有纹理地址

GPU linear algebra framework 的第一步不是写 solver，而是设计表示。一个 GPU 上的单浮点值看似可以直接作为 uniform，但许多操作会生成并再次读取标量；把它放在单元素纹理里，才能和向量、矩阵共享同一套访问路径。

一维向量可以重排成二维纹理，让相邻数据落在相近的 texel。为了减少内部表示大小，原章进一步把连续的 2×2 元素打包到一个 RGBA texel。这个布局不是视觉上的“把四个颜色混在一起”，而是一个可逆的地址协议：给定逻辑索引，就能算出 texel 坐标和 channel。

texture representation 还决定边界行为：长度不是 2 的倍数时要处理 padding，矩阵对角线或稀疏条带要保持索引一致，临时纹理要在下一次操作前可复用。表示层一旦稳定，上面的数值算子才能组合。

2. 按非零结构选择矩阵布局

矩阵经常比向量稀疏。dense 矩阵几乎每个元素都非零；banded 矩阵只在若干条对角带上有值；sparse 矩阵的非零项分布更零散。三者不能只换一个类型名，因为它们的 texture fetch、地址计算、内存占用和输出写法都不同。

对由 PDE 离散得到的邻域算子，banded 布局通常比完整 dense 矩阵更合适；对通用输入或小矩阵，dense 可能更简单。原章的框架支持多种矩阵和向量类型，目标是让上层 solver 不必知道每个非零值究竟存在哪一张纹理里。

3. matrix-vector：让位置和坐标承担索引

matrix-vector product 是共轭梯度和波方程更新的核心。对 x = A · b，输出 xᵢ 等于第 i 行与 b 的点积。GPU 实现可以让点的位置编码矩阵列索引，让纹理坐标编码向量取样位置；渲染这些点时，结果自然落到输出向量的正确位置。

// 抽象 API：隐藏矩阵的实际纹理布局
clMatrix->matrixVectorOp(CL_NULL, p, nullptr, q);
// q = A * p
clR->subtractVector(q, r, 1.0f);
// r = r - q（具体系数由 solver 控制）

这里的 API 名称体现框架目标，而不是要求读者照抄旧版 C++。关键是每个算子有明确的输入、输出和布局契约：不能因为换成 banded 矩阵就改变上层 matrixVectorOp 的语义。

4. reduction 和可复用的临时内存

共轭梯度需要 dot(r, r)、dot(p, q) 这类 reduction。逐元素乘法仍然适合 fragment unit，但把 n 个结果收成一个标量需要反复缩小纹理：相邻元素先合并，再对中间结果继续合并，直到只剩一个值。

框架因此需要一个虚拟内存管理层，为 reduction 和中间向量共享临时纹理。复用的好处不是只省一次分配，而是让多轮 solve 不必不断创建和销毁 render target；代价是生命周期必须清楚，任何仍被后续算子读取的纹理都不能提前归还。

5. 共轭梯度：把基础算子组成稳定迭代

conjugate gradient 从初始解 x 开始计算残差 r = b − A x，令搜索方向 p = r。每轮先算 q = A p，用两个 reduction 得到步长 α 和方向系数 β，然后更新 x、r、p。每个变量都可以对应框架里的向量对象，矩阵只通过 matrix-vector 算子参与。

void solveIteration() {
    A->matrixVectorOp(CL_NULL, P, nullptr, Q); // Q = A * P
    P->reduceAdd(Q, Temp);                     // temp = dot(P, Q)
    Rho->div(Temp, Alpha);                    // alpha = rho / temp
    X->addVector(P, X, 1.0f, Alpha);          // X += alpha * P
    R->subtractVector(Q, R, 1.0f, Alpha);     // R -= alpha * Q
    R->reduceAdd(R, NewRho);                  // newrho = dot(R, R)
}

循环的停止条件可以是 residual 小于阈值，也可以是达到最大迭代次数。它并不保证任意矩阵都快速收敛：矩阵条件数、初始猜测、浮点精度和 reduction 误差都会影响实际表现。GPU 框架的价值在于把这些变量放到同一个可测量的算子图里。

先预测：把 banded 切换为 dense，残差曲线会改变，还是只有 storage count 改变？把 RGBA packing 打开时，solver 的数学结果应该保持不变；如果结果变了，优先检查地址和 channel 映射，而不是怀疑共轭梯度公式。

6. 从线性系统到 2D 波方程

Crank-Nicolson scheme 用当前和下一时刻的空间差分平均来离散 2D 波方程。显式方案可以直接从 current 和 last 算出 next，代码短且每次只做一轮矩阵-向量操作；但时间步过大就会不稳定。

隐式方案把下一帧未知高度收进 A · x(next) = b(current)，因此每个时间步都要调用一次线性系统 solver。这里共轭梯度不再是抽象的 benchmark，而是让水面或流体模拟在更宽松的时间步下保持稳定的工具。代价是每帧需要若干轮 matrix-vector、reduction 和向量更新。

一个可审计的实现顺序是：先用 CPU 验证 central difference 的离散式，再用 GPU framework 验证单次 matrix-vector，然后比较显式与隐式同一时间步的结果，最后才优化 texture packing 和临时内存。这样可以把“数学错误”“布局错误”和“稳定性错误”分开定位。

小结

• GPU 线性代数框架先统一标量、向量和矩阵的纹理表示。
• dense、banded、sparse 的非零结构决定不同的存储与 matrix-vector 路径。
• reduction 是点积和 residual 的基础，需要分层归约与临时纹理管理。
• 共轭梯度由 matrix-vector、归约和向量更新组合而成。
• Crank-Nicolson 通过求解线性系统换取波方程的条件稳定性改善。