解线性方程组大概是线性代数里最实用的一块内容而高斯-约当Gauss-Jordan消元法是我个人在写数值计算代码时最常用、也最“省脑”的方法。它的步骤非常规整构造增广矩阵、选主元、归一化、消掉其他行最后左边变成单位矩阵右边就是解。相比经典高斯消元还要回代Gauss-Jordan把整个求解过程压缩成一套统一循环特别适合写成通用函数。这篇博文不讲花架子直接给出一份能跑的C实现同时把每一步设计背后的原因、容易踩的坑以及工程上什么时候该自己写、什么时候该用库都尽量说清楚。适合正在学数值方法、准备算法竞赛或者需要在嵌入式、无依赖环境里解方程的同学参考。1. 先搞懂算法本身它和高斯消元的真正区别1.1 从增广矩阵说起线性方程组可以写成 Ax b 的形式。A 是系数矩阵x 是未知数向量b 是常数向量。高斯-约当消元法的思路是把 A 和 b 拼成一个增广矩阵也就是在 A 右边直接多加一列常数项。比如x 2y z 82x y z 7x y 3z 12增广矩阵就是[ 1 2 1 | 8 ][ 2 1 1 | 7 ][ 1 1 3 | 12 ]中间的竖线只是给人看的程序里不需要单独存这一列矩阵宽度直接加 1 就行。解方程组的过程本质上是对这个矩阵做三种行变换交换两行、某行整体乘以非零常数、某行加上另一行的倍数。这三种操作都是“行等价变换”不会改变方程组的解集所以随便怎么折腾最后得到的矩阵和原方程是同一个解集的不同写法。我们的目标是把系数部分化成行最简形Reduced Row Echelon Form简称 RREF每个主元所在的列除了主元本身是 1 之外其他位置全是 0。一旦化成这种形态方程组的解相当于直接印在最后一列上。1.2 为什么我常用高斯-约当而不是经典高斯消元经典高斯消元是两步走先通过消元把系数矩阵化成上三角然后再从最后一个方程开始往上回代。高斯-约当则是“一步到位”在消元过程中不仅消掉主元下方的元素还把主元上方的元素也全消掉最后系数矩阵直接变成单位矩阵。为了方便对比我列个表对比项经典高斯消元高斯-约当消元最终形态上三角矩阵行最简形RREF求解方式回代直接读解乘加次数约 n³/3约 n³/2代码复杂度和边界情况回代部分容易写错统一循环逻辑简单求逆矩阵扩展需额外处理右侧直接在右侧拼单位阵计算量上高斯-约当确实比高斯消元多一点大概多 50% 的乘加操作。但很多场景里我们不只求一个方程组的解而是要求多个右端项甚至直接求矩阵的逆。这个时候高斯-约当有个天然优势把 A 和单位矩阵 I 拼成增广矩阵跑完一遍之后右侧就是 A⁻¹。一步到位没有额外代码。而且从工程角度看代码少一个回代环节就少一类维护边界条件的麻烦。对于 n 在几千以内的中小规模问题这点性能差距几乎可以忽略。可维护性和逻辑统一性反而更重要。1.3 无解、唯一解、无穷解怎么一眼看出行最简形化完以后解的情况可以通过三条规则直接判断无解某一行出现了“左侧全零、右侧非零”相当于数学上的 0 非零数矛盾。无穷解主元个数小于未知数个数说明存在自由变量。唯一解主元个数等于未知数个数每个变量都被唯一确定解就是最后一列的数字。例如某个方程组行最简形是[ 1 0 | 2 ][ 0 1 | 3 ]主元数是 2未知数也是 2唯一解 x2, y3。如果出现[ 1 0 0 | 1 ][ 0 0 0 | 2 ]第二行左侧全零、右侧是 2这就是无解。把这三条写进代码并不复杂真正复杂的是怎么在浮点数的世界里判断“全零”。这就要聊到 EPS 和主元策略了。2. C实现前必须想清楚的三个设计问题2.1 矩阵用什么存二维数组还是嵌套 vector我建议优先用vectorvectordouble不要自己管理二维裸数组。原因很实际行交换方便。C 的swap(mat[row1], mat[row2])可以直接交换两个 vector底层只是交换指针开销极小。矩阵行数、列数可以动态确定不用在编译期写死函数适配性更好。打印和调试直观循环起来就是一个标准的二维结构。有同学会担心性能。说实话几千阶以内的矩阵vectorvectordouble完全够用。如果你确实要处理超大矩阵又不想引入第三方库可以用一维 vector 自行模拟二维索引mat[i * cols j]。这样内存连续缓存友好度更高。不过付出的代价是代码可读性下降行交换需要手动实现内存拷贝。另外提醒一个小坑代码里涉及到矩阵行列数时最好强制转成int不要用size_t直接参与比较。否则for (int i 0; i mat.size() - 1; i)这种写法在mat为空时会爆炸因为size_t无符号减出一个巨大正整数。2.2 主元选择为什么不能偷懒消元的第一步是确定当前列用哪一行作为基准行这一行的主元列元素叫“主元”。最朴素的做法是从当前行往下找第一个非零元素把它当主元。这个思路在理论上没问题但在浮点运算里可能翻车。假设主元是 1e-300归一化时要让这行除以 1e-300结果就是把这个极小的数放大成 1同时它原本携带的相对误差也被放大到天文数字级别。后续消元操作会把误差扩散到其他行最终解出来的结果可能完全不能用。标准解法是“部分主元法”在当前列从当前行往下扫描找绝对值最大的那个元素所在行把它交换到当前行。这个策略的收益只是 O(n) 次比较但数值稳定性提升非常明显。你用 double 写代码如果不做部分主元遇到稍微病态一点的矩阵就等着看诡异结果吧。完整的主元法还要交换列实现复杂度高实际应用中除非做符号计算否则很少用到不推荐自己造轮子。2.3 浮点误差EPS 到底怎么取判断浮点数是不是 0绝不能直接if (x 0.0)。因为经过一系列加减乘除后原本该是 0 的位置很可能是 1e-15 这种微小残差。我的习惯是定义const double EPS 1e-9;所有需要判断为 0 的地方都用fabs(x) EPS。这个值不是拍脑袋定的double 大约有 15 位有效十进制数字1e-9 作为判断阈值比较平衡既能容忍合理误差又不容易把真正该保留的小数误判成 0。如果你的问题数值本身很大比如元素动辄 1e10 级别那 EPS 可以适当放大到 1e-7。反过来如果是严谨的数值实验更推荐long double或者干脆用分数类做精确有理数运算后面第五章会提到。3. 完整代码实现把 Gauss-Jordan 写成通用函数3.1 函数接口怎么设计我习惯用一个整数返回值表达最终状态返回 0唯一解。返回 1无穷多解。返回 -1无解。同时传入的矩阵会被原地修改成行最简形解通过引用参数solution输出。这样设计的好处是调用方可以根据返回值走不同分支同时还能拿到化简后的矩阵方便调试和后续处理。如果只需要判断一次“这方程组有没有解”你可以把返回值当状态码用如果还需要具体解的值就从solution里取。3.2 核心循环逐段拆解整个代码的骨架是一个两层循环外层遍历列内层做“找主元、交换、归一化、消去”四件事。找主元的逻辑从第row行开始往下扫描col列找到绝对值最大的元素记录所在行pivotRow。如果最大绝对值小于EPS说明这列在当前剩余行里全是 0这个变量是自由变量跳过这一列。把pivotRow和row交换让主元跑到当前行。归一化那一步是把第row行从第col列开始的所有元素都除以pivot。为什么从col开始而不是从 0 开始因为col左边的列都已经是格式化完毕的主元列再动它们也不会影响正确性还省了计算量。最后是消去操作。和高斯消元只消下边行不同这里要遍历所有行把主元列的元素全部消成 0。先判断factor的绝对值是否小于EPS如果本身已经接近 0 就跳过避免多余的浮点乘加把误差越滚越大。3.3 完整可运行代码下面是完整代码依赖只有标准库任何支持 C11 的编译器都能直接编。#include iostream #include vector #include cmath #include iomanip using namespace std; const double EPS 1e-9; using Matrix vectorvectordouble; void printMatrix(const Matrix mat, int precision 6) { for (const auto row : mat) { for (double val : row) { cout setw(12) fixed setprecision(precision) val ; } cout endl; } } // 返回值 // 0 - 唯一解 // 1 - 无穷多解 // -1 - 无解 int gaussJordan(Matrix mat, vectordouble solution) { int m (int)mat.size(); if (m 0) return -1; int n (int)mat[0].size() - 1; // 未知数个数 列数 - 1 int row 0; for (int col 0; col n row m; col) { // 1. 部分主元在当前列中找绝对值最大的元素所在行 int pivotRow row; double maxAbs fabs(mat[row][col]); for (int i row 1; i m; i) { if (fabs(mat[i][col]) maxAbs) { maxAbs fabs(mat[i][col]); pivotRow i; } } // 如果这一列从 row 行起全是 0跳过这一列 if (maxAbs EPS) { continue; } // 2. 把主元行交换到当前 row 行 if (pivotRow ! row) { swap(mat[pivotRow], mat[row]); } // 3. 归一化让 mat[row][col] 变成 1 double pivot mat[row][col]; for (int j col; j n; j) { mat[row][j] / pivot; } // 4. 消去其余所有行的 col 列元素 for (int i 0; i m; i) { if (i row) continue; double factor mat[i][col]; if (fabs(factor) EPS) continue; for (int j col; j n; j) { mat[i][j] - factor * mat[row][j]; } } row; } // 检查无解某一行系数全为 0但常数项非 0 for (int i row; i m; i) { bool allZero true; for (int j 0; j n; j) { if (fabs(mat[i][j]) EPS) { allZero false; break; } } if (allZero fabs(mat[i][n]) EPS) { return -1; } } // 主元数量 row 小于未知数数量 n意味着有自由变量 if (row n) { return 1; } // 唯一解解就在最后一列 solution.assign(n, 0.0); for (int i 0; i n; i) { solution[i] mat[i][n]; } return 0; } int main() { cout 示例1唯一解 endl; Matrix a1 { {1, 2, 1, 8}, {2, 1, 1, 7}, {1, 1, 3, 12} }; vectordouble sol; int ret gaussJordan(a1, sol); printMatrix(a1); if (ret 0) { cout 唯一解; for (int i 0; i (int)sol.size(); i) { cout x i 1 sol[i] ; } cout endl; } else if (ret 1) { cout 无穷多解 endl; } else { cout 无解 endl; } cout \n 示例2无解 endl; Matrix a2 { {1, 1, 1}, {1, 1, 3} }; ret gaussJordan(a2, sol); printMatrix(a2); cout 返回码 ret - (ret -1 ? 无解 : 有解) endl; cout \n 示例3无穷解 endl; Matrix a3 { {1, 1, 1, 1}, {2, 2, 2, 2} }; ret gaussJordan(a3, sol); printMatrix(a3); cout 返回码 ret - (ret 1 ? 无穷多解 : 其他) endl; return 0; }这段代码我实测过可以直接跑。几个细节说一下swap(mat[pivotRow], mat[row])交换的是整个 vector不是逐元素交换效率高。归一化和消去都从col开始因为更左边的列已经处理完毕动了也不影响但没必要浪费时间。消去循环里先判断factor绝对值避免把小误差当作有效数去乘加能减少不少无谓运算。4. 测试运行与结果分析验证代码到底靠不靠谱4.1 唯一解案例三个方程三个未知数x 2y z 82x y z 7x y 3z 12理论上解是 x1, y2, z3。代码跑完增广矩阵会变成[ 1 0 0 | 1 ][ 0 1 0 | 2 ][ 0 0 1 | 3 ]打印出来如果看到对角线附近是 1、其他位置是接近 0 的小数比如 1e-15 级别就说明归一化和消去过程正确。浮点运算不可能精确得到 0允许出现微小残差。4.2 无解案例两个方程x y 1x y 3增广矩阵是[ 1 1 | 1 ][ 1 1 | 3 ]代码会先把第一行归一化成 [1, 1, 1]然后用它消第二行得到[ 1 1 | 1 ][ 0 0 | 2 ]第二行对应的方程是 0 2无解。程序走到第 4.1 步检查时发现左侧全零、右侧非零直接返回 -1。这种情况在现实里常见于数据采集有矛盾的过定方程组解出来没有数学意义硬算只会得到一个“最小二乘意义”的近似解。这也是为什么状态码设计要单独区分无解而不是简单给个空 vector。4.3 无穷解案例方程组x y z 12x 2y 2z 2第二个方程其实是第一个的两倍信息完全冗余。代码跑完第二行会被消成全零。主元数量是 1未知数是 3返回无穷多解。如果想要进一步描述这个解集可以从行最简形中读出一个可行解x 1 - y - z其中 y 和 z 是自由变量。代码层面如果确实需要把通解形式也输出那就要额外记录自由变量索引和主元变量索引这部分内容我觉得大多数场景用不上就没有写进函数里。你需要的时候可以在返回1的分支里继续扩展。4.4 多元组测试时的小建议我测试的时候习惯把三组样例写进一个 main用printMatrix看每一步结果。你复制代码运行时如果看到某些位置出现 -0.000000这不是错误是浮点运算里的负零视觉上吓人而已。介意的可以在打印函数里加一行判断如果fabs(val) EPS就输出 0。这块我后来踩过一次小坑明明逻辑没问题看到-0.000000总觉得是不是代码哪里写错浪费了十分钟排查。所以打印函数里把近零数统一成 0是很有必要的调试友好性优化。5. 躲开这些坑精度、效率与工程替代方案5.1 浮点误差的典型症状与对策最常见的症状是明明手算整数解程序解出来却是 0.999999999 或者 1.000000001。这属于 normal 现象不是代码逻辑错误。如果你发现程序结果和手算结果差得离谱优先检查三点有没有做部分主元没做的话遇到主元接近 0 的矩阵很容易爆炸。EPS 设置是否合理太小会导致“接近 0 但被当成非 0”太大会把有效小数误杀。数据本身是否病态比如 Hilbert 矩阵 H[i][j] 1 / (i j 1)n 到 10 左右用 double 做高斯-约当得到的解可能惨不忍睹。这是数值线性代数里非常经典的问题单靠选主元只能缓解不能根治。遇到病态矩阵我的建议是不要死磕自写实现直接上库。5.2 数据规模与性能优化建议这个自写实现的时间复杂度是 O(m × n²)m 是方程数n 是未知数个数。如果 m 和 n 都在几千以内运行时间一般可以接受再大就要认真考虑性能问题了。几个优化方向把vectorvectordouble改成一维数组保证内存连续利用 CPU 缓存。如果只求一组解可以用经典高斯消元加回代比高斯-约当省一点乘加操作。如果同一个系数矩阵要配十几个不同的右端项 b那么 Gauss-Jordan 反而有优势一次消元做完行最简形所有右侧向量同时搞定。开启编译优化-O2对循环优化效果非常明显。5.3 什么时候该换库嵌入式、算法竞赛笔试、教学演示、依赖受限的场景手写实现完全够用。但如果你在公司项目里处理真实数值计算我更推荐直接用现成库EigenC 头文件库PartialPivLU、FullPivLU接口清爽数值稳定性有保障。Armadillo语法接近 MATLAB适合快速原型。LAPACK老牌 Fortran 库dgesv性能极其稳定适合大规模矩阵。这些库内部不只是简单的高斯消元还包含矩阵条件数估计、迭代改进、分块算法等优化。遇到明显病态的问题自己写代码和大厂级的数值库相比稳定性差距不是同一个量级。需要注意的是我说这些不是为了贬低手写实现恰恰相反面试、讲课、理解原理时手写一遍高斯-约当才是最能加深理解的方式。哪怕你以后长期用 Eigen你也得知道FullPivLU到底在干什么才能判断什么时候该用部分主元什么时候该用完整主元。根据我个人写代码和调 bug 的习惯最后再分享一个小改造思路把代码里的double换成自定义分数类实现加减乘除、取绝对值、比较大小这些运算符就可以做精确有理数消元。我拿这套改进版去验证普通 double 版本算出来的结果很多精度问题立刻现出原形。高斯-约当本身不复杂真正体现功力的地方在于你对边界情况的敬畏和对数值误差的理解。 高斯-约当消元法详解：从原理到C++实现
解线性方程组大概是线性代数里最实用的一块内容而高斯-约当Gauss-Jordan消元法是我个人在写数值计算代码时最常用、也最“省脑”的方法。它的步骤非常规整构造增广矩阵、选主元、归一化、消掉其他行最后左边变成单位矩阵右边就是解。相比经典高斯消元还要回代Gauss-Jordan把整个求解过程压缩成一套统一循环特别适合写成通用函数。这篇博文不讲花架子直接给出一份能跑的C实现同时把每一步设计背后的原因、容易踩的坑以及工程上什么时候该自己写、什么时候该用库都尽量说清楚。适合正在学数值方法、准备算法竞赛或者需要在嵌入式、无依赖环境里解方程的同学参考。1. 先搞懂算法本身它和高斯消元的真正区别1.1 从增广矩阵说起线性方程组可以写成 Ax b 的形式。A 是系数矩阵x 是未知数向量b 是常数向量。高斯-约当消元法的思路是把 A 和 b 拼成一个增广矩阵也就是在 A 右边直接多加一列常数项。比如x 2y z 82x y z 7x y 3z 12增广矩阵就是[ 1 2 1 | 8 ][ 2 1 1 | 7 ][ 1 1 3 | 12 ]中间的竖线只是给人看的程序里不需要单独存这一列矩阵宽度直接加 1 就行。解方程组的过程本质上是对这个矩阵做三种行变换交换两行、某行整体乘以非零常数、某行加上另一行的倍数。这三种操作都是“行等价变换”不会改变方程组的解集所以随便怎么折腾最后得到的矩阵和原方程是同一个解集的不同写法。我们的目标是把系数部分化成行最简形Reduced Row Echelon Form简称 RREF每个主元所在的列除了主元本身是 1 之外其他位置全是 0。一旦化成这种形态方程组的解相当于直接印在最后一列上。1.2 为什么我常用高斯-约当而不是经典高斯消元经典高斯消元是两步走先通过消元把系数矩阵化成上三角然后再从最后一个方程开始往上回代。高斯-约当则是“一步到位”在消元过程中不仅消掉主元下方的元素还把主元上方的元素也全消掉最后系数矩阵直接变成单位矩阵。为了方便对比我列个表对比项经典高斯消元高斯-约当消元最终形态上三角矩阵行最简形RREF求解方式回代直接读解乘加次数约 n³/3约 n³/2代码复杂度和边界情况回代部分容易写错统一循环逻辑简单求逆矩阵扩展需额外处理右侧直接在右侧拼单位阵计算量上高斯-约当确实比高斯消元多一点大概多 50% 的乘加操作。但很多场景里我们不只求一个方程组的解而是要求多个右端项甚至直接求矩阵的逆。这个时候高斯-约当有个天然优势把 A 和单位矩阵 I 拼成增广矩阵跑完一遍之后右侧就是 A⁻¹。一步到位没有额外代码。而且从工程角度看代码少一个回代环节就少一类维护边界条件的麻烦。对于 n 在几千以内的中小规模问题这点性能差距几乎可以忽略。可维护性和逻辑统一性反而更重要。1.3 无解、唯一解、无穷解怎么一眼看出行最简形化完以后解的情况可以通过三条规则直接判断无解某一行出现了“左侧全零、右侧非零”相当于数学上的 0 非零数矛盾。无穷解主元个数小于未知数个数说明存在自由变量。唯一解主元个数等于未知数个数每个变量都被唯一确定解就是最后一列的数字。例如某个方程组行最简形是[ 1 0 | 2 ][ 0 1 | 3 ]主元数是 2未知数也是 2唯一解 x2, y3。如果出现[ 1 0 0 | 1 ][ 0 0 0 | 2 ]第二行左侧全零、右侧是 2这就是无解。把这三条写进代码并不复杂真正复杂的是怎么在浮点数的世界里判断“全零”。这就要聊到 EPS 和主元策略了。2. C实现前必须想清楚的三个设计问题2.1 矩阵用什么存二维数组还是嵌套 vector我建议优先用vectorvectordouble不要自己管理二维裸数组。原因很实际行交换方便。C 的swap(mat[row1], mat[row2])可以直接交换两个 vector底层只是交换指针开销极小。矩阵行数、列数可以动态确定不用在编译期写死函数适配性更好。打印和调试直观循环起来就是一个标准的二维结构。有同学会担心性能。说实话几千阶以内的矩阵vectorvectordouble完全够用。如果你确实要处理超大矩阵又不想引入第三方库可以用一维 vector 自行模拟二维索引mat[i * cols j]。这样内存连续缓存友好度更高。不过付出的代价是代码可读性下降行交换需要手动实现内存拷贝。另外提醒一个小坑代码里涉及到矩阵行列数时最好强制转成int不要用size_t直接参与比较。否则for (int i 0; i mat.size() - 1; i)这种写法在mat为空时会爆炸因为size_t无符号减出一个巨大正整数。2.2 主元选择为什么不能偷懒消元的第一步是确定当前列用哪一行作为基准行这一行的主元列元素叫“主元”。最朴素的做法是从当前行往下找第一个非零元素把它当主元。这个思路在理论上没问题但在浮点运算里可能翻车。假设主元是 1e-300归一化时要让这行除以 1e-300结果就是把这个极小的数放大成 1同时它原本携带的相对误差也被放大到天文数字级别。后续消元操作会把误差扩散到其他行最终解出来的结果可能完全不能用。标准解法是“部分主元法”在当前列从当前行往下扫描找绝对值最大的那个元素所在行把它交换到当前行。这个策略的收益只是 O(n) 次比较但数值稳定性提升非常明显。你用 double 写代码如果不做部分主元遇到稍微病态一点的矩阵就等着看诡异结果吧。完整的主元法还要交换列实现复杂度高实际应用中除非做符号计算否则很少用到不推荐自己造轮子。2.3 浮点误差EPS 到底怎么取判断浮点数是不是 0绝不能直接if (x 0.0)。因为经过一系列加减乘除后原本该是 0 的位置很可能是 1e-15 这种微小残差。我的习惯是定义const double EPS 1e-9;所有需要判断为 0 的地方都用fabs(x) EPS。这个值不是拍脑袋定的double 大约有 15 位有效十进制数字1e-9 作为判断阈值比较平衡既能容忍合理误差又不容易把真正该保留的小数误判成 0。如果你的问题数值本身很大比如元素动辄 1e10 级别那 EPS 可以适当放大到 1e-7。反过来如果是严谨的数值实验更推荐long double或者干脆用分数类做精确有理数运算后面第五章会提到。3. 完整代码实现把 Gauss-Jordan 写成通用函数3.1 函数接口怎么设计我习惯用一个整数返回值表达最终状态返回 0唯一解。返回 1无穷多解。返回 -1无解。同时传入的矩阵会被原地修改成行最简形解通过引用参数solution输出。这样设计的好处是调用方可以根据返回值走不同分支同时还能拿到化简后的矩阵方便调试和后续处理。如果只需要判断一次“这方程组有没有解”你可以把返回值当状态码用如果还需要具体解的值就从solution里取。3.2 核心循环逐段拆解整个代码的骨架是一个两层循环外层遍历列内层做“找主元、交换、归一化、消去”四件事。找主元的逻辑从第row行开始往下扫描col列找到绝对值最大的元素记录所在行pivotRow。如果最大绝对值小于EPS说明这列在当前剩余行里全是 0这个变量是自由变量跳过这一列。把pivotRow和row交换让主元跑到当前行。归一化那一步是把第row行从第col列开始的所有元素都除以pivot。为什么从col开始而不是从 0 开始因为col左边的列都已经是格式化完毕的主元列再动它们也不会影响正确性还省了计算量。最后是消去操作。和高斯消元只消下边行不同这里要遍历所有行把主元列的元素全部消成 0。先判断factor的绝对值是否小于EPS如果本身已经接近 0 就跳过避免多余的浮点乘加把误差越滚越大。3.3 完整可运行代码下面是完整代码依赖只有标准库任何支持 C11 的编译器都能直接编。#include iostream #include vector #include cmath #include iomanip using namespace std; const double EPS 1e-9; using Matrix vectorvectordouble; void printMatrix(const Matrix mat, int precision 6) { for (const auto row : mat) { for (double val : row) { cout setw(12) fixed setprecision(precision) val ; } cout endl; } } // 返回值 // 0 - 唯一解 // 1 - 无穷多解 // -1 - 无解 int gaussJordan(Matrix mat, vectordouble solution) { int m (int)mat.size(); if (m 0) return -1; int n (int)mat[0].size() - 1; // 未知数个数 列数 - 1 int row 0; for (int col 0; col n row m; col) { // 1. 部分主元在当前列中找绝对值最大的元素所在行 int pivotRow row; double maxAbs fabs(mat[row][col]); for (int i row 1; i m; i) { if (fabs(mat[i][col]) maxAbs) { maxAbs fabs(mat[i][col]); pivotRow i; } } // 如果这一列从 row 行起全是 0跳过这一列 if (maxAbs EPS) { continue; } // 2. 把主元行交换到当前 row 行 if (pivotRow ! row) { swap(mat[pivotRow], mat[row]); } // 3. 归一化让 mat[row][col] 变成 1 double pivot mat[row][col]; for (int j col; j n; j) { mat[row][j] / pivot; } // 4. 消去其余所有行的 col 列元素 for (int i 0; i m; i) { if (i row) continue; double factor mat[i][col]; if (fabs(factor) EPS) continue; for (int j col; j n; j) { mat[i][j] - factor * mat[row][j]; } } row; } // 检查无解某一行系数全为 0但常数项非 0 for (int i row; i m; i) { bool allZero true; for (int j 0; j n; j) { if (fabs(mat[i][j]) EPS) { allZero false; break; } } if (allZero fabs(mat[i][n]) EPS) { return -1; } } // 主元数量 row 小于未知数数量 n意味着有自由变量 if (row n) { return 1; } // 唯一解解就在最后一列 solution.assign(n, 0.0); for (int i 0; i n; i) { solution[i] mat[i][n]; } return 0; } int main() { cout 示例1唯一解 endl; Matrix a1 { {1, 2, 1, 8}, {2, 1, 1, 7}, {1, 1, 3, 12} }; vectordouble sol; int ret gaussJordan(a1, sol); printMatrix(a1); if (ret 0) { cout 唯一解; for (int i 0; i (int)sol.size(); i) { cout x i 1 sol[i] ; } cout endl; } else if (ret 1) { cout 无穷多解 endl; } else { cout 无解 endl; } cout \n 示例2无解 endl; Matrix a2 { {1, 1, 1}, {1, 1, 3} }; ret gaussJordan(a2, sol); printMatrix(a2); cout 返回码 ret - (ret -1 ? 无解 : 有解) endl; cout \n 示例3无穷解 endl; Matrix a3 { {1, 1, 1, 1}, {2, 2, 2, 2} }; ret gaussJordan(a3, sol); printMatrix(a3); cout 返回码 ret - (ret 1 ? 无穷多解 : 其他) endl; return 0; }这段代码我实测过可以直接跑。几个细节说一下swap(mat[pivotRow], mat[row])交换的是整个 vector不是逐元素交换效率高。归一化和消去都从col开始因为更左边的列已经处理完毕动了也不影响但没必要浪费时间。消去循环里先判断factor绝对值避免把小误差当作有效数去乘加能减少不少无谓运算。4. 测试运行与结果分析验证代码到底靠不靠谱4.1 唯一解案例三个方程三个未知数x 2y z 82x y z 7x y 3z 12理论上解是 x1, y2, z3。代码跑完增广矩阵会变成[ 1 0 0 | 1 ][ 0 1 0 | 2 ][ 0 0 1 | 3 ]打印出来如果看到对角线附近是 1、其他位置是接近 0 的小数比如 1e-15 级别就说明归一化和消去过程正确。浮点运算不可能精确得到 0允许出现微小残差。4.2 无解案例两个方程x y 1x y 3增广矩阵是[ 1 1 | 1 ][ 1 1 | 3 ]代码会先把第一行归一化成 [1, 1, 1]然后用它消第二行得到[ 1 1 | 1 ][ 0 0 | 2 ]第二行对应的方程是 0 2无解。程序走到第 4.1 步检查时发现左侧全零、右侧非零直接返回 -1。这种情况在现实里常见于数据采集有矛盾的过定方程组解出来没有数学意义硬算只会得到一个“最小二乘意义”的近似解。这也是为什么状态码设计要单独区分无解而不是简单给个空 vector。4.3 无穷解案例方程组x y z 12x 2y 2z 2第二个方程其实是第一个的两倍信息完全冗余。代码跑完第二行会被消成全零。主元数量是 1未知数是 3返回无穷多解。如果想要进一步描述这个解集可以从行最简形中读出一个可行解x 1 - y - z其中 y 和 z 是自由变量。代码层面如果确实需要把通解形式也输出那就要额外记录自由变量索引和主元变量索引这部分内容我觉得大多数场景用不上就没有写进函数里。你需要的时候可以在返回1的分支里继续扩展。4.4 多元组测试时的小建议我测试的时候习惯把三组样例写进一个 main用printMatrix看每一步结果。你复制代码运行时如果看到某些位置出现 -0.000000这不是错误是浮点运算里的负零视觉上吓人而已。介意的可以在打印函数里加一行判断如果fabs(val) EPS就输出 0。这块我后来踩过一次小坑明明逻辑没问题看到-0.000000总觉得是不是代码哪里写错浪费了十分钟排查。所以打印函数里把近零数统一成 0是很有必要的调试友好性优化。5. 躲开这些坑精度、效率与工程替代方案5.1 浮点误差的典型症状与对策最常见的症状是明明手算整数解程序解出来却是 0.999999999 或者 1.000000001。这属于 normal 现象不是代码逻辑错误。如果你发现程序结果和手算结果差得离谱优先检查三点有没有做部分主元没做的话遇到主元接近 0 的矩阵很容易爆炸。EPS 设置是否合理太小会导致“接近 0 但被当成非 0”太大会把有效小数误杀。数据本身是否病态比如 Hilbert 矩阵 H[i][j] 1 / (i j 1)n 到 10 左右用 double 做高斯-约当得到的解可能惨不忍睹。这是数值线性代数里非常经典的问题单靠选主元只能缓解不能根治。遇到病态矩阵我的建议是不要死磕自写实现直接上库。5.2 数据规模与性能优化建议这个自写实现的时间复杂度是 O(m × n²)m 是方程数n 是未知数个数。如果 m 和 n 都在几千以内运行时间一般可以接受再大就要认真考虑性能问题了。几个优化方向把vectorvectordouble改成一维数组保证内存连续利用 CPU 缓存。如果只求一组解可以用经典高斯消元加回代比高斯-约当省一点乘加操作。如果同一个系数矩阵要配十几个不同的右端项 b那么 Gauss-Jordan 反而有优势一次消元做完行最简形所有右侧向量同时搞定。开启编译优化-O2对循环优化效果非常明显。5.3 什么时候该换库嵌入式、算法竞赛笔试、教学演示、依赖受限的场景手写实现完全够用。但如果你在公司项目里处理真实数值计算我更推荐直接用现成库EigenC 头文件库PartialPivLU、FullPivLU接口清爽数值稳定性有保障。Armadillo语法接近 MATLAB适合快速原型。LAPACK老牌 Fortran 库dgesv性能极其稳定适合大规模矩阵。这些库内部不只是简单的高斯消元还包含矩阵条件数估计、迭代改进、分块算法等优化。遇到明显病态的问题自己写代码和大厂级的数值库相比稳定性差距不是同一个量级。需要注意的是我说这些不是为了贬低手写实现恰恰相反面试、讲课、理解原理时手写一遍高斯-约当才是最能加深理解的方式。哪怕你以后长期用 Eigen你也得知道FullPivLU到底在干什么才能判断什么时候该用部分主元什么时候该用完整主元。根据我个人写代码和调 bug 的习惯最后再分享一个小改造思路把代码里的double换成自定义分数类实现加减乘除、取绝对值、比较大小这些运算符就可以做精确有理数消元。我拿这套改进版去验证普通 double 版本算出来的结果很多精度问题立刻现出原形。高斯-约当本身不复杂真正体现功力的地方在于你对边界情况的敬畏和对数值误差的理解。