简介：经典代数多重网格（AMG）演示代码，面向数值计算、科学工程计算与偏微分方程数值解学习者，帮助理解如何以代数方式构造粗网格，加速稀疏线性方程组求解。资源为系统开源MATLAB实现，压缩包共20个文件，其中11个m脚本对应AMG核心模块（V循环、强连接依赖矩阵、粗网格生成、粗问题构造、光滑迭代与插值），4个mat文件存放泊松方程有限元离散后的系数矩阵与右端项，2个png图片展示FEM网格及AMG与GMG的对比，另有README说明，整体仅309KB，结构清晰便于按模块阅读。已有1227人学习。代码刻意保持简洁易读，可直接调整平滑器前后次数等参数开展实验；附带五点/九点差分模板生成器以及均匀三角形单元FEM算例，并与几何多重网格（GMG）进行对比，可帮助读者体会代数粗化与几何粗化的差异。生成粗网格环节可能稍慢，且未包含Classic AMG的第二遍粗化过程，适合侧重理解AMG主流程的入门者与教学场景，并附Yousef Saad等参考文献线索便于深入探索。

1. 为什么 MATLAB 里还要自己写 AMG：当内置求解器开始吃内存时

经典代数多重网格（AMG）这个 Demo 解决的是一个很具体的工程痛感：当稀疏线性方程组规模跨过某个量级，直接分解法开始把内存翻倍吃掉，而 MATLAB 自带的 pcg 迭代法在无预条件时收敛慢到让人怀疑人生。Classic_AMG_Demo 用一小段可读性极好的 MATLAB 代码，把 AMG 的两层网格构建、插值算子、粗网格修正和一整套 V 循环完整演示了一遍。它不追求硬件极限，而是把求解思路拆开给你看——适合刚接触多重网格、想在理解公式的同时看到向量化实现的工程师和学生。这套代码是开源的，下载后直接改参数就能跑，是快速建立对 AMG 直觉的最短路径。

2. 经典 AMG 的算法骨架：从矩阵到分层网格的四步流程

AMG 的核心思路是：几何多重网格靠网格分层处理高低频误差，而代数多重网格只靠系数矩阵本身来分层，所以任何稀疏矩阵哪怕没有网格背景都能套用。Classic_AMG_Demo 实现的是经典 Ruge-Stüben 路线，完整流程可以拆成四个步骤：光滑器、强耦合矩阵、C/F 分裂、插值算子。这四步的顺序在代码里不能乱，因为每一步的产物都是下一步的输入。下面我用一个小例子来解释每一步为什么存在，以及参数调错时你会看到什么。

2.1 光滑器为什么选加权 Jacobi：收敛速度与代码复杂度的平衡

多重网格的第一步是光滑（smoothing）。它的目的不是直接求解，而是把误差里的高频分量压下去，留给粗网格去处理低频误差。Gauss-Seidel 收敛速度通常更好，但 MATLAB 里逐点更新的循环写起来绕，而加权 Jacobi 可以直接矩阵向量化。

% 加权 Jacobi 光滑器，w 取 0.67 到 0.8 之间 function x = smooth_wjacobi(A, x, b, nu, w) d = diag(A); for k = 1:nu x = x + w * (b - A * x) ./ d; end end

这段代码的核心是b - A*x残差计算，./ d相当于乘上对角阵的逆，w 是阻尼系数。为什么取 0.67 而不是 1？因为纯 Jacobi 迭代（w=1）虽然实现最简单，却经常不收敛，尤其在对角占优但不绝对占优的五点差分矩阵上，谱半径可能接近 1，迭代上百步高频误差依然压不掉。加权以后，光滑器的谱半径被压缩到 0.8 以下，稳定性好很多。

参数上，nu 表示每个 V-cycle 前后各做几次光滑。Demo 里默认给 [2,2]，意思是进入粗网格前做 2 次，回来后再做 2 次。如果你把 nu 调成 [1,1]，收敛会变慢；调成 [4,4]，每轮代价上升但未必换来回合数的等比下降，因为低频误差在细网格上本来就不是光滑器的职责。还有一个容易被忽略的点：diag(A)在进入循环前取一次就够了，别写在 for 里面，否则大规模稀疏矩阵每次都要重取对角元，白浪费一截时间。

对比 Gauss-Seidel，它的收敛因子通常比加权 Jacobi 小，但每步需要按行更新，难以向量化。在 MATLAB 里通常用A\\b或mldivide，可那又回到了直接法。所以演示代码选 Jacobi 是刻意为之：把概念说清楚，同时让代码能在中等规模矩阵上跑出可观察的收敛曲线。真实项目里，如果矩阵是好的 M 矩阵，你也可以把光滑器换成 Gauss-Seidel 或 SSOR，这会在第 6 章讲复用技巧时再展开。

2.2 C/F 分裂与强耦合阈值：决定粗网格质量的第一道关口

粗网格选择是整个 AMG 的黑匣子部分——选哪些点留在粗网格（C 点）、哪些点只能靠插值（F 点），直接决定收敛速度。

要选 C 点，首先得知道谁和谁强耦合。经典做法是计算强耦合矩阵 S，规则是：当 |a_ij| 大于等于当前行最大非对角元的 theta 倍时，认为 j 是 i 的强连接。

% 构建强耦合矩阵 S，theta 为阈值 function S = strong_connections(A, theta) n = size(A, 1); [ii, jj, aa] = find(A); % 拿到全部非零元 max_row = full(max(abs(A), [], 2)); threshold = theta * max_row(ii); % 每个非零元对应的行阈值 keep = abs(aa) >= threshold; S = sparse(ii(keep), jj(keep), 1, n, n); S(1:n+1:end) = 0; % 去掉对角，避免自耦合 end

这里先取所有非零元的行列索引，再按行阈值过滤，最后用 sparse 一次性拼出 S。比逐行判断快得多，而且不会因为稀疏矩阵被转成稠密而爆内存。theta 的取值在经典 AMG 里通常取 0.25，但你要意识到这不是什么物理常数，而是一个经验阈值。取太小，强连接集合过大，粗网格上点太多，稀疏结构被撑大；取太大，强连接集合过于稀疏，F 点和 C 点之间失去代表关系，插值做出来不收敛。

C/F 分裂的过程用一句话说就是：从强连接度最大的未标记点开始，把它标记为 C，然后把它的强连接邻居全部标记为 F；重复这个过程，直到所有点都被染色。

% 贪心 C/F 分裂：一次遍历，适合 Demo 理解 deg = full(sum(S, 2)); [~, order] = sort(deg, 'descend'); CF = zeros(N, 1); for j = 1:length(order) idx = order(j); if CF(idx) == 0 CF(idx) = 1; % 标记为 C neighbors = find(S(:, idx)); CF(neighbors(CF(neighbors) == 0)) = -1; % 强邻居标为 F end end

这里CF数组 1 代表 C，-1 代表 F。关键动作是CF(neighbors(CF(neighbors) == 0)) = -1，它只把尚未标记的强邻居标成 F，已经标过的保持不动。这个朴素实现的好处是一遍跑完，代码好懂；缺点是分裂结果不是最优的，某些 F 点可能缺少足够的 C 邻居支撑，导致插值矩阵某一行几乎没非零元。碰到这种情况时，收敛曲线会出现很诡异的平台-跌落交替，后面第 5 章会专门讲怎么排查。

2.3 插值算子构建：经典 Ruge-Stüben 插值的核心公式

插值算子 P 的意义是：把粗网格上的误差修正量映射回细网格。限制算子 R 反过来，把细网格残差压到粗网格，通常取 R = P'。

经典插值公式基于矩阵行方程推导。对每个 F 点 i，设 C_i 是 i 强连接的 C 点集合，D_i 是强连接的 F 点集合，弱连接部分被吸收进对角元。用直接插值版本表达就是：

% 经典直接插值：F 点通过其强连接的 C 邻居取值 alpha = A(i, i) + sum(A(i, D_i)); P(i, C_i) = -A(i, C_i) / alpha;

这里的数学含义是：把 i 点方程里所有非 C 强邻居的贡献都折算到对角上，得到修正后的对角 alpha，然后用这个 alpha 归一化对 C 邻居的耦合强度。alpha 越小，插值权重越大，表示该 F 点越依赖 C 点信息，这是符合直觉的。完整的 Ruge-Stüben 插值还会把 D_i 内 F 点对 C 点的间接贡献再累加进去，公式上更严谨，但演示代码采用直接插值版本是够用的——至少在均匀网格和高斯型矩阵上，两者收敛步数差别不大。

向量化实现时，逐行循环当然可以，但如果矩阵是 1e5 阶，MATLAB 的 for 循环会明显拖时间。常见做法是先确认所有 C 点的全局编号，再为每个 F 点建立稀疏行，代码示例放在下一章，因为搭建过程容易出错，单独拿出来讲更有条理。

关于 R = P' 为什么成立，粗网格算子 A_c = RAP，如果 R 不是 P'，那么 A_c 不保证对称，后面的 CG 类外层迭代就不能用。演示代码里把这个关系做成硬约束，是值得保持的好习惯。

2.4 V-cycle 与算子组合：粗网格修正的完整路径

有了 P 和 R，粗网格算子 A_c = R * A * P。这一步在 MATLAB 里就是稀疏矩阵乘法，但要注意稀疏矩阵乘完以后会变得更稠密，这是正常的，只是提醒你别在循环里反复做这个乘法。

function x = vcycle(A, P, R, Ac, x, b, nu) x = smooth(A, x, b, nu(1)); r = b - A * x; rc = R * r; if size(Ac, 1) < 50 ec = Ac \ rc; else ec = vcycle(Ac, Pc, Rc, Acc, zeros(size(rc)), rc, nu); end x = x + P * ec; x = smooth(A, x, b, nu(2)); end

注意递归终止条件是粗网格小到可以直接分解。这个阈值取多少，取决于你的内存和耐心。我一般会把它设成 100 以下，因为粗网格矩阵即使只有 50 阶，直接分解也只是一瞬间的事。还有一个细节是递归调用时必须传入该层对应的 Pc、Rc、Ac，而不是复用高层变量；Demo 代码里如果只建了两层，那第一层递归就直接落到底，少一层也无所谓。

V-cycle 名字的由来就是从细到粗、再回到细的路径像一个 V 字。更激进的 W-cycle 会在粗层多做一次循环，收敛因子更小，但代价也更高。演示代码里只用 V-cycle，已经足够展示 AMG 的核心行为。

如果你把 V-cycle 的收敛曲线画出来，会看到一条特别干净的直线，从第一个循环开始就以稳定速率下降。如果前几步下降很快、后面突然平了，问题几乎一定出在插值算子对低频误差的表示能力不足，也就是 C/F 分裂里 C 点太少了。

提示：调参时只记录最终迭代数不够，把整条残差曲线存下来，结合前三步斜率一起判断，才能定位是光滑不足还是粗网格修正不足。

3. Demo 代码的关键实现：逐段拆解 Classic_AMG_Demo.m

这一章我们把 Classic_AMG_Demo.m 按执行顺序拆开。实际上源码本身规模不大，主体大概几个函数：生成测试矩阵、强耦合矩阵、C/F 分裂、构建插值、V-cycle、主脚本。理解顺序比理解单行更关键。

3.1 测试矩阵怎么生成：五点差分离散化与稀疏矩阵构造

AMG 最经典的测试床是二维泊松方程 -Δu = f 在单位正方形上经过五点差分离散化得到的矩阵。在 MATLAB 里可以直接调用 gallery 函数，也可以手动拼：

% 生成二维 Poisson 五点差分矩阵，每维 n 个内部网格点 n = 64; e = ones(n, 1); A1 = spdiags([-e 4*e -e], [-1 0 1], n, n); % 一维三对角 I = speye(n); A = kron(I, A1) + kron(A1, I); % 二维五点差分 N = size(A, 1); % N = n^2

这个 Kronecker 积构造法的原理是：一维 n 点 Laplacian 的三对角矩阵记为 A1，二维网格的拉普拉斯可以写成 A = I ⊗ A1 + A1 ⊗ I，其中第一个 kron 处理 y 方向，第二个 kron 处理 x 方向。这样生成的矩阵只包含中心点和上下左右五个非零元，对角线是 4，四条副对角线是 -1，完全符合五点差分格式。

这里有一个容易踩的小坑：如果你在spdiags里把-e 4*e -e的顺序写反，得到的矩阵符号反了，对角线变成 -4，AMG 插值公式全部失效。所以生成之后立刻用max(abs(A*A - A*A))检查一下对称性，或者用spy(A)看一下非零元分布是否符合五点差分形状。

另外，右端项 b 的建议取 b = A * ones(N,1)，初始解 x0 取 zeros(N,1)。这样设置的精确解就是全 1 向量，V-cycle 每步迭代都能算真实误差，而不仅是残差。残差小不等于误差小，在多重网格里这两者可以差好几个数量级。

3.2 粗网格选择实现：从强耦合矩阵到 C/F 标记

演示代码里 C/F 分裂通常用一个两步标记完成。前面 2.2 节给的是一次遍历版本，这里再补充一个更接近工程实现的二次修正版：

% 一次遍历后统计孤立 F 点 isolated = find(CF == -1 & sum(S .* (CF == 1)') == 0); for i = isolated' CF(i) = 1; % 强制升级为 C 点 end

这段代码的作用是修复一种经典翻车现场：贪心分裂后有个别 F 点周围全是 F，没有 C 点可以直接取值的对象，插值矩阵那行会变成全零。强制升级为 C 点后，粗网格多一个点，但换来了插值矩阵每行至少有一个非零元，收敛曲线会立刻恢复正常。

为什么要强调这一步？因为很多从教材拷出来的演示代码会忽略它，而真实工程矩阵和规则网格不同，孤立 F 点出现的概率并不低。你在跑 Demo 时可以故意注释掉这段二次修正，对比一下收敛历史，就能直观体会每行至少一个 C 邻居对 AMG 有多重要。

% 一次遍历后统计孤立 F 点 isolated = find(CF == -1 & sum(S .* (CF == 1)') == 0);

这里的sum(S .* (CF == 1)')是在算每个 F 点的强邻居中有多少个 C 点，如果结果是 0，说明它孤立了。这个矩阵乘法的写法比逐行 find 快，但要注意CF == 1是逻辑向量，转置后和 S 做逐元素乘法，会产生一个稀疏矩阵，sum 之后变成向量。

3.3 插值与限制算子的 MATLAB 向量化实现

构建插值算子的逐行循环写法，在小网格上没问题，网格一大就让人着急。向量化的一种方式是先确定所有 F 点的 C 邻居索引，然后把 A 的非零元分散到 P 的对应位置：

c_list = find(CF == 1); NC = length(c_list); P = sparse(N, NC); for i = find(CF == -1)' Ci = find(CF == 1 & logical(S(:, i))); % i 强连接的 C 点 Di = setdiff(find(CF == -1 & logical(S(:, i))), i); % 强连接的 F 点，排除自身 alpha = A(i, i) + sum(A(i, Di)); P(i, Ci) = -A(i, Ci) / alpha; end R = P'; Ac = R * A * P;

这段代码在演示级别已经足够清晰。要注意Ci和Di分别对应公式里的 C_i 与 D_i，顺序不能反。Di里显式排除自身是必须的，否则 alpha 会把对角线自己算进去一次，导致权重偏小、插值偏保守，收敛速度下降。

P 初始化为 sparse(N, NC)，NC 是 C 点总数。构建完成后nnz(P)应该约等于 F 点数的 2～4 倍，如果某个行完全没有非零元，说明刚才的孤立 F 点修复没起作用。R 和 P 的转置关系保证了粗网格算子的对称性，这是后续 CG 类算法能用的前提。Ac 的稠密度通常比 A 高 20% 左右，属于正常现象；如果高出一倍以上，就得回到 theta 和 C/F 分裂去查。

注意：构建 P 的时候顺手加一句assert(size(P,2) == length(c_list))，能在索引错位时第一时间报错，而不是等整个 V-cycle 跑飞了才回头找。

3.4 主循环：V-cycle 迭代与收敛历史记录

主脚本的最后一步是驱动 V-cycle 迭代若干次，并记录残差范数。注意不要把残差和真实误差混了，有精确解时两个都算最好。

x = zeros(N, 1); res(1) = norm(b - A * x); maxit = 20; tol = 1e-8; for k = 1:maxit x = vcycle(A, P, R, Ac, x, b, [2 2]); res(k+1) = norm(b - A * x); if res(k+1) < tol * res(1) break; end end semilogy(0:k, res);

semilogy 画出来的收敛曲线，如果是一条斜率为常数的直线，说明 AMG 的收敛因子稳定。如果曲线在后期变平甚至上翘，说明插值或 C/F 分裂有问题。还有一个判断技巧：把 y 轴范围固定住，比如 [1e-12, 1]，这样不同参数下的曲线才能在同一起跑线上对比，否则每次自动缩放会掩盖掉平台期的细节。

主循环里还有一个细节是 maxit 别设太大，我一般给 20 就够。AMG 的收敛因子通常在 0.1 左右，20 个循环足够把残差压 12 个数量级。如果你的残差 20 步都压不下去，去调参数比加大迭代次数有意义得多。

4. 参数怎么调、结果怎么读：数值实验与对比维度

Demo 代码不是拿来当黑盒跑一遍就完事的。下载回来第一件事是改参数、看收敛历史。这一章给出一组可复现的实验方法，帮你判断自己的矩阵适不适合这套经典 AMG。

4.1 网格规模与迭代次数的关系：从 32x32 到 128x128

理论上，AMG 的收敛速度应该与网格规模无关，这是它比 CG 加对角预条件强的地方。你可以这样验证：

for n = [32, 64, 128] A = gallery('poisson', n); [P, R, Ac, CF] = setup_amg(A, 0.25); x = zeros(size(A,1), 1); b = A * ones(size(A,1), 1); [iter, res] = run_vcycle(A, P, R, Ac, x, b, 20); fprintf('n=%4d, iter=%3d, final_res=%.2e\n', n, iter, res(end)); end

下面是我在经典矩阵上看到的一组合格参考值：

实际项目中如果迭代次数随 n 明显增加，先别怀疑 AMG 不行，多半是 C/F 分裂阈值取得不对，或者插值权重漏掉了弱连接吸收项。测试时固定 theta=0.25、w=0.67、nu=[2,2] 这三组基准参数，然后再逐个调。

这个实验的另一个价值是帮你看 setup 时间。AMG 的 setup（强耦合矩阵、C/F 分裂、P 和 Ac 的构建）只做一次，后续每个 V-cycle 代价都很低。n=128 时 setup 可能占掉总时间的 40%，n=512 时这个比例往往降到 20% 以下，说明规模越大，AMG 摊销 setup 成本的能力越强。

4.2 光滑器权重 w 的影响：什么值最容易翻车

w 太低，比如 0.5，光滑效果不足，高频误差没压下去，V-cycle 整体收敛变慢。w 太高，比如 0.9，在 M 矩阵上容易导致 Jacobi 迭代本身振荡，残差曲线会出现锯齿。经典理论推荐区间 [0.6, 0.8]，实际我用 0.67 起步，然后再试 0.7、0.75。

这里有个容易被忽略的细节：加权 Jacobi 里的 w 是乘在修正量上的，不是乘在松弛后的值上。写错位置的话，w=1 变成纯 Jacobi，w>1 直接发散。在 Demo 中检查一下x = x + w * (b - A*x) ./ d这行的括号位置。

做参数扫描时有个小技巧：不要只记录最终迭代次数，还要记录收敛曲线的前三步斜率。三步斜率能反映光滑器高频误差压制能力，最终迭代次数则反映插值与粗网格的综合表现。两个指标一起看，才能区分是光滑不够还是粗网格修正不够：

这张表的规律是：w 太大时前三步可能看起来不错，但后续会现出原形。所以只盯着最终迭代数是不够的，要把整条曲线保存下来看形状。

4.3 与 MATLAB 自带 pcg/ilu 的对比：什么场景下 AMG 值得用

对二维泊松这种理想问题，直接用 pcg 加 incomplete Cholesky 预条件，速度可能比 AMG 还快，因为问题太规整。AMG 的价值在缺口更大的地方：电阻网络、非结构网格、含强各向异性系数的扩散问题。

对比实验建议这么做：同一矩阵，分别用pcg(A, b, tol, maxit)、pcg(A, b, tol, maxit, ichol(A))和自己封装的 AMG 循环各跑一遍，记录时间和迭代数。你会发现 AMG 的 setup cost 比较高，但每步迭代代价低，所以当迭代次数超过 20 时，AMG 优势就开始显现。如果矩阵是三维问题，这个临界点会更早出现。

% 对比三种求解方式，注意用 tic/toc 分别记录 setup 和迭代时间 tol = 1e-8; maxit = 100; tic; x1 = pcg(A, b, tol, maxit); t_pcg = toc; tic; x2 = pcg(A, b, tol, maxit, ichol(A)); t_ichol = toc; tic; [x3, iter] = run_amg(A, b, tol); t_amg = toc; fprintf('pcg(%.2fs) ichol+pcg(%.2fs) amg(%.2fs, %d iter)\n', ... t_pcg, t_ichol, t_amg, iter);

跑这个对比时要小心一件事：ichol(A)对非正定矩阵会直接报错，而经典 AMG 只要矩阵非奇异都能跑。所以如果矩阵是从对流扩散方程来的非对称矩阵，pcg/ichol 这条线可能根本不适用，反而是 AMG 的天然主场。这也解释了为什么在 CFD 和电路仿真领域，AMG 的出场率远比通用科学计算里高。

5. 避坑与常见问题排查：动手跑 Demo 时最容易踩的五个坑

我前前后后把这个 Demo 在七八种矩阵上反复跑过，下面这五个坑是出现频率最高的。每一条都按现象、原因、解决三步说透。AMG 的收敛性有时候挺玄学的，但你遇到的绝大多数问题其实都出在代码逻辑上，而不是算法本身。

5.1 现象：粗网格矩阵越来越稠密，甚至比细网格还大

现象：nnz(Ac)比nnz(A)还高，内存直接翻倍。

原因：C/F 分裂选的 C 点太少，导致 P 每行向粗网格投射大量非零元，RAP 乘出来自然稠密。根本原因又回到 theta 取值上，theta 太小则强连接集合过大，贪心分裂容易选出很多 F 点，C 点比例太低。

解决：把 theta 调大，比如从 0.25 调到 0.3 或 0.35。theta 越大，强连接越少，C 点越多，粗网格越小。注意不要超过 0.5，否则强连接集合变得太稀疏，插值失去意义。调完以后用length(find(CF==1))/N看一下 C 点比例，经验上这个比例应该在 0.3～0.5 之间，偏离太远就该回来改 theta。

5.2 现象：插值算子拼接出错导致解完全发散

现象：残差第一次迭代就大到 1e10，semilogy 图上直接冲出边界。

原因：P 的列索引和 C 点在粗网格中的编号不是一一对应。演示代码通常用find(CF==1)生成 C 点列表，但给 P 赋值时忘了建立"原始节点号 → 粗网格列号"的映射。如果直接拿全局节点号当列号用，P 的列数会大于 NC，稀疏矩阵的赋值会把多余的列自动扩展出来，而 Ac = RAP 的维度又恰好对得上，错误被隐藏到最后一步才爆发。

解决：显式存一个c_list = find(CF == 1)，并让P(:, j)的 j 对应到c_list(j)，而不是原始的全局节点编号。构建 P 的时候可以顺手加一句断言assert(size(P,2) == length(c_list))，防的就是这种索引错位。

5.3 现象：V-cycle 不收敛，残差曲线平坦

现象：残差在 1e-2 附近挂着，怎么迭代都下不去，曲线像一条水平线。

原因：多半是插值里的贡献项都来自 F 点而不是 C 点，或者 D_i 集合里把 i 自己算进去了，造成 alpha 被多减了一块。另一个可能是 R 和 P 不是转置关系，导致粗网格算子不对称，V-cycle 的收敛理论全部失效。

解决：先用spy(P)看插值矩阵的非零结构，确认每行至少有一个非零元；再打印 P 的前 10 行，对照 CF 数组检查取值来源是否都是 C 点。最后验证 R 和 P 的转置关系：assert(isequal(R, P'))。这几步走完，90% 的平坦曲线都能救回来。

5.4 现象：MATLAB 内存不足，稀疏矩阵被意外转为稠密

现象：运行到Ac = R * A * P时 MATLAB 直接报内存错误，或者在主循环里突然卡住十几秒。

原因：可能是 R 或 P 中间变量被存成了全矩阵。常见误操作是把diag(A)的结果直接赋值给d，如果 A 是 sparse，diag(A)返回的是稀疏列向量，没问题；但一旦你在构造 S 时用了类似full(S)的调试语句又忘了删，后面所有数组都被带成稠密，内存就爆了。

解决：养成习惯，所有与 A 同规模的向量都用 sparse 或按稀疏操作处理；检查变量后加whos确认类型是 sparse 而不是全矩阵。还有一个排查技巧：在赋值语句周围写disp(nnz(P))，如果 nnz 突然从几千跳到几十万，那一定有什么东西被转成稠密了。

5.5 现象：网格规模改变后收敛性骤降

现象：n=32 时 7 步收敛，n=256 时十几步还停不下来，甚至出现不收敛。

原因：C/F 分裂的贪心顺序对点数敏感，特别是阈值边界附近的点稍一变化，分裂结果就差异很大。这不是代码 bug，而是经典 AMG 的已知敏感点。规则网格相对还好，一旦矩阵来自非结构网格，这个现象会放大很多倍。

解决：改用更鲁棒的粗网格选择策略，比如把度数排序改成"先按强连接度再按对角元大小"的加权排序；或者在固定网格上先用参数扫描确定稳定区间，再放大网格。还有一招是从 V-cycle 换成 W-cycle，粗层多做一次循环，可以在不改动 C/F 分裂的前提下把收敛因子拉回来。

6. 把 Demo 改造成自己的求解器：从演示代码到可复用函数的三个技巧

演示代码写得好，但直接拿进工程里还是太单薄。这三个技巧能让你的 AMG 从"能跑"变成"能用在其他问题上"。

第一个技巧是参数化光滑器。不要写死某一种迭代格式，而是把光滑器当函数句柄传进 V-cycle：

smoother = @(A,x,b,nu) smooth_wjacobi(A,x,b,nu,0.67); x = vcycle(A, P, R, Ac, x, b, nu, smoother);

换一种光滑器只是换一行函数句柄，改造成本趋近于零。第二个技巧是用 spdiags 写一个接受系数函数句柄的测试矩阵生成器，这样你能快速生成各向异性扩散矩阵、变系数矩阵来测试 AMG 的鲁棒性，而不是每次手拼 kron。第三个技巧是养成跑 profile 的习惯：profile on; run_demo; profile viewer能直接看到哪一行偷偷把稀疏矩阵转稠密，往往这类隐藏转换比算法本身还吃内存。

从那以后，我每次拿到别人的 AMG 代码，第一件事就是跑一遍 profile，把耗时前五行指给他们看，再聊算法。先跑通、再优化、最后再谈收敛性，这个顺序能帮你省掉大量无效调参时间。希望这个习惯也帮到你。

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