F. 小北买文具¶

分数：150 分

LU 分解¶

小北是个爱学习的学生，放学后经常去文具店买铅笔和钢笔。

有一天他翻出了之前的购物账单，发现上面的单价已经看不清了，只剩下购买数量和总价：

2 支铅笔和 1 支钢笔一共 7 元 1 支铅笔和 3 支钢笔一共 11 元

小北决定把单价算回来。设铅笔单价为 $x_1$ 元，钢笔单价为 $x_2$ 元，列方程：

$\begin{cases} 2x_1 + x_2 = 7 \\ x_1 + 3x_2 = 11 \end{cases}$

代入消元，很快就搞定了。

不过随着小北买的文具种类越来越多，而且他有个强迫症——买 $n$ 种文具就非要凑够 $n$ 张账单，问题就没那么简单了。用更一般的符号来写：

第 $i$ 张账单中，第 $j$ 种商品的数量记为 $a_{ij}$ （可以是任意实数，double 类型）
每种商品单价记为 $x_j$

第 $i$ 张账单对应：

$a_{i1}x_1 + a_{i2}x_2 + \cdots + a_{in}x_n = b_i$

写成完整的方程组就是：

$\begin{cases} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n = b_1 \\ a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n = b_2 \\ \vdots \\ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n \end{cases}$

把系数排成表格的话：

$\begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_n \end{pmatrix}$

这里出现了三个概念：

系数矩阵 $A$ ：所有方程的系数排成的表
未知量向量 $\mathbf{x}$ ：要求的那些单价
常数向量 $\mathbf{b}$ ：账单上写的总价

一长串方程压缩成一行：

$A\mathbf{x} = \mathbf{b}$

目标是把 $\mathbf{x}$ 算出来。

回想解二元方程时做的那些事——交换方程顺序、用一个方程的若干倍去消掉另一个方程里的未知数——本质上都没有改变解，只是换了个写法。消元允许的操作就三种：

交换两行
用非零常数乘某一行
把一行的倍数加到另一行上

反复做这些行变换，主对角线下方的元素会被逐步消掉，矩阵变成上三角形式，这时候回代就能求解了。这就是高斯消元。

消元过程中，为了消掉第 $i$ 列下面的元素，我们会做：

$\text{第 } j \text{ 行} \;\leftarrow\; \text{第 } j \text{ 行} - l_{ji} \times \text{第 } i \text{ 行}$

这里的 $l_{ji}$ 就是消元时用到的倍数。把这些倍数收集起来，它们本身构成一个单位下三角矩阵 $L$ ；消元结束后剩下的上三角矩阵记为 $U$ 。于是有：

$A = LU$

$U$ ：消元的结果
$L$ ：消元过程中用过的倍数

但实际计算中会遇到一个问题：主元位置上的值可能是 0 或者非常小，直接拿来做除法会引入很大的数值误差。解决办法是交换行，在当前列里选一个绝对值更大的元素做主元。（如果从主元往下整列都接近 0，那说明矩阵本身就是奇异的，不是消元策略能解决的。本题保证不会出现这种情况。）

把行交换也纳入进来，就需要一个置换矩阵 $P$ ：

$PA = LU$

这就是带部分选主元的 LU 分解（PLU 分解）。

总结一下：

高斯消元：一步步把矩阵化简
LU 分解：把消元过程记录成 $L$ 和 $U$
置换矩阵 $P$ ：记录行交换

实现上，LU 分解可以原地完成，不需要额外分配整个矩阵的空间。下发的 solver_naive.c里提供了一个朴素实现，其中mydgetrf 接受 $n \times n$ 矩阵 $A$ ：

对 $A$ 做 LU 分解，结果直接存回 $A$ （下三角存 $L$ ，上三角存 $U$ ）
行交换通过长度为 $n$ 的整数数组 ipiv 记录，相当于隐式存储了 $P$

有了 LU 分解，解 $Ax = b$ 的步骤是：

改写为 $PAx = Pb$
代入得 $LUx = Pb$
先解下三角方程 $Ly = Pb$
再解上三角方程 $Ux = y$

两次回代就搞定了，而且解直接覆盖到向量里，不用额外开空间。

LU 分解和朴素高斯消元的复杂度都是 $O(n^3)$ ，那 LU 分解有什么好处？考虑这样一个场景： $A$ 不变， $b$ 换了好多次。如果每次都从头消元，那每次都是 $O(n^3)$ ；但如果先分解一次 $A = LU$ ，以后每次只要做两次回代，复杂度就降到了 $O(n^2)$ 。

不过朴素实现的性能不太行。 $n=100$ 的时候瞬间出结果， $n=32768$ 的时候就明显慢了——CPU 占用率拉满，但 GFLOPS 上不去，说明计算效率很低。

另外还有一个问题：现在的服务器都是多核的，朴素的单线程实现只用了其中一个核，剩下的全浪费了。

提示¶

OpenMP 是常用的共享内存并行编程工具，加几行 #pragma omp parallel for 就能把循环分给多个线程跑。但光加 OpenMP 不一定快——如果内存访问模式不好，内存带宽会成为瓶颈。

分块（Blocking/Tiling） 是改善局部性的关键手段。把大矩阵切成小块，每次把一小块加载到 Cache 里算完再换下一块，尽量让 CPU 在 Cache 里干活而不是反复去主存取数据。

两者结合起来：分块把数据留在 Cache 里，OpenMP 让多个核心并行算 Cache 里的数据，CPU 的算力才能真正发挥出来。

任务描述¶

你的任务是编写 solver.c，实现基于分块策略的 LU 分解算法，并用 OpenMP 做多线程并行加速。

具体要求：

分块算法：把大矩阵拆成小块，利用 Cache 局部性提高效率。
并行化：用 OpenMP 充分利用多核。
尽量写对编译器友好的代码。

1. 接口规范¶

数据交互用二进制文件。你不需要写 main 函数，只需要实现核心求解逻辑。

Driver 程序：我们提供了 driver.c，负责读数据、计时、调你的函数、写结果。不要改 driver.c。
Solver 接口：你要在 solver.c 里实现：

// n: 矩阵大小（最大 32768
// A: 系数矩阵 (n*n)，行优先存储
// b: 常数向量 (n)
// 任务：解 Ax=b，把解存到 b 里。
// 注意：可以（也应该）直接在 A 上做 LU 分解，破坏原内容没关系。
void my_solver(int n, double *A, double *b);

2. 编译与运行环境¶

编译命令：¶

gcc -O3 -march=native -fopenmp driver.c solver.c -o solver -lm

硬件：纯 CPU 题，评测机是 Intel(R) Xeon(R) Platinum 8358 CPU @ 2.60GHz 。
线程数：评测时 OMP_NUM_THREADS 设为 32（物理核心数），并绑定核心。

运行命令：¶

export OMP_NUM_THREADS=32
export OMP_PLACES=cores
export OMP_PROC_BIND=close
./solver input.bin output.bin

容器环境¶

镜像：crmirror.lcpu.dev/hpcgame/base:latest

3. 数据规模与评分¶

矩阵维度 $N$ 最大32768，双精度浮点。
正确性：计算相对残差 $r = \frac{||Ax - b||_2}{||b||_2}$
合法性：
用下发的 Makefile 测试，不能用外部数学库（OpenBLAS 、 MKL 、 LAPACK 、 Eigen 等的线性代数函数如 dgetrf、dgemm 都不行）。自己写。
允许用 C/C++ 标准库（math.h、stdlib.h 等）和 OpenMP 。
性能评分：用 bash run_judge.sh 本地测试，最终以集群评测为准。

分数计算方式：

程序运行时间为 t，满分时间为t_0，满分为u，则得分为min(u, u^(t_0/t)) 。

测试点	矩阵规模 $N$	满分成绩	时间限制	分值
1	2049	0.1s	0.5s	10
2	4096	0.8s	3s	20
3	8192	5.3s	10s	20
4	16384	36.5s	75s	20
5	32768	268s	500s	30

Wish you good luck!

附件¶

handout-lu