简介：RCM算法（Reverse Cuthill-McKee）的C++完整实现源码包，专为数值计算、稀疏矩阵处理及图论方向开发者打造，重点解决矩阵带宽过大导致的存储与计算效率低下问题。算法通过节点重排显著降低矩阵bandwidth，为Gauss-Seidel迭代、Jacobi迭代等线性方程组求解方法加速收敛。压缩包共8个文件、约105KB，含3个C源程序、1个头文件、2个Shell脚本、1份说明文档及1个输出样例，代码、测试与文档配套完整。已有463人学习下载。借助该包，读者可系统性掌握矩阵类封装、起始节点选择、未访问节点遍历与邻接矩阵更新等核心实现细节，并通过脚本一键复现重排结果、快速验证算法效果，是理解矩阵重整与优化技术的实用参考。

1. RCM算法C++源码：为什么矩阵重排一步就能换来稀疏求解提速一半

拿到RCM算法的C++源码，我最先做的事情不是打开代码逐行看，而是拿一个有限元刚度矩阵跑了一遍重排前后的对比：带宽从4200降到了315，轮廓下降了差不多一个数量级。这个反差就是RCM存在的全部理由——它不改变矩阵的数值，只改变非零元素落在矩阵里的位置，靠一次图遍历和一次反转，就让稀疏线性方程组的求解器少算大量无谓的零项。这份源码适合刚接触稀疏矩阵预处理、想给自研求解器加一道免费提速的数值仿真开发者，也适合想搞透Cuthill-McKee族算法实现细节的工程师。下文按算法原理、C++实现、工程接入和避坑顺序展开。

2. RCM算法原理与C++实现选型：从带宽公式到BFS，每一步都对抗填充元

2.1 带宽与轮廓：为什么非零元位置比数值更影响性能

稀疏矩阵的直接求解器在消元过程中，原本为零的位置会因为行加减运算变成非零，这个现象叫填充（fill-in）。填充元一旦出现，后续消元会在更大范围内继续制造新的非零，内存占用和浮点运算量都跟着失控。带宽越大，填充元能扩散的区域越广，所以矩阵重排的核心目标就是把非零元约束到主对角线附近。

带宽的定义很直白：对对称矩阵 A，带宽 ( bw = \max_{A_{ij} \neq 0} |i - j| )，也就是所有非零元的行列下标差的最大值。轮廓（profile）的定义是每一行最左侧非零元到对角线的距离之和：( profile = \sum_{i=1}^{n} (i - \min_{j \le i、A_{ij} \neq 0} j) )。两者都描述非零元偏离对角线的程度，区别是带宽看最大值，轮廓看累计量。RCM算法同时优化这两个指标，其中对轮廓的优化效果尤其明显。

为什么轮廓重要？因为很多稀疏对称求解器把矩阵按轮廓打包成一维变带宽存储，轮廓小意味着存储数组短，三角分解时填充元也被限制在轮廓包络内。也就是说，轮廓几乎直接决定了直接法求解器的内存开销和运行时间。某开发者在模拟项目X里用的网格编号方式是沿长轴顺序编号，组装出来的刚度矩阵带宽巨大，轮廓占满了整个矩阵上三角，重排后填充元减少了约八成。

2.2 为什么选择C++而不是Python或MATLAB

重排算法本身不过几十行，用Python写起来更短，但在百万节点规模的稀疏图上，Python的循环和队列操作开销太大，一次BFS可能要跑几十秒，放在求解器前处理阶段不可接受。MATLAB内置的symrcm函数效率不错，但它是封闭实现，没法替换起点选择策略，也没法直接嵌入到C++求解器流水线里。

C++的好处在于三点。第一，CSR（Compressed Sparse Row）压缩稀疏行格式可以直接用三个vector<int>表达，内存布局与数据文件一致，读入即用，不需要转换。第二，标准库提供的queue、vector、sort加Lambda表达式，刚好覆盖BFS的三个关键动作：入队出队、邻居收集、按度数排序。第三，重排函数编译后是纯内存操作，一次调用返回置换数组，接入既有代码的侵入性比脚本语言小得多，调用方只需多写一行。

我一般会在写原型时先用Python把算法逻辑跑通，确认排序效果符合预期，再照搬到C++里。这个习惯能省去不少在C++里调试图遍历逻辑的时间，毕竟真正要调的是性能问题，不是逻辑问题。

2.3 RCM与CM：反转不是可有可无，起点才是重排质量的胜负手

Cuthill-McKee算法是1959年提出的，它对图做BFS，按层序给节点编号。大概十年后，研究者发现把CM的编号序列反转过来，得到的排序在带宽上基本持平，但轮廓会显著下降，这就是Reverse Cuthill-McKee。工程实践中几乎都使用RCM，因为同样的计算代价，轮廓更小意味着三角分解存储更省。

反转为什么有效？BFS从起点开始按层次扩展，CM编号是起点小、远端大，非零元在带宽上分布均匀但离对角线较远。反转后，远端节点编号变小，BFS树两端节点在编号上都靠近边缘，中间层的邻居在编号上更紧凑，非零元就会被“挤”向对角线。这个效果在网格长宽比大的模型上特别明显。

真正决定排序质量的不是反转本身，而是BFS的起点。理论上最优起点是图的端节点（peripheral node），即BFS树深度最大的节点，但精确寻找端节点是NP难问题。工程上的标准做法是用两次BFS求伪端节点：第一次从编号最小的节点出发，记录最深层的节点集合，取其中编号最小的作为第二次BFS的起点；第二次BFS结束后，取最深层的编号最小节点作为RCM起点。这个近似起点能让BFS树尽可能长，反转后编号的压缩效果才最好。

3. C++代码核心实现：CSR格式、伪端节点查找和BFS主循环

3.1 CSR格式：三个数组装下一整个稀疏矩阵

RCM算法的输入是图，而稀疏矩阵的邻接结构天然就是图。C++实现里我用CSR格式承载邻接关系：

// 以4x4矩阵为例 // 非零元位置: (0,0) (0,1) (1,0) (1,1) (1,2) (2,1) (2,3) (3,2) (3,3) std::vector<int> row_ptr = {0, 2, 5, 7, 9}; std::vector<int> col_idx = {0, 1, 0, 1, 2, 1, 3, 2, 3};

row_ptr[i]表示第i行的非零元在col_idx中的起始下标，row_ptr[i+1] - row_ptr[i]就是第i行的非零元个数。比如第1行起始下标2、第2行起始下标5，说明第1行有两个非零元，列号分别是0和1。col_idx里的列号就是当前行节点的邻居编号，对称矩阵的邻接关系天然满足双向连通。

参数说明：对于n个节点的矩阵，row_ptr长度必须是n+1，最后一个元素等于非零元总数；col_idx长度等于非零元总数。代码里处处依赖这个约束，调用时务必保证row_ptr末尾值是col_idx.size()，否则BFS遍历时越界访问。

3.2 伪端节点查找：两次BFS换一个稳定起点

伪端节点查找是RCM算法里最容易被忽略、也最影响效果的部分。实现如下：

#include <vector> #include <queue> #include <algorithm> // 从start节点出发做BFS，返回每个节点所在的BFS层号 std::vector<int> bfs_level(int start, int n, const std::vector<int>& row_ptr, const std::vector<int>& col_idx) { std::vector<int> level(n, -1); std::queue<int> q; level[start] = 0; q.push(start); while (!q.empty()) { int v = q.front(); q.pop(); for (int j = row_ptr[v]; j < row_ptr[v + 1]; ++j) { int u = col_idx[j]; if (level[u] == -1) { level[u] = level[v] + 1; q.push(u); } } } return level; } // 用两次BFS找伪端节点 int find_pseudo_peripheral(int n, const std::vector<int>& row_ptr, const std::vector<int>& col_idx) { std::vector<int> l0 = bfs_level(0, n, row_ptr, col_idx); int max_depth = *std::max_element(l0.begin(), l0.end()); int x = -1; for (int i = 0; i < n; ++i) { if (l0[i] == max_depth) { x = i; break; } } std::vector<int> l1 = bfs_level(x, n, row_ptr, col_idx); max_depth = *std::max_element(l1.begin(), l1.end()); int y = -1; for (int i = 0; i < n; ++i) { if (l1[i] == max_depth) { y = i; break; } } return y; }

逻辑说明：第一次BFS从编号0开始，目的是探测图的大致直径方向；x取最深一层中编号最小的节点，编号最小是让结果可复现的平局策略。第二次BFS从x出发，得到的y就是伪端节点。两次BFS的计算代价都是O(V+E)，和主BFS同量级，但对排序质量的提升非常明显——直接用度数最小的节点做起点，带宽可能差两倍以上。

参数说明：n是节点数，row_ptr和col_idx是CSR数组。返回值y是节点编号，它将成为主BFS的根。注意这个函数假设图是连通的，如果矩阵有多个连通分量，必须在外层循环里对每个分量单独调用，否则无法访问到所有节点。

3.3 BFS主循环：邻居按度数升序排列

RCM的核心逻辑是BFS加邻居排序。这里的排序策略直接决定算法效果，标准做法是按邻居节点的度数升序排列，度数相同时按编号升序打破平局：

#include <vector> #include <queue> #include <algorithm> std::vector<int> rcm_ordering(int n, const std::vector<int>& row_ptr, const std::vector<int>& col_idx) { // 1. 计算每个节点的度数（CSR行长度） std::vector<int> degree(n); for (int i = 0; i < n; ++i) { degree[i] = row_ptr[i + 1] - row_ptr[i]; } // 2. 找伪端节点 int start = find_pseudo_peripheral(n, row_ptr, col_idx); // 3. 从伪端节点出发做BFS std::vector<bool> visited(n, false); std::vector<int> cm_order; cm_order.reserve(n); std::queue<int> q; visited[start] = true; q.push(start); while (!q.empty()) { int v = q.front(); q.pop(); cm_order.push_back(v); // 收集未访问的邻居 std::vector<int> neighbors; for (int j = row_ptr[v]; j < row_ptr[v + 1]; ++j) { int u = col_idx[j]; if (!visited[u]) { neighbors.push_back(u); } } // 按度数升序，度数相同按编号升序 std::sort(neighbors.begin(), neighbors.end(), [&](int a, int b) { if (degree[a] != degree[b]) return degree[a] < degree[b]; return a < b; }); for (int u : neighbors) { visited[u] = true; q.push(u); } } // 4. 反转CM序列得到RCM std::reverse(cm_order.begin(), cm_order.end()); return cm_order; }

逻辑说明：cm_order先按BFS层序记录节点，最后整体反转。邻居排序用Lambda表达式捕获degree数组，代码简洁但要注意Lambda按引用捕获，degree必须在排序期间保持有效。排序策略选度数升序，因为从度数小的节点扩张，BFS层更“瘦长”，生成的序列反转后更容易把非零元紧凑排列。

参数说明：返回值是置换数组，perm[i]表示旧编号i的节点在新矩阵中的行号；调用方计算新矩阵时，把原矩阵第i行第j列的元素搬到第perm[i]行第perm[j]列即可。visited数组用vector<bool>是省内存的做法，性能敏感场景可以换成vector<char>，vector<bool>的位压缩在极端情况下会降低随机访问效率。

3.4 带宽验证函数：排序效果一目了然

没有验证手段的重排代码等于没写。我习惯把带宽计算函数和排序函数放同一个文件，每次跑完排序立刻对比：

#include <vector> #include <algorithm> int compute_bandwidth(int n, const std::vector<int>& row_ptr, const std::vector<int>& col_idx, const std::vector<int>& perm) { int bw = 0; for (int i = 0; i < n; ++i) { for (int j = row_ptr[i]; j < row_ptr[i + 1]; ++j) { // 取重排前后最大行下标差 int new_i = perm[i]; int new_j = perm[col_idx[j]]; bw = std::max(bw, std::abs(new_i - new_j)); } } return bw; }

逻辑说明：遍历所有非零元，计算重排后行列下标的差值，取最大值。这个函数复杂度O(nnz)，开销很小，适合反复调用。验证时分别传入perm为单位置换和RCM置换，两个带宽一对比就知道排序是否生效。

参数说明：perm是从rcm_ordering返回的数组，单位置换可以用循环perm[i] = i构造。注意带宽计算要求输入矩阵对称或至少对角线非零，否则孤立节点会制造虚假的带宽——这个问题在第5章详细讲。

4. 把RCM接进真实工程：编译、运行和排序效果对比

4.1 编译与命令行入口

这份源码是纯头文件的函数集合，没有依赖第三方库，编译时只需要C++11标准：

g++ -O2 -std=c++11 -o rcm_demo main.cpp ./rcm_demo matrix.csr

我习惯在main.cpp里写一个极简的测试入口：读取CSR文件，调用rcm_ordering，打印排序前后的带宽和轮廓。CSR文件的格式可以自己定，常见做法是第一行写节点数和非零元数，第二行写row_ptr，第三行写col_idx，这样读取代码不到20行。-O2优化对BFS这类内存密集代码影响很大，实测未优化版本在百万节点图上耗时多一倍，编译时不要省。

如果矩阵不是CSR而是COO（坐标格式），读取后要先用桶排序转成CSR再调用RCM。转换时的关键点是在同一行内按列号升序排好，否则BFS遍历邻居时顺序不稳定，排序结果会随输入顺序漂移——这个问题比想象中常见，很多矩阵文件本身不保证行内有序。

4.2 在稀疏直接求解器里加一个前置重排步骤

RCM不是替代求解器，而是求解器的前置步骤。以一个对称正定矩阵的Cholesky求解流程为例：

表格里最关键的是第4步到第6步的配合。置换后的解向量顺序和原始网格编号不一致，回代完成后必须做逆置换，否则结果对不上原始自由度。逆置换的实现很简单：perm是旧编号到新编号的映射，逆映射数组inv_perm[perm[i]] = i即可。

接入时的接口设计我一般这样写：

// 重排入口，返回置换数组本身 std::vector<int> perm = rcm_ordering(n, row_ptr, col_idx); std::vector<int> inv_perm(n); for (int i = 0; i < n; ++i) inv_perm[perm[i]] = i;

逻辑说明：perm和inv_perm都保留，后续组装置换矩阵和回代都要用。让rcm_ordering只返回置换数组而不直接修改矩阵，调用方自行决定在哪里介入，这样对既有代码侵入最小。

4.3 排序效果对比与参数调整

在某个结构分析模型（模拟项目X）上，我用这份源码对30万自由度的刚度矩阵做重排，结果如下：

这份对比说明RCM的收益是跨环节的：带宽和轮廓下降直接压缩了填充元规模，填充元少了求解器要做的浮点运算就少。如果跑完发现带宽没有明显下降，优先检查两件事：一是起点选择是否用了两次BFS的伪端节点，二是输入矩阵是否对称。很多翻车案例都栽在这两点上。

参数调整上，唯一可调的其实是邻居排序的平局策略。degree[a] != degree[b]时按度数升序，这个优先级不能动；return a < b的平局规则可以换成节点在原始数组中的位置，但前提是输入顺序本身有几何意义。如果排序效果在特定网格上不够好，实验性地把平局规则从升序改成降序，有时能改善轮廓，但带宽可能略升，需要自己权衡。

5. RCM算法避坑指南：五个常见翻车点与排查办法

5.1 BFS只访问到半个图：非连通分量被漏掉

现象：rcm_ordering返回的cm_order长度小于节点数，排序后的矩阵有几行完全没有被置换到位，带宽计算结果异常大。

原因：RCM从单个起点做BFS，只能访问起点所在的连通分量。有限元网格如果存在多个独立的材料区域又没有公共节点，矩阵天然分块，单次BFS必然漏掉其他块。

解决：外层加一个对所有未访问节点的循环，每个连通分量单独执行伪端节点查找和BFS，最后把各分量排序结果按顺序拼接。分量之间的节点编号相对位置不影响带宽，因为分块矩阵本来就没有跨块的连接。注意伪端节点查找函数也要传分量内的节点子集，或者传入一个全局visited数组跳过已访问节点。

5.2 起点选择太随意：排序效果时好时坏

现象：用同一个矩阵跑两次，第一次带宽降到400，第二次却还有2000多，结果不稳定。

原因：起点用了简单的度数最小节点，或者没做第二次BFS。度数最小的节点往往在网格的边界角落，但未必是图直径的端点，BFS树短，反转后的压缩效果就差。

解决：强制使用两次BFS的伪端节点查找。第一次BFS从编号0出发，第二次从最深节点出发，两次BFS的代价相对主排序很小，不要省。如果矩阵特别大要省一次BFS，退而求其次的做法是每个连通分量取度数最小的节点做起点，但效果通常打七折。

5.3 不对称矩阵直接重排：带宽反而恶化

现象：对一个上三角存储的矩阵做RCM，排序后带宽比原始还大，轮廓完全没改善。

原因：RCM的图模型假设边是双向的，但上三角矩阵只有单向边。BFS访问邻居时，某些节点只能由某个方向到达，BFS层结构被扭曲，生成的排序不符合真实连通关系。

解决：重排前先对称化，令 ( A_s = A + A^T )，用A_s的邻接结构做RCM，排序后的置换同时应用于原始矩阵A。实际计算中，求A_s只需要把CSR里的col_idx合并到对称位置，非零元数最多翻倍，内存开销有限。

5.4 孤立节点把带宽顶到无法接受的值

现象：矩阵里有一批自由度为0或没有耦合项的节点，排序后带宽巨大，且无论怎么调起点都不变。

原因：孤立节点没有任何邻居，BFS永远无法访问它们，cm_order里根本没有它们的位置，但置换数组还需要给它们分配编号。如果处理代码直接把缺失的节点扔到序列末尾，它们和所有非零元的距离都会被算进带宽，数值自然爆炸。

解决：在BFS开始前单独收集度数为0的节点，放到序列的最前面或最后面，不参与BFS。固定放在开头还是末尾取决于求解器偏好，关键是它们的相对顺序无所谓，但必须集中存放，不能让一个孤立节点夹在两个密集节点中间。

5.5 邻居排序写法不严谨：结果随STL版本漂移

现象：同一份代码在一台机器上排序正常，换个编译器跑出来的置换数组和原来不一样。

原因：std::sort是不稳定排序，如果Lambda比较函数对两个等价节点返回false（即胜负未分），排序后的顺序依赖实现。虽然对本问题结果影响可能不大，但在调试时会造成“换了环境结果就变”的假象。

解决：比较函数里加完整的平局策略。我的习惯是度数优先、编号升序打平，让每个比较都有明确结论。排查时如果怀疑排序不稳定，可以直接在Lambda函数里输出比较对，确认每个a和b都有严格的大小关系。

6. 让RCM进入常规流程：多分量拼接、与嵌套剖分配合和性能热点检查

真正把RCM用顺之后，我会做三件额外的事。第一是写一个统一的重排入口，把连通分量检测、伪端节点查找、BFS、孤立节点处理全部封装进去，调用方只需要传CSR数组和两个输出指针。这样排序逻辑只维护一处，不会出现某个调用点漏掉孤立节点处理的低级问题。

第二是把RCM和嵌套剖分（Nested Dissection）组合用。RCM适合中规模矩阵（几十万自由度以内），更大规模时单靠RCM无法充分压制二维三维模型的填充元。常见做法是先做RCM重排，把矩阵压缩到更集中的带形，然后做一次粗粒度划分，对划分出的子块递归剖分。RCM在这里的作用是让剖分前矩阵更“紧凑”，降低粗化阶段的负载不均。

第三是每次排序后打印一个热点检查报告：带宽、轮廓、每行最左侧非零元的位置分布。分布如果不均匀，比如前半段密集后半段稀疏，说明起点选择有问题或者分量拼接顺序不合理，这比盯着带宽一个数字更能定位问题。

我从那以后每次接入新的稀疏求解器，都强制走一遍“对称化—RCM重排—带宽轮廓对比”这三个步骤，确认数字有数量级级别的改善才继续后续开发。很多玄学级别的求解性能问题，最后查出来就是少了这一道重排。希望帮到你。

本文还有配套的精品资源，点击获取