前置讲 P1 · 从 Python 到 C：指针、数组与内存布局 · From Python to C: Pointers, Arrays & Memory Layout

这是一篇前置讲，不对应原课的某一次课。原课的先修要求写着”会 Python 以及 C / C++ / Java 之一”，第 2 讲起就开始写 CUDA 内核——而 CUDA 是 C/C++ 的扩展。原课第 3 讲 slides 上的第一段代码是 cudaMalloc(&dA, n * sizeof(int))，第 4 讲的第一个内核里有 a[row * N + k]。这篇为只会 Python 的读者，把读懂、写出这些代码所需的 C 概念从头讲清楚。 依据说明：本讲讲的是通用的 C 语言与计算机知识（类型、指针、栈与堆、行优先存储），参考 Kernighan & Ritchie《C 程序设计语言》、NVIDIA CUDA C++ Programming Guide 的内存管理接口、PMPP（Programming Massively Parallel Processors，第 4 版）第 2 章的向量加法例子、NumPy 文档里关于 strides 的说明，以及原课 HW1 / miniTorch 的张量数据结构（index_to_position、to_index、broadcast_index）。所有实验都是本讲义自己设计的。 配套代码 Code：code/prereq1/（3 个脚本 + 2 个模块 + 2 个 C 源文件；需要 numpy 和一个 C 编译器；每个脚本 1–2 秒） 预计阅读时间：3 小时（含运行代码） 建议阅读时机：第 2 讲之前（写 CUDA 内核必需）

0. 本讲导航 · How to read this lecture

0.1 为什么要懂这些

下面是一个完整的 CUDA 程序的骨架（向量加法，原课第 2–3 讲和 PMPP 第 2 章的标准例子；本讲义的环境没有 NVIDIA GPU，这段代码只用来”读”）：

__global__ void vecAdd(const float* A, const float* B, float* C, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;   // 这个线程负责第 i 个元素
    if (i < n) C[i] = A[i] + B[i];                   // 越界检查 + 计算
}

int main(void) {
    int n = 1 << 20;
    size_t size = n * sizeof(float);
    float* A_h = (float*)malloc(size);               // CPU 内存（host）
    /* …… 给 A_h、B_h 填数据 …… */
    float* A_d;                                      // 将来指向 GPU 显存（device）
    cudaMalloc((void**)&A_d, size);                  // 为什么是 &A_d？为什么要转成 void**？
    cudaMemcpy(A_d, A_h, size, cudaMemcpyHostToDevice);
    vecAdd<<<(n + 255) / 256, 256>>>(A_d, B_d, C_d, n);
    cudaMemcpy(C_h, C_d, size, cudaMemcpyDeviceToHost);
    cudaFree(A_d);  free(A_h);
    return 0;
}

如果你只会 Python，这里几乎每一行都有问题：float* 是什么？sizeof(float) 是多少？为什么 cudaMalloc 要传 &A_d、还要转成 (void**)？A[i] 越界了会怎样？malloc 和 free 为什么要成对出现？<<<…>>> 之外的部分其实都是普通的 C。这篇前置讲把这些概念讲清楚，§7 会回到这段代码，逐行解释。

0.2 学习目标 Learning goals

1. 说清楚编译型语言和解释型语言的区别；知道一个 C 程序从源代码到可执行文件经过哪几步；能用实验解释 Python 循环为什么慢，以及”调用一次 C 函数”本身的开销。
2. 知道 C 的基本类型各占几个字节，一个数在内存里是哪几个字节（字节序）。
3. 理解指针：取地址 &、取值 *、指针运算、数组与指针的关系、按指针传参、指向指针的指针（void**）。
4. 理解栈和堆：局部变量的生命周期、malloc / free、内存泄漏和悬空指针；知道 GPU 显存的 cudaMalloc / cudaFree / cudaMemcpy 是同一套思路。
5. 理解多维数组的行优先存储和步长：能算出任意元素的一维位置，知道转置、切片、广播为什么是”零拷贝的视图”、什么时候必须复制；能写出原课 HW1 里的 index_to_position、to_index、broadcast_index。
6. 会用 ctypes 从 Python 调用 C 函数、传递 NumPy 数组和结构体，并知道为什么”一次处理一大批数据”比”每个元素调用一次”快得多。

0.3 本讲地图 Map

§1  编译型 vs 解释型：C 程序的一生；Python 为什么慢；求和实验；跨语言调用的开销      ← 代码 01
§2  类型与字节：sizeof；字节序；读懂 CUDA 所需的最少 C 语法                         ← 代码 02 第 1–2 部分
§3  指针：& 与 *；指针运算；越界；按指针传参；void**                                ← 代码 02 第 3–5 部分
§4  内存的生命周期：栈与堆；malloc / free；泄漏；GPU 上的对应                        ← 代码 02 第 6 部分
§5  多维数组：行优先、步长、转置、切片、广播；HW1 式下标函数                          ← 代码 03 第 1–4 部分
§6  结构体与 Python ↔ C：Tensor2D、ctypes、共享库；为什么还是要用库                  ← 代码 03 第 5–6 部分
§7  回到开头：逐行读懂一个 CUDA 程序
§8  小结   §9  常见误区   §10 练习   §11 延伸阅读   §12 下一篇预告

0.4 先修

会写 Python 函数、循环、列表；知道 NumPy 数组是什么。不需要任何 C 的经验。代码需要一个 C 编译器：macOS 上运行一次 xcode-select --install，Linux 上安装 gcc（sudo apt install build-essential）。

1. 编译型 vs 解释型 · Compiled vs interpreted

1.1 一个 C 程序的一生

直觉。 Python 程序像是请了一位同声传译：你每说一句，他当场翻译一句，下次再说同样的话，他还得再翻译一遍。C 程序像是一本事先译好的书：翻译只做一次，之后读者直接读译文。

C 的”翻译”由编译器完成，分四步。01 号脚本第 1 部分把 hello_pointers.c（一个 58 行的小程序）按步骤拆开编译，统计了每一步的产物（本机 Apple clang 21）：

平时一条 cc -O2 -o hello hello_pointers.c 就把四步一起做了。-O2 是优化级别：让编译器花更多时间，生成更快的机器码。

这些工作都在运行之前一次做完。运行时，CPU 直接执行机器指令，不需要再”理解”源代码。CUDA 也一样：.cu 文件先用 NVIDIA 的编译器 nvcc 编译，CPU 部分变成普通机器码，GPU 部分（内核）变成 GPU 的机器码；运行时 CPU 把内核”发射”到 GPU 上执行。

1.2 Python 为什么慢

Python 是解释型的：源代码先被翻译成字节码（一种虚拟机指令），解释器再一条一条地执行字节码，每执行一条都要检查操作数的类型、决定调用哪段 C 代码。01 号脚本第 2 部分用 dis 模块看了一个求和循环：

def add_all(xs):
    s = 0.0
    for x in xs:
        s += x
    return s

每循环一次要执行 6 条字节码（Python 3.13）：FOR_ITER → STORE_FAST → LOAD_FAST_LOAD_FAST → BINARY_OP → STORE_FAST → JUMP_BACKWARD。BINARY_OP 要先看 s 和 x 是什么类型（float？int？自定义类？），再调用对应的加法，结果还要装进一个新的 float 对象。

数据的存放方式也不一样：

Python 的灵活（变量可以随时换类型、整数可以任意大）是用这些开销换来的。

1.3 实验：把一百万个数加起来

01 号脚本第 3 部分用五种写法把 100 万个 [0, 1) 之间的随机 float32 加起来（⏱ 本机实测，取 5 次中最快的一次）：

（“精确值”用 float64 求和得到：500377.2696。）

这张表里有三个值得细想的地方，其中两个和我写代码前的预期不一样：

(1) C 只比 Python 快 14 倍左右，开 -O2 也只多快 1.2 倍。 我原以为编译器优化会带来明显的加速。可是这个循环有一个依赖链：第 i 次加法必须等第 i−1 次加法的结果出来才能开始。CPU 做一次浮点加法要几个时钟周期，于是每个元素至少要等这么久。要打破依赖链，就得把数组拆成几段、各自累加、最后再合起来——但这会改变加法的顺序，而浮点加法不满足结合律（P2），结果会变。C 编译器默认不允许自己改变结果，所以不会这样做。

(2) NumPy 比我们的 C 循环还快约 5 倍。 NumPy 的求和恰好做了编译器不敢做的事：用 SIMD 指令（一条指令同时处理好几个数）维护多个部分和，最后两两相加（成对求和，pairwise summation）。这就是并行归约（parallel reduction）：GPU 上的 reduce 内核（原课 HW1 的 Problem 3、第 3 讲）是同一个思路，只是并行度高得多——几千个线程各自加一段，再一层层合并。

(3) 结果不一样。 用 32 位 float 顺序累加，差了 3.83：当部分和长到几十万时，再加一个不到 1 的数，舍入误差就不可忽略了。用 double 累加没有可见误差；NumPy 的成对求和把每一步的部分和都控制在相近的大小，误差只有 0.02。换一种加法顺序，结果就在最后几位上不同——这是 GPU 程序的输出常常不能逐位复现的根源，也是 GPU 上做归约时要用高精度累加器的原因。P2 会专门讲。

1.4 跨语言调用的开销

如果 C 这么快，能不能在 Python 循环里对每个元素调用一次 C 函数？01 号脚本第 4 部分测了一下（⏱）：

跨语言调用一次约 143 纳秒，相当于 C 里做约 270 次加法。“每个元素调用一次 C 函数”会把 C 的速度全部吃掉；正确的做法是一次调用处理一大批数据。这条原则在 GPU 上更重要：CPU 发射一个 GPU 内核要几微秒，如果内核只处理很少的数据，时间几乎全花在发射上。第 10–11 讲的”算子融合”要把多个小内核合成一个，原因之一就在这里。

1.5 Python + C：今天的标准分工

PyTorch、NumPy 的核心都是 C/C++（和 CUDA）写的，Python 只是”指挥官”：它决定调用哪个算子、按什么顺序，真正的计算全在编译好的代码里。原课的 miniTorch 也是这样：张量和自动微分用 Python 写，CUDA 内核编译成共享库，再从 Python 调用（HW1）。本讲的 c_lab.py 用 ctypes 做的是同一件事的最小版本（§6.2）。

2. 类型与字节 · Types and bytes

2.1 C 的基本类型

C 的每个变量都要先声明类型，类型决定它占几个字节、怎样解释这些字节。02 号脚本第 1 部分让 C 自己报告了 sizeof（本机 64 位 macOS）：

几点说明：

• long 的大小因平台而异。需要确定位数时用 <stdint.h> 里的 int32_t、int64_t、uint8_t 等，它们在任何平台上都是固定的位数。CUDA 和 PyTorch 的源代码里到处是这些类型。
• 指针是 8 字节：在 64 位机器上，地址本身就是一个 64 位整数。
• size_t 是”大小”的专用类型（无符号）：malloc(n * sizeof(float))、cudaMalloc(&p, size) 的 size 都是它。
• 对照 Python：int 1 占 28 字节、2**100 占 40 字节、float 占 24 字节——Python 的整数可以任意大，代价是每个数都是一个对象。

2.2 一个数在内存里是哪几个字节

内存可以想象成一条很长的”字节街”，每个字节有一个门牌号（地址）。一个 4 字节的 int 占 4 个相邻的门牌。02 号脚本第 2 部分把几个数的原始字节读了出来：

• 本机是小端（little-endian）：低位字节放在低地址，所以 42 = 0x2a 存成 2a 00 00 00。几乎所有 CPU 和 GPU 都是小端；网络协议常用大端。读写二进制文件（例如模型权重的 .bin / .safetensors）时，格式里会约定字节序。
• −1 是四个 ff：这是补码表示（P2 §1）。
• 0.1 在 float 和 double 里都不是精确的 0.1，只是最接近 0.1 的二进制小数（P2 §2）。

2.3 读懂 CUDA 内核所需的最少 C 语法

对照着 Python 看：

最常让 Python 程序员出错的三个地方：整数除法会取整（1 / 2 是 0）；0.1 默认是 double，和 float 运算时会悄悄提升精度（在 GPU 上会变慢）；变量不初始化就使用，值是随机的垃圾（C 不会自动清零）。

3. 指针 · Pointers

3.1 地址与指针

直觉。 变量住在内存里，每个变量都有一个”门牌号”——地址。指针就是一个专门用来存门牌号的变量。

两个运算符：

• &x：取 x 的地址（“x 住在哪里”）。
• *p：取指针 p 指向的那个变量（“去 p 写的那个门牌号，看看里面是什么”）。

hello_pointers.c 的开头：

int x = 42;
int* px = &x;      // px 是一个"指向 int 的指针"，存的是 x 的地址
*px = 7;           // 通过 px 修改 x
printf("%d", x);   // 打印 7

运行时它打印 &x = 0x16d8e282c 之类的地址——每次运行都不同，因为操作系统每次把程序放在随机的位置（地址空间布局随机化，ASLR，一种安全措施）。我们关心的从来不是地址的绝对值，而是地址之间的差。

int* px 的写法读作”px 的类型是 int*“，也就是”指向 int 的指针”。float* 指向 float，void* 是”不知道指向什么类型”的通用指针（只能存地址，不能直接取值）。

3.2 数组与指针运算

C 的数组就是一段连续的内存，数组名可以当作首元素的地址来用。02 号脚本第 3 部分建了一个 float a[5]，从 C 里读回每个元素：

两条规则：

1. a[i] 就是 *(a + i) 的简写：先算地址，再取值。
2. 指针运算以”元素”为单位：p + 3 不是地址加 3，而是加 3 × sizeof(元素)。02 号脚本测得 float* 加 3 前进 12 字节，double* 加 3 前进 24 字节。指针”知道”自己指向什么类型，这正是它需要类型的原因。

所以 GPU 内核里的 A[i]，意思就是”从 A 指向的地址往后数 i 个 float，读出那个数”。

3.3 C 不检查越界

Python 的列表越界会抛出 IndexError；C 的 a[5]、a[-1] 都能编译，运行时会读写到别人的数据——可能读到垃圾、可能悄悄改坏别的变量、也可能让程序崩溃（“段错误” segmentation fault）。C 不检查，是为了速度：每次访问都检查边界是要花时间的。

这就是 GPU 内核开头那句 if (i < n) 的来历：线程是按块（例如每块 256 个）发射的，n 不一定是 256 的倍数，最后一块里会有一些线程的 i ≥ n，它们必须什么都不做。线程下标越界是 CUDA 程序里最常见的 bug。（02 号脚本没有真的做越界访问——那是”未定义行为”，结果不可预测，不适合当演示。）

3.4 按值传参与按指针传参

C 的函数参数总是复制一份传进去（按值传递）。02 号脚本第 4 部分：

void swap_by_value(int a, int b)   { int t = a; a = b; b = t; }      // 交换的是副本
void swap_by_pointer(int* a, int* b) { int t = *a; *a = *b; *b = t; } // 通过地址改原来的变量

指针传参还有一个作用：避免复制大数据。把一个 100 万元素的数组传给函数时，传的只是 8 字节的首地址，函数通过地址直接读写原数组。CUDA 内核没有返回值，结果全部通过”传进来的指针”写回显存，是同一个机制。

3.5 指向指针的指针：cudaMalloc(&dA, size)

现在可以回答开头的问题了。cudaMalloc 要做两件事：报告成功还是失败，以及把新分配的显存地址交给调用者。返回值被错误码占用了（cudaError_t，0 表示成功），地址只能写进调用者的变量里。要修改调用者手里的指针变量 dA，就必须拿到 dA 的地址——dA 本身是 float*，它的地址就是 float**，再统一转成通用的 void**。

02 号脚本第 5 部分用一个普通的 malloc 模仿了这个接口：

int my_malloc(void** out, size_t nbytes) {
    void* p = malloc(nbytes);
    if (p == NULL) return 2;     // 失败：返回错误码
    memset(p, 0, nbytes);
    *out = p;                    // 把地址写回调用者的变量
    return 0;                    // 成功
}

从 Python 调用：先建一个空指针变量 p（NULL），把它的地址传进去；调用后返回 0，p 变成了一个地址；再把这个地址当作 float* 交给 fill_squares 写入 0, 1, 4, 9, …，读回来验证；最后 my_free(p)。

和真正的 cudaMalloc 唯一的区别：cudaMalloc 得到的地址指向 GPU 显存，CPU 不能直接读写它，必须用 cudaMemcpy 在主机内存和显存之间搬运（§4.3）。

3.6 const 和 restrict

两个经常出现在内核签名里的修饰词：

• const float* A：承诺”不会通过 A 去修改它指向的数据”。编译器会检查（写了就报错），也可以据此做优化。只读的输入通常都标 const。
• float* __restrict__ C（CUDA）/ restrict（C99）：承诺”这个指针指向的内存不会和其他指针重叠”。没有这个承诺，编译器担心写 C[i] 会改到 A 的内容，就不敢把 A 的值缓存在寄存器里。性能敏感的内核里常见。

4. 内存的生命周期：栈与堆 · Stack and heap

4.1 栈：自动管理的短期内存

函数里声明的普通局部变量（int x;、float a[4];）放在栈（stack）上：进入函数时自动分配，函数返回时自动回收。栈很快，但很小（通常几 MB），而且大小在编译时就要确定。

一个经典的 bug：

float* make_array(void) {
    float a[100];      // 局部数组，在栈上
    return a;          // 错！函数返回后 a 的内存就被回收了
}

调用者拿到的是一个悬空指针（dangling pointer）：它指向的内存已经不属于你，随时会被别的函数覆盖。编译器通常会对这种写法给出警告。

4.2 堆：手动管理的长期内存

运行时才知道要多大、或者要在函数返回后继续使用的内存，从堆（heap）上用 malloc 申请，用完必须 free：

float* a = (float*)malloc(n * sizeof(float));   // 申请 n 个 float
if (a == NULL) { /* 申请失败 */ }
/* …… 使用 a …… */
free(a);       // 归还
a = NULL;      // 好习惯：避免之后误用

堆上的三种经典错误：忘了 free（内存泄漏）、free 之后还在用（use-after-free）、同一块 free 两次（double free）。Python 有垃圾回收，这些都替你管了；C 没有。

02 号脚本第 6 部分实测了泄漏（⏱）：分配 40 次 × 8 MiB 并写满，

在一个长时间运行的服务里，泄漏会一直涨，直到内存耗尽、进程被杀掉。

4.3 GPU 上的对应

GPU 有自己的内存（显存，device memory），和 CPU 的内存（host memory）是两块不同的物理内存，通过 PCIe 或 NVLink 相连。CUDA 提供了一一对应的接口：

两条要记住的规则：设备指针和主机指针都是”一个地址”，类型上看不出区别，混用（在 CPU 上解引用设备指针）会崩溃；主机和设备之间的拷贝很贵（PCIe 4.0 x16 约 32 GB/s，而 H100 的显存带宽是 3.35 TB/s，第 1 讲 §4.7），所以要尽量让数据留在 GPU 上。显存比内存小得多，“谁持有这块显存、什么时候释放”在推理服务里是核心问题（第 24 讲 vLLM 管理 KV 缓存）。

5. 多维数组：行优先与步长 · Row-major layout and strides

5.1 一维内存里的二维矩阵

内存只有一维，矩阵却有两维。最常见的排法是行优先（row-major，C、NumPy、PyTorch 的默认）：第 0 行的元素依次排完，接着排第 1 行……03 号脚本第 1 部分看了一个 3×4 的 float32 矩阵：

  a =   0  1  2  3        内存里的顺序： 0 1 2 3 | 4 5 6 7 | 8 9 10 11
        4  5  6  7
        8  9 10 11

元素 $(i, j)$ 的一维位置是 $i \times 4 + j \times 1$。这两个系数叫步长（strides）：行号加 1 要跳过 4 个元素，列号加 1 跳过 1 个。NumPy 的 a.strides 是 (16, 4)——它以字节为单位（每个 float32 4 字节），除以 4 就是元素为单位的 (4, 1)。03 号脚本用一个只拿到”首地址 + 两个步长”的 C 函数 get2d(base, i, j, 4, 1) 读出了和 NumPy 完全相同的元素。

5.2 步长的一般公式

$n$ 维张量的元素 $(i_0, i_1, \dots, i_{n-1})$ 的位置：

行优先连续存放时，最后一维的步长是 1，往前每一维的步长等于后面所有维度大小的乘积。例如形状 (2, 3, 4) 的步长是 (12, 4, 1)：第 0 维加 1 要跳过一整个 3×4 的”面”。tensor_layout.py 的 strides_for 就是这个计算。

5.3 转置 = 交换步长

直觉。 转置不需要搬动任何数据：只要换一种”读法”。03 号脚本第 2 部分：

a.T 和 a 共用同一块内存，只是形状和步长交换了。C 函数用同一个指针、步长 (1, 4) 读出的 4×3 矩阵正好等于 a.T。

代价是 a.T 不再连续：按行读它时，相邻元素在内存里相隔 4 个位置。很多 C / CUDA 函数假设输入是连续的，这时要么先复制成连续的（ascontiguousarray、PyTorch 的 .contiguous()），要么把步长也传进去。

什么时候必须复制？ 把 a.T 展平成 12 个元素，顺序是 0 4 8 1 5 9 2 6 10 3 7 11——没有任何一组步长能表达这个顺序，所以 a.T.reshape(12) 只能复制一份（03 号脚本验证了结果与 a 不共享内存）。PyTorch 的 view() 遇到这种情况会直接报错，reshape() 会悄悄复制——很多人第一次写 Transformer 时都在这里踩过坑（多头注意力里的 transpose 之后接 view，第 6 讲）。

5.4 切片与广播也是视图

• 切片 = 新的起点偏移 + 新的形状 + 新的步长；隔一个取一个就是步长翻倍。写切片就是写原数组（03 号脚本里改 s1[0, 0]，a[0, 1] 也变了）。
• 广播 = 某一维的步长为 0：沿着那一维走，地址不动，同一个元素被重复使用。偏置加法 $y = xW + b$ 就是这样把长度为 d 的向量 b 加到每一行上的，不需要把 b 真的复制成矩阵（广播视图只占 v 本身的 16 字节，而真正的 3×4 数组要 48 字节）。

5.5 原课 HW1 的下标函数

原课 HW1 的提示里专门有一节”数据布局与步长”：要在 CUDA 里实现 map、zip、reduce，每个线程首先要知道自己负责哪个元素、它的输入在哪里。miniTorch 用三个小函数完成这件事（下面是讲解用的参考版本，原课作业要你自己写）：

• to_index(ordinal, shape)：按行优先顺序数的第 ordinal 个元素 → 多维下标。从最后一维开始，反复”除以这一维的大小，余数就是这一维的下标”。例：shape (2, 3, 4) 的第 17 个元素 → (1, 1, 1)，因为 17 = 1×12 + 1×4 + 1。GPU 线程拿到的是一维编号 i = blockIdx.x * blockDim.x + threadIdx.x，第一步就是 to_index。
• index_to_position(index, strides)：多维下标 → 存储里的位置（§5.2 的公式）。
• broadcast_index(big_index, big_shape, shape)：广播后的大张量里的下标 → 原来的小张量里的下标：小张量和大张量的最后几维对齐，大小为 1 的维度下标取 0，多出来的前几维丢掉。例：把 (3, 1) 广播到 (2, 3, 4) 时，大下标 (1, 2, 3) 对应小下标 (2, 0)。

03 号脚本第 4 部分把它们和 NumPy 逐个对了答案：4 种形状下 to_index 再 index_to_position 回到原编号（71/71）；(2, 3, 4) 的 6 种维度排列下按步长读取等于 NumPy 的 transpose（144/144）；广播加法 (2, 1, 4) + (3, 1) → (2, 3, 4) 的每个输出元素，用 broadcast_index 找到两个输入后相加，等于 NumPy 的结果（24/24）。

于是一个”逐元素 zip”的 GPU 内核在逻辑上就是这样（伪代码）：

线程 i:
    if i >= out_size: return
    out_index = to_index(i, out_shape)
    a_index   = broadcast_index(out_index, out_shape, a_shape)
    b_index   = broadcast_index(out_index, out_shape, b_shape)
    out[index_to_position(out_index, out_strides)] =
        f(a[index_to_position(a_index, a_strides)], b[index_to_position(b_index, b_strides)])

第 3 讲会把它真正写出来（在 SIMT 模拟器上运行）。

6. 结构体与 Python ↔ C · Structs and calling C from Python

6.1 结构体：把相关的数据打包

结构体（struct）把几个字段按顺序排在一起。描述一个二维张量视图，需要一个数据指针、形状和步长：

typedef struct {
    float* data;      // 指向第一个元素
    long rows, cols;
    long stride_row;  // 行号加 1 时跳过多少个元素
    long stride_col;  // 列号加 1 时跳过多少个元素
} Tensor2D;

float tensor_get(const Tensor2D* t, long i, long j) {
    return t->data[i * t->stride_row + j * t->stride_col];   // t->x 等价于 (*t).x
}

03 号脚本第 5 部分用同一个数据指针构造了两个 Tensor2D，交给同一个 C 函数 tensor_row_sums：

换一组步长，“按行求和”就变成了”按列求和”，一个字节都没有复制。miniTorch 的 TensorData（存储 + 形状 + 步长）就是这个结构体的 Python 版本；PyTorch 的张量内部也是”存储 + 尺寸 + 步长 + 偏移”。（这个结构体占 40 字节：1 个指针 + 4 个 long。）

6.2 用 ctypes 从 Python 调用 C

c_lab.py 演示了完整的流程：

1. 编译成共享库：cc -O2 -shared -fPIC -o build/p1lib.dylib csrc/p1lib.c（Linux 上是 .so）。共享库是”可以被别的程序在运行时加载的函数集合”。
2. 加载：lib = ctypes.CDLL(路径)。
3. 声明签名：lib.sum_floats.argtypes = [POINTER(c_float), c_long]，lib.sum_floats.restype = c_float。ctypes 默认把所有参数当 int、返回值当 int；不声明的话，传 float 或返回 double 时会得到莫名其妙的数字（甚至崩溃）。
4. 传 NumPy 数组：arr.ctypes.data_as(POINTER(c_float)) 取出数组数据区的地址（不复制）。前提是 dtype 对得上（float32 ↔ float），并且 C 用完之前数组不能被释放。
5. 传”变量的地址”：ctypes.byref(x) 相当于 C 里的 &x（§3.4、§3.5）。
6. 传结构体：Python 侧定义一个字段顺序和类型完全相同的 ctypes.Structure。

PyTorch 的 C++ 扩展、miniTorch 加载 CUDA 内核，用的是更高级的工具（pybind11 等），但原理一样：Python 拿着数据的地址和几个整数，交给编译好的代码去算。

6.3 为什么还是要用库

03 号脚本第 6 部分最后做了两个实验（⏱ 本机）：

• 同样是 $2n^3$ 次运算，专业的 BLAS 库比朴素三重循环快两个多数量级：它做了分块（让数据留在缓存里）、向量化，并用上了专用的矩阵运算单元。第 4 讲会在 GPU 上一步步重做这些优化，你会看到为什么实际工作中总是调用 cuBLAS。（这一组输入上两种实现的结果碰巧逐位相同；一般来说，加法顺序不同会带来最后几位的差异。）
• 按列遍历读的是同样多的数，却慢了 5.7 倍：CPU 每次从内存取来的是一整块（缓存行，本机 128 字节 = 32 个 float），按行读时整块都用上了，按列读时每次只用其中一个。内存访问的顺序决定速度——P4 会讲存储层次，第 4 讲会讲 GPU 上的”合并访存”（coalescing），本质是同一件事。

7. 回到开头：逐行读懂一个 CUDA 程序

现在可以逐行读 §0.1 的那段代码了：

唯一还没讲的是 __global__、<<<…>>>、blockIdx、threadIdx——它们是 CUDA 在 C 之上加的扩展，正是第 2 讲的内容。

8. 小结 · Summary

• 编译型 vs 解释型：C 在运行前一次翻译成机器码（预处理 → 汇编 → 目标文件 → 链接）；Python 边解释边执行，每条字节码都要判断类型、每个数都是对象。实测求和时 C 比 Python 快约 14 倍，NumPy 快约 71 倍——NumPy 用多个部分和打破了依赖链，这正是并行归约的思想。换一种加法顺序，浮点结果就会不同。
• 跨语言调用有固定开销（约 140 纳秒一次）：要一次处理一大批数据，GPU 内核同理。
• 类型决定字节数和解释方式：int 4 字节、float 4 字节、double 8 字节、指针 8 字节；需要确定位数时用 int32_t 等。本机是小端。
• 指针存地址：& 取地址、* 取值；a[i] ≡ *(a + i)，指针运算以元素为单位；C 不检查越界，所以内核要写 if (i < n)；要修改调用者的变量就传它的地址；cudaMalloc(&ptr, size) 用 void** 把地址写回调用者。
• 栈上的局部变量自动回收，堆上的内存要 malloc / free 成对；忘了 free 就泄漏（实测 40 × 8 MiB）。GPU 显存有对应的 cudaMalloc / cudaFree / cudaMemcpy，主机和设备之间的拷贝很贵。
• 张量 = 一维存储 + 形状 + 步长：位置 = Σ 下标 × 步长。转置、切片、广播（步长为 0）都是零拷贝的视图；视图不连续时，有的操作（展平、某些内核）必须复制。HW1 的 to_index、index_to_position、broadcast_index 就是把线程编号翻译成内存位置的全部工具。
• 结构体把”指针 + 形状 + 步长”打包交给 C；ctypes 让 Python 直接调用 C。专业库（BLAS / cuBLAS）比朴素实现快几百倍；内存访问的顺序决定速度。

9. 常见误区 · Common pitfalls

1. “C 一定比 Python 快一百倍。” 取决于 C 那一侧写得怎样。本讲的顺序求和 C 只快 14 倍，因为它被依赖链卡住了；NumPy 的成对求和又比它快 5 倍。
2. “打印出来的地址每次都一样，可以用来判断两个变量是不是同一个。” 绝对地址每次运行都不同（ASLR）。判断两个数组是否共享内存，用 np.shares_memory 或比较同一次运行里的地址。
3. “p + 1 让地址加 1。” 加的是 1 个元素的大小：float* 加 4 字节，double* 加 8 字节。
4. “函数里改了参数，调用者的变量就变了。” C 按值传参，改的是副本；要改调用者的变量，必须传它的地址。
5. “转置会复制数据。” 转置只交换步长；但之后如果需要连续的数据（展平、调用要求连续输入的内核），才会发生复制。PyTorch 里 x.transpose(1, 2).view(...) 会报错，要用 .contiguous().view(...) 或 .reshape(...)。
6. “广播会把小张量复制成大张量。” 广播是步长为 0 的视图，不占额外内存；只有在你真的把结果写出来时，输出才占满整块。
7. “1 / 2 等于 0.5。” 在 C 里两个整数相除结果是 0（取整）。(n + 255) / 256 正是利用这一点实现”向上取整”。
8. “设备指针和主机指针可以混用。” 它们都只是地址，类型上看不出区别；在 CPU 上解引用设备指针会崩溃，必须用 cudaMemcpy 搬运。
9. “ctypes 能自动猜出参数类型。” 不能。不声明 argtypes / restype，float 和 double 会被当成 int，得到错误的结果。

10. 练习 · Exercises

练习 1（字节数）。 (a) 一个有 100 万个 double 的数组占多少字节？(b) 一个指向它的指针占多少字节？(c) 一个 70 亿参数的模型，权重存成 float32 要多少 GB？

参考答案：(a) $10^6 \times 8 = 8 \times 10^6$ 字节 = 8 MB。(b) 8 字节（64 位机器上任何指针都是 8 字节，和它指向多大的数组无关）。(c) $7 \times 10^9 \times 4 = 28 \times 10^9$ 字节 = 28 GB（这就是为什么要用 16 位甚至更低的精度，P2、第 19 讲）。

练习 2（指针运算）。 假设 float* p 的值是 0x1000。p + 5 的值是多少？如果 p 是 double* 呢？(char*)p + 5 呢？

参考答案：float*：0x1000 + 5×4 = 0x1014；double*：0x1000 + 5×8 = 0x1028；char*：0x1005（char 是 1 字节）。这也是 02 号脚本里 ptr_diff_bytes 先把指针转成 char* 再相减的原因。

练习 3（按值传参）。 下面的函数想把数组每个元素乘 2，并把数组长度改成 0 表示”处理完了”。哪一部分生效、哪一部分不生效？怎样修改？

void process(float* a, int n) {
    for (int i = 0; i < n; i++) a[i] *= 2;
    n = 0;
}

参考答案：乘 2 生效（通过指针 a 修改的是调用者的数组）；n = 0 不生效（n 是副本）。要让调用者的长度变成 0，参数要改成 int* n，函数里写 *n = 0，调用时传 &n。

练习 4（void**）。 仿照 cudaMalloc，写一个 C 函数 int alloc_floats(float** out, int n)：申请 n 个 float、全部置 1、把地址写进 *out、成功返回 0。再写出调用它的代码。

参考答案：

int alloc_floats(float** out, int n) {
    float* p = (float*)malloc(n * sizeof(float));
    if (p == NULL) return 1;
    for (int i = 0; i < n; i++) p[i] = 1.0f;
    *out = p;
    return 0;
}
/* 调用：*/
float* buf = NULL;
if (alloc_floats(&buf, 100) == 0) { /* 使用 buf …… */ free(buf); }

练习 5（行优先位置）。 形状 (2, 3, 4) 的连续张量，元素 (1, 2, 3) 的一维位置是多少？反过来，位置 23 对应哪个下标？

参考答案：步长 (12, 4, 1)，位置 = 1×12 + 2×4 + 3×1 = 23；to_index(23, (2, 3, 4)) = (1, 2, 3)。（23 是最后一个元素，因为总共 24 个。）

练习 6（维度排列）。 对形状 (2, 3, 4)、步长 (12, 4, 1) 的张量做 permute((2, 0, 1))，新的形状和步长是多少？它连续吗？

参考答案：新的第 k 维 = 原来的第 order[k] 维：形状 (4, 2, 3)，步长 (1, 12, 4)。连续张量的步长应该是 (6, 3, 1)，所以不连续。

练习 7（广播）。 (a) 形状 (4, 1) 和 (1, 5) 广播后是什么形状？(b) np.broadcast_to(v, (4, 5))（v 的形状是 (5,)、float32）的 .strides 是多少（字节）？(c) 广播到 (4, 5) 后，大下标 (3, 2) 对应 (4, 1) 张量的哪个下标、(1, 5) 张量的哪个下标？

参考答案：(a) (4, 5)。(b) (0, 4)：第 0 维步长 0（每一行复用同一个 v），第 1 维 4 字节。(c) (4, 1) 张量：(3, 0)；(1, 5) 张量：(0, 2)。

练习 8（找 bug）。 下面的函数有什么问题？怎样修改？

float* ones(int n) {
    float a[1024];
    for (int i = 0; i < n; i++) a[i] = 1.0f;
    return a;
}

参考答案：两个问题。① 返回了栈上局部数组的地址，函数返回后这块内存就被回收，调用者拿到的是悬空指针。② 当 n > 1024 时会越界写。修改：用 float* a = (float*)malloc(n * sizeof(float)); 在堆上分配（并检查是否为 NULL），由调用者负责 free。

练习 9（泄漏）。 一个推理服务每处理一个请求就 malloc 一块 2 MB 的缓冲区存中间结果，但忘了 free。每秒 100 个请求，一台 64 GB 内存的机器大约多久会耗尽内存？

参考答案：每秒泄漏 200 MB，64 GB / 0.2 GB/s = 320 秒，约 5 分钟。泄漏的可怕之处在于它在测试时（几个请求）看不出来，上线后才爆发。

练习 10（求和的精度）。 01 号脚本里，用 float 顺序累加 100 万个 [0, 1) 的随机数差了 3.83，NumPy 只差 0.02。用你自己的话解释为什么，并说出两种不换成 double 也能提高精度的办法。

参考答案：顺序累加时，部分和很快长到几十万，而 float 只有大约 7 位有效数字，几十万量级的数再加上一个零点几的数，最后几位要被舍入掉，误差一步步累积。NumPy 的成对求和让每次相加的两个数大小相近，舍入误差小得多。不换 double 的办法：① 成对 / 分块求和（树形归约，GPU 上的 reduce 就是这样）；② Kahan 补偿求和（用一个额外的变量记住每一步被舍掉的部分，P2）。

练习 11（局部性）。 03 号脚本里按列求和比按行慢 5.7 倍。如果你确实需要一个行优先矩阵的每一列的和，怎样写才快？

参考答案：仍然按行遍历，但维护一个长度为 cols 的累加数组：对每一行，把这一行的元素逐个加到 col_sum[j] 上。这样内存访问仍然是连续的，结果就是每一列的和。（另一种办法是先转置成连续的再按行求和，但要多一次复制。GPU 上做按列归约时也会用类似的”换一种遍历顺序”的技巧。）

练习 12（HW1 预习）。 用 §5.5 的伪代码，手算一个 zip 内核的一个线程：输出形状 (2, 3)、输入 a 的形状 (2, 1)（连续）、b 的形状 (3,)（连续），线程 i = 4。它读 a 和 b 的哪个位置、写输出的哪个位置？

参考答案：to_index(4, (2, 3)) = (1, 1)。a：broadcast_index((1, 1), (2, 3), (2, 1)) = (1, 0)，步长 (1, 1)，位置 1；b：broadcast_index((1, 1), (2, 3), (3,)) = (1,)，位置 1；输出步长 (3, 1)，位置 1×3 + 1 = 4。即 out[4] = f(a[1], b[1])。

11. 延伸阅读 · Further reading

• C 语言：Kernighan & Ritchie,《C 程序设计语言》（第 2 版）第 5 章（指针与数组）、第 6 章（结构体）；Bryant & O’Hallaron, Computer Systems: A Programmer’s Perspective（CS:APP）第 2–3 章（数据的机器级表示）。
• CUDA 的内存管理接口：NVIDIA CUDA C++ Programming Guide 的 “Programming Interface → CUDA Runtime → Device Memory” 一节；PMPP 第 2 章（向量加法的完整例子）。
• 步长：NumPy 文档 “Internal memory layout of an ndarray”；Edward Z. Yang 的博客文章 PyTorch internals（讲 PyTorch 张量的存储、尺寸、步长和视图）。
• 原课 HW1：作业说明里”Data layout and strides”一节（https://llmsystem.github.io/llmsystemhomework/assignment_1/ ）。
• ctypes：Python 官方文档 ctypes — A foreign function library for Python。

12. 下一篇预告 · Next

本讲多次碰到”浮点数的舍入误差”：float 顺序累加差了 3.83、0.1 不能精确表示、换一种加法顺序结果就变。前置讲 P2（数字的表示） 会从比特讲起：补码与整数溢出，IEEE 754 浮点数的符号、指数、尾数，FP32、TF32、FP16、BF16、FP8 各自的范围与精度，浮点运算为什么不满足结合律，混合精度训练为什么需要 FP32 主权重和损失缩放，以及量化（第 19–20 讲）的基本思路。