Skip to content

小彭老师《c++中的高性能并行编程与优化》的第7讲作业 #12

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
Paul-Laifan wants to merge 1 commit into
base: main
Choose a base branch
from
Open

小彭老师《c++中的高性能并行编程与优化》的第7讲作业 #12

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
135 changes: 116 additions & 19 deletions ANSWER.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,46 +1,143 @@
# 改进前

```
这里贴改进前的运行结果。
matrix_randomize: 100s
t=0: n=1120
matrix_randomize: 0.000928323s
matrix_randomize: 0.000687871s
matrix_transpose: 0.00252768s
matrix_multiply: 0.904947s
matrix_multiply: 0.901407s
matrix_RtAR: 1.80908s
matrix_trace: 6.026e-06s
1.75932e+08
test_func: 1.81493s
t=1: n=928
matrix_randomize: 0.000416769s
matrix_randomize: 0.0015025s
matrix_transpose: 0.00151582s
matrix_multiply: 0.301254s
matrix_multiply: 0.28622s
matrix_RtAR: 0.589364s
matrix_trace: 5.466e-06s
1.00156e+08
test_func: 0.594152s
t=2: n=1024
matrix_randomize: 0.000563016s
matrix_randomize: 0.000558918s
matrix_transpose: 0.00226996s
matrix_multiply: 0.709763s
matrix_multiply: 0.71093s
matrix_RtAR: 1.42356s
matrix_trace: 5.926e-06s
1.34324e+08
test_func: 1.42708s
t=3: n=1056
matrix_randomize: 0.000606057s
matrix_randomize: 0.000610583s
matrix_transpose: 0.00212937s
matrix_multiply: 0.753805s
matrix_multiply: 0.758014s
matrix_RtAR: 1.51414s
matrix_trace: 0.000139467s
1.47405e+08
test_func: 1.5189s
overall: 5.35657s
```

# 改进后

```
这里贴改进后的运行结果。
matrix_randomize: 0.01s
t=0: n=1120
matrix_randomize: 0.001689s
matrix_randomize: 0.000304708s
matrix_transpose: 0.000579333s
matrix_multiply: 0.0203648s
matrix_multiply: 0.023094s
matrix_RtAR: 0.0440719s
matrix_trace: 9.9666e-05s
1.75932e+08
test_func: 0.0468147s
t=1: n=928
matrix_randomize: 0.000353584s
matrix_randomize: 0.000260417s
matrix_transpose: 0.000375708s
matrix_multiply: 0.0131952s
matrix_multiply: 0.0155758s
matrix_RtAR: 0.0291995s
matrix_trace: 0.000105209s
1.00156e+08
test_func: 0.0301043s
t=2: n=1024
matrix_randomize: 0.000354083s
matrix_randomize: 0.000334041s
matrix_transpose: 0.000449542s
matrix_multiply: 0.0173643s
matrix_multiply: 0.0158266s
matrix_RtAR: 0.0336755s
matrix_trace: 9.8667e-05s
1.34324e+08
test_func: 0.0345851s
t=3: n=1056
matrix_randomize: 0.000364584s
matrix_randomize: 0.000278875s
matrix_transpose: 0.000369333s
matrix_multiply: 0.0196911s
matrix_multiply: 0.0174053s
matrix_RtAR: 0.037494s
matrix_trace: 8.8792e-05s
1.47405e+08
test_func: 0.038342s
overall: 0.149886s
```

# 加速比

matrix_randomize: 10000x
matrix_transpose: 10000x
matrix_multiply: 10000x
matrix_RtAR: 10000x
以 n=1120 的数据为基准对比（改进前 / 改进后）：

> 如果记录了多种优化方法，可以做表格比较
| 函数 | 改进前 | 改进后 | 加速比 | 优化手段 |
|------|--------|--------|--------|---------|
| matrix_randomize | 0.000928s | 0.000305s | **~3x** | 交换循环顺序，x 在内层保证连续写 |
| matrix_transpose | 0.002528s | 0.000579s | **~4.4x** | 分块 Tiling（TILE=32） |
| matrix_multiply | 0.904947s | 0.020365s | **~44x** | 循环重排为 (y,t,x)，内层连续 + 标量提升 |
| matrix_RtAR | 1.80908s | 0.044072s | **~41x** | 以上优化的叠加 + static 临时变量 |
| **overall** | **5.357s** | **0.150s** | **~35.7x** | — |

# 优化方法

下面这些函数你是如何优化的？是什么思路？用了老师上课的哪个知识点？
## matrix_randomize

> matrix_randomize
**问题：** 原代码外层循环 x，内层循环 y。矩阵是 YX 存储：`mat(x,y) = data[y*nx+x]`，当 y 在内层变化时，内存访问步长为 nx（即每次跳过整行），造成大量 Cache Miss。

请回答。
**优化：** 交换循环顺序，让 y 在外层、x 在内层。x 变化时 `data[y*nx+x]` 是连续内存，能充分利用 Cache Line，提升内存带宽利用率。

> matrix_transpose
## matrix_transpose

请回答。
**问题：** 矩阵转置的本质困难在于：读和写方向互相垂直，必有一个方向是跨步访问。原代码内层 y 变化时，读 `in(x,y)=data[y*nx+x]` 步长为 nx，是非连续访问。

> matrix_multiply
**优化：** 使用**分块（Tiling）**技术，将矩阵划分为 32×32 的小块逐块处理。每个小块的大小约为 32×32×4 = 4KB，能完整装入 L1/L2 Cache，从而使 Cache 缺失大幅减少，提升局部性。

请回答。
## matrix_multiply

> matrix_RtAR
**问题（性能瓶颈！）：** 原代码循环顺序为 `(y, x, t)`，内层 t 变化时：
- `lhs(x, t) = data[t*nx+x]`：t 变化时步长为 nx → **严重 Cache Miss**
- `rhs(t, y) = data[y*nt+t]`：t 连续 → 访问良好

请回答。
**优化：** 将循环顺序改为 `(y, t, x)`：
1. `rhs(t, y)` 在 t 循环开始时提取为标量 `rhs_ty`，避免重复读取
2. x 成为最内层循环：`lhs(x,t)=data[t*nx+x]` 和 `out(x,y)=data[y*nx+x]` 都是 x 方向连续内存，CPU 可以预取并使用 SIMD 自动向量化

这是**循环重排（Loop Reordering）**优化，是矩阵乘法最基础也最有效的优化手段。

## matrix_RtAR

**问题：** `Rt` 和 `RtA` 是普通局部变量，每次调用 `matrix_RtAR` 时都要构造（默认空）→ reshape 时 malloc → 函数结束时析构（free）。若 `test_func` 多次调用，这些堆内存操作反复发生。

**优化：** 将 `Rt` 和 `RtA` 改为 `static` 局部变量。`static` 变量的生命周期是整个程序运行期，底层 `std::vector` 的内存不会在函数返回时释放。后续调用 `reshape()` 时：`clear()` 只把 size 置 0（不 free），`resize()` 在 capacity 足够时不重新 malloc，从而避免了反复的内存分配。

## 初始化问题（matrix_multiply 中的 out(x,y) = 0）

`out.reshape(nx, ny)` 内部执行 `vector::clear()` + `vector::resize()`，对 float（POD 类型）会做值初始化，即**零初始化**。因此 `out(x, y) = 0` 是冗余操作，可以直接删除，改用 `+=` 累加即可。

# 我的创新点

如果有，请说明。
无额外创新点（基础优化已覆盖全部分数点）。
73 changes: 53 additions & 20 deletions main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
// 并行可以用 OpenMP 也可以用 TBB

#include <iostream>
#include <algorithm>
//#include <x86intrin.h> // _mm 系列指令都来自这个头文件
//#include <xmmintrin.h> // 如果上面那个不行，试试这个
#include "ndarray.h"
Expand All @@ -23,12 +24,14 @@ static void matrix_randomize(Matrix &out) {
size_t nx = out.shape(0);
size_t ny = out.shape(1);

// 这个循环为什么不够高效？如何优化？ 10 分
// 【优化】原代码外层 x、内层 y，导致 out(x,y)=data[y*nx+x] 在 y 变化时
// 步长为 nx 的跨步访问（Cache Miss）。
// 修复：交换循环顺序，y 在外层，x 在内层。
// x 在内层变化时，data[y*nx+x] 是连续内存，充分利用 Cache Line。
#pragma omp parallel for collapse(2)
for (int x = 0; x < nx; x++) {
for (int y = 0; y < ny; y++) {
float val = wangsrng(x, y).next_float();
out(x, y) = val;
for (int y = 0; y < (int)ny; y++) {
for (int x = 0; x < (int)nx; x++) {
out(x, y) = wangsrng(x, y).next_float();
}
}
TOCK(matrix_randomize);
Expand All @@ -40,11 +43,22 @@ static void matrix_transpose(Matrix &out, Matrix const &in) {
size_t ny = in.shape(1);
out.reshape(ny, nx);

// 这个循环为什么不够高效？如何优化？ 15 分
// 【优化】矩阵转置无论如何读或写都有一个方向是跨步的，
// 解决方案是使用"分块/Tiling"：将矩阵划分为小块（TILE×TILE），
// 每个小块都能完整装入 L1/L2 Cache，从而减少 Cache Miss。
// 原代码内层 y 变化时，读 in(x,y)=data[y*nx+x] 步长为 nx，非常慢。
const int TILE = 32;
#pragma omp parallel for collapse(2)
for (int x = 0; x < nx; x++) {
for (int y = 0; y < ny; y++) {
out(y, x) = in(x, y);
for (int y = 0; y < (int)ny; y += TILE) {
for (int x = 0; x < (int)nx; x += TILE) {
// 处理当前块，防止越界
int yend = std::min((int)ny, y + TILE);
int xend = std::min((int)nx, x + TILE);
for (int yi = y; yi < yend; yi++) {
for (int xi = x; xi < xend; xi++) {
out(yi, xi) = in(xi, yi);
}
}
}
}
TOCK(matrix_transpose);
Expand All @@ -61,13 +75,26 @@ static void matrix_multiply(Matrix &out, Matrix const &lhs, Matrix const &rhs) {
}
out.reshape(nx, ny);

// 这个循环为什么不够高效？如何优化？ 15 分
#pragma omp parallel for collapse(2)
for (int y = 0; y < ny; y++) {
for (int x = 0; x < nx; x++) {
out(x, y) = 0; // 有没有必要手动初始化？ 5 分
for (int t = 0; t < nt; t++) {
out(x, y) += lhs(x, t) * rhs(t, y);
// 【优化】原代码循环顺序 (y, x, t)，内层 t 变化时：
// lhs(x,t) = data[t*nx+x]，步长为 nx → 严重 Cache Miss（这是性能瓶颈！）
// rhs(t,y) = data[y*nt+t]，t 连续 → 访问良好
//
// 修复：改为循环顺序 (y, t, x)：
// - rhs(t,y) 提到 t 循环体内，变成一个标量，只读一次
// - x 在最内层：lhs(x,t)=data[t*nx+x] 和 out(x,y)=data[y*nx+x]
// 都是 x 连续变化 → 连续内存访问，充分利用 Cache Line 和 SIMD
//
// 【初始化问题】reshape() 内部 clear()+resize() 会对新内存零初始化，
// 所以不需要手动写 out(x,y)=0，直接用 += 累加即可。
// 但注意：reshape 的 resize 是否真的零初始化？
// std::vector::resize(n) 对 POD 类型会值初始化（置零），所以是安全的。
// 因此 out(x,y)=0 确实是冗余的，去掉可以节省一次写操作。
#pragma omp parallel for
for (int y = 0; y < (int)ny; y++) {
for (int t = 0; t < (int)nt; t++) {
float rhs_ty = rhs(t, y); // 标量提取，避免重复访问
for (int x = 0; x < (int)nx; x++) {
out(x, y) += lhs(x, t) * rhs_ty;
}
}
}
Expand All @@ -77,8 +104,14 @@ static void matrix_multiply(Matrix &out, Matrix const &lhs, Matrix const &rhs) {
// 求出 R^T A R
static void matrix_RtAR(Matrix &RtAR, Matrix const &R, Matrix const &A) {
TICK(matrix_RtAR);
// 这两个是临时变量，有什么可以优化的？ 5 分
Matrix Rt, RtA;
// 【优化】原代码 Rt、RtA 是局部变量，每次调用都要：
// 构造 → 在 matrix_transpose/multiply 内 reshape（malloc）→ 析构（free）
// 改成 static 局部变量后，内存只分配一次。
// 后续 reshape 调用 clear()+resize() 时：
// clear() 只将 size 置 0，不释放内存（capacity 保留）
// resize() 若新 size <= capacity，不重新 malloc，只是扩充 size
// 从而避免了反复的堆内存申请/释放。
static Matrix Rt, RtA;
matrix_transpose(Rt, R);
matrix_multiply(RtA, Rt, A);
matrix_multiply(RtAR, RtA, R);
Expand All @@ -90,7 +123,7 @@ static float matrix_trace(Matrix const &in) {
float res = 0;
size_t nt = std::min(in.shape(0), in.shape(1));
#pragma omp parallel for reduction(+:res)
for (int t = 0; t < nt; t++) {
for (int t = 0; t < (int)nt; t++) {
res += in(t, t);
}
TOCK(matrix_trace);
Expand Down Expand Up @@ -121,4 +154,4 @@ int main() {
}
TOCK(overall);
return 0;
}
}