简介经典代数多重网格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 迭代w1虽然实现最简单却经常不收敛尤其在对角占优但不绝对占优的五点差分矩阵上谱半径可能接近 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 的强连接。% 构建强耦合矩阵 Stheta 为阈值 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:n1: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的顺序写反得到的矩阵符号反了对角线变成 -4AMG 插值公式全部失效。所以生成之后立刻用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 点数的 24 倍如果某个行完全没有非零元说明刚才的孤立 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(k1) norm(b - A * x); if res(k1) 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矩阵阶数V-cycle 达到 tol1e-8 的迭代次数321024约 7644096约 712816384约 8实际项目中如果迭代次数随 n 明显增加先别怀疑 AMG 不行多半是 C/F 分裂阈值取得不对或者插值权重漏掉了弱连接吸收项。测试时固定 theta0.25、w0.67、nu[2,2] 这三组基准参数然后再逐个调。这个实验的另一个价值是帮你看 setup 时间。AMG 的 setup强耦合矩阵、C/F 分裂、P 和 Ac 的构建只做一次后续每个 V-cycle 代价都很低。n128 时 setup 可能占掉总时间的 40%n512 时这个比例往往降到 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 是乘在修正量上的不是乘在松弛后的值上。写错位置的话w1 变成纯 Jacobiw1 直接发散。在 Demo 中检查一下x x w * (b - A*x) ./ d这行的括号位置。做参数扫描时有个小技巧不要只记录最终迭代次数还要记录收敛曲线的前三步斜率。三步斜率能反映光滑器高频误差压制能力最终迭代次数则反映插值与粗网格的综合表现。两个指标一起看才能区分是光滑不够还是粗网格修正不够w前三步残差比总迭代次数tol1e-8现象0.50.612收敛慢但稳定0.670.457推荐基准0.80.46收敛最快但接近振荡边界0.950.659前三步锯齿后期开始振荡这张表的规律是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) icholpcg(%.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(CF1))/N看一下 C 点比例经验上这个比例应该在 0.30.5 之间偏离太远就该回来改 theta。5.2 现象插值算子拼接出错导致解完全发散现象残差第一次迭代就大到 1e10semilogy 图上直接冲出边界。原因P 的列索引和 C 点在粗网格中的编号不是一一对应。演示代码通常用find(CF1)生成 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 是 sparsediag(A)返回的是稀疏列向量没问题但一旦你在构造 S 时用了类似full(S)的调试语句又忘了删后面所有数组都被带成稠密内存就爆了。解决养成习惯所有与 A 同规模的向量都用 sparse 或按稀疏操作处理检查变量后加whos确认类型是 sparse 而不是全矩阵。还有一个排查技巧在赋值语句周围写disp(nnz(P))如果 nnz 突然从几千跳到几十万那一定有什么东西被转成稠密了。5.5 现象网格规模改变后收敛性骤降现象n32 时 7 步收敛n256 时十几步还停不下来甚至出现不收敛。原因C/F 分裂的贪心顺序对点数敏感特别是阈值边界附近的点稍一变化分裂结果就差异很大。这不是代码 bug而是经典 AMG 的已知敏感点。规则网格相对还好一旦矩阵来自非结构网格这个现象会放大很多倍。解决改用更鲁棒的粗网格选择策略比如把度数排序改成先按强连接度再按对角元大小的加权排序或者在固定网格上先用参数扫描确定稳定区间再放大网格。还有一招是从 V-cycle 换成 W-cycle粗层多做一次循环可以在不改动 C/F 分裂的前提下把收敛因子拉回来。6. 把 Demo 改造成自己的求解器从演示代码到可复用函数的三个技巧演示代码写得好但直接拿进工程里还是太单薄。这三个技巧能让你的 AMG 从能跑变成能用在其他问题上。第一个技巧是参数化光滑器。不要写死某一种迭代格式而是把光滑器当函数句柄传进 V-cyclesmoother (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把耗时前五行指给他们看再聊算法。先跑通、再优化、最后再谈收敛性这个顺序能帮你省掉大量无效调参时间。希望这个习惯也帮到你。本文还有配套的精品资源点击获取