C++实现高斯消元法：从数学原理到工业级求解器
1. 项目概述：从线性方程组到代码实现
高斯消元法，这个名字对于任何学过线性代数或者接触过算法竞赛的人来说，都再熟悉不过了。它不仅是求解线性方程组最经典、最核心的算法，更是理解矩阵运算、线性空间乃至整个计算数学的一块基石。我第一次在工程中真正用上它，是在做一个简单的物理引擎碰撞检测模块时，需要快速求解一个3x3的线性系统来计算碰撞后的速度变化。当时我尝试了各种库，但总觉得不够“透明”，于是决定自己动手实现一遍。这一实现，才发现这个看似教科书式的算法，在实际编码中藏着不少门道——从浮点数精度这个“老冤家”，到各种特殊矩阵的处理策略，每一步都考验着程序员对数学和计算机科学的双重理解。
今天，我们就来彻底拆解高斯消元法，并用C++从头实现一个工业级可用的版本。我们不止步于“能跑”，更要追求“稳健”和“高效”。无论你是正在学习线性代数的学生，还是需要在图形学、机器学习或科学计算中处理线性系统的开发者，这篇文章都将带你从原理到实践，从基础实现到高级优化，完整地走一遍。你会发现，实现一个正确的高斯消元程序，远比调用numpy.linalg.solve()要有趣和深刻得多。
2. 核心原理：高斯消元法的数学骨架
在动手写代码之前，我们必须把算法背后的数学逻辑吃透。高斯消元法的目标非常明确：将一个普通的线性方程组，通过一系列行变换，转化成一个上三角矩阵（行阶梯形），然后通过回代求解。这个过程本质上是在做矩阵的初等行变换。
2.1 算法步骤分解
高斯消元法通常分为两个阶段：前向消元和回代。
前向消元的目标，是将系数矩阵A化为上三角矩阵。对于一个有n个未知数的方程组，这个过程需要进行n-1步。在第k步（k从0开始），我们的任务是消除第k列中，第k行以下所有元素。具体操作是：对于i从k+1到n-1，我们用第i行减去第k行的某个倍数，使得A[i][k]变为0。这个倍数就是factor = A[i][k] / A[k][k]。这里，A[k][k]被称为第k步的主元。
回代过程则从最后一行开始，自底向上求解。对于上三角矩阵U，最后一行方程是U[n-1][n-1] * x[n-1] = b[n-1]，可以直接解出x[n-1] = b[n-1] / U[n-1][n-1]。然后，将已知的x[n-1]代入倒数第二行方程，解出x[n-2]，依此类推，直到求出所有未知数。
2.2 选主元：算法稳健性的关键
细心的你可能已经发现了一个潜在的问题：如果主元A[k][k]恰好为0或者非常接近0，那么计算factor时就会导致除零错误或者引入巨大的舍入误差，使得结果完全失真。这就是朴素高斯消元法的最大缺陷。
为了解决这个问题，我们必须引入选主元技术。最常用的有两种：
- 部分选主元：在第
k步，我们并不直接使用A[k][k]作为主元，而是在第k列中，从第k行到第n-1行，寻找绝对值最大的元素所在的行pivot_row。然后，交换第k行和第pivot_row行。这样能确保主元的绝对值尽可能大，显著提高数值稳定性。 - 完全选主元：不仅在第
k列中找，而是在右下角的子矩阵（从第k行第k列开始）中寻找绝对值最大的元素作为主元。找到后，需要同时交换行和列。完全选主元稳定性最高，但因为它改变了未知数的顺序（列交换对应未知数重排），需要在回代后对解向量进行重排，实现更复杂。对于绝大多数工程问题，部分选主元已经足够稳健。
核心经验：永远不要使用无选主元的高斯消元法。这是我在早期项目中踩过的一个大坑。一个看起来良性的矩阵，可能因为一个极小的主元而导致整个解向量溢出或产生荒谬的结果。部分选主元是性价比最高的选择。
2.3 算法复杂度分析
高斯消元法的时间复杂度是 O(n³)，其中 n 是方程个数（也是未知数个数）。具体来说，前向消元过程是一个三重循环，主导了计算成本。空间复杂度，如果我们原地修改矩阵A和向量b，则是 O(n²)。这个复杂度决定了它适用于中小规模（n 在几千以内）的稠密线性系统。对于更大规模或稀疏的系统，需要迭代法（如共轭梯度法）或专门的稀疏矩阵求解器。
3. C++实现：从零搭建稳健求解器
理解了原理，我们开始用C++实现。我们的目标是编写一个函数vector<double> gaussianElimination(vector<vector<double>>& A, vector<double>& b)，它接收一个 n x n 的系数矩阵A和一个长度为 n 的常数向量b，返回解向量x。我们将实现部分选主元，并仔细处理数值精度问题。
3.1 基础数据结构与函数签名
首先，我们确定使用std::vector来表示矩阵和向量，这是C++中最灵活和通用的选择。为了避免拷贝开销，我们使用引用传递。
#include <vector> #include <cmath> #include <algorithm> #include <stdexcept> #include <iostream> // 定义一个全局的精度容忍度，用于判断浮点数是否为0 const double EPS = 1e-10; /** * @brief 使用部分选主元的高斯消元法求解线性方程组 Ax = b * @param A n x n 的系数矩阵 (将被修改) * @param b 右侧常数向量，长度为 n (将被修改) * @return 解向量 x * @throws std::runtime_error 如果矩阵是奇异的（无解或无穷多解） */ std::vector<double> gaussianElimination(std::vector<std::vector<double>>& A, std::vector<double>& b) { int n = A.size(); // 基础校验 if (n == 0) return {}; if (A[0].size() != n) { throw std::invalid_argument("Matrix A must be square."); } if (b.size() != n) { throw std::invalid_argument("Vector b size must match matrix dimension."); } std::vector<double> x(n, 0.0); // 主消元过程将在后续填充 return x; }3.2 前向消元与部分选主元实现
这是算法的核心部分。我们需要一个外层循环k遍历每一列（作为主元列），在每一列内部，先进行选主元，然后进行消元。
std::vector<double> gaussianElimination(std::vector<std::vector<double>>& A, std::vector<double>& b) { int n = A.size(); // ... 参数校验 ... // 前向消元过程 for (int k = 0; k < n; ++k) { // --- 部分选主元 --- int pivot_row = k; double max_val = std::fabs(A[k][k]); for (int i = k + 1; i < n; ++i) { if (std::fabs(A[i][k]) > max_val) { max_val = std::fabs(A[i][k]); pivot_row = i; } } // 如果找到的主元绝对值太小，视为奇异矩阵 if (max_val < EPS) { throw std::runtime_error("Matrix is singular or nearly singular. No unique solution exists."); } // 交换当前行 k 和主元行 pivot_row if (pivot_row != k) { std::swap(A[k], A[pivot_row]); std::swap(b[k], b[pivot_row]); // 注意：在实际的数值线性代数库中，我们会记录行交换的次数， // 以确定行列式的符号。这里我们只关心解。 } // --- 消去第 k 列下方所有元素 --- // 将主元行归一化（可选，但能提高后续计算的数值稳定性） // double pivot = A[k][k]; // 此时主元已是非零最大值 // for (int j = k; j < n; ++j) { A[k][j] /= pivot; } b[k] /= pivot; // 注意：我们更常见的做法是不显式归一化，而是在计算乘数时做除法。 for (int i = k + 1; i < n; ++i) { // 计算乘数：第 i 行要减去主元行的多少倍，才能使 A[i][k] 为 0 double factor = A[i][k] / A[k][k]; // 如果乘数已经为0，可以跳过该行以节省计算（微优化） if (std::fabs(factor) < EPS) continue; // 消去第 i 行第 k 列的元素，并更新该行后面的列和对应的 b // 从第 k 列开始消去，因为前面的列已经是0了 for (int j = k; j < n; ++j) { A[i][j] -= factor * A[k][j]; } b[i] -= factor * b[k]; // 显式地将 A[i][k] 置为 0，避免浮点误差积累导致非零（可选，但清晰） // A[i][k] = 0.0; } } // 前向消元结束，此时 A 已成为上三角矩阵关键细节解析：
- 选主元的时机：必须在计算乘数
factor之前完成。我们寻找的是第k列中从第k行到末尾的绝对值最大者。 - 奇异矩阵判断：当所有候选主元的绝对值都小于
EPS时，意味着矩阵的秩不足，方程组要么无解，要么有无穷多解。我们的简单实现选择抛出异常。更完善的实现可以返回一个状态码或尝试计算最小二乘解。 - 交换操作：
std::swap(A[k], A[pivot_row])交换的是两个vector<double>，即两行数据，效率很高。同时别忘了交换常数向量b中对应的元素。 - 消元循环的起始列：内层循环
for (int j = k; ...)从j = k开始。因为经过前k-1步消元后，第i行第0到k-1列的元素理论上已经是0了（尽管可能有微小误差），从k开始计算可以节省一半的操作。这是算法从 O(n³) 优化到 (2/3)n³ 的关键。
3.3 回代过程实现
前向消元结束后，我们得到了一个上三角矩阵。现在从最后一行开始，反向求解。
// ... 前向消元代码之后 ... // 回代过程 for (int i = n - 1; i >= 0; --i) { // 先计算等式右边：b[i] 减去已知解与对应系数的乘积之和 double sum = 0.0; for (int j = i + 1; j < n; ++j) { sum += A[i][j] * x[j]; } // 求解 x[i] x[i] = (b[i] - sum) / A[i][i]; // 一个可选的稳健性检查：解出后再次检查分母 if (std::fabs(A[i][i]) < EPS) { // 虽然前向消元已做检查，但浮点运算可能使本应非零的主元变得极小 throw std::runtime_error("Zero pivot encountered during back substitution."); } } return x; } // 函数结束回代要点：
- 循环顺序：必须从下往上 (
i = n-1to0)。 - 求和变量：对于每个
x[i]，需要累加A[i][j] * x[j]，其中j从i+1到n-1，这些都是已经求出的解。 - 除法检查：尽管前向消元保证了主元非零，但浮点舍入可能使其变得极小。在回代时再次检查
A[i][i]是一个良好的防御性编程习惯。
4. 处理边界情况与精度问题
一个健壮的实现必须考虑各种边界情况和数值陷阱。
4.1 奇异与病态矩阵
我们的代码通过检查主元是否接近零来判断奇异矩阵。但还有一种更隐蔽的情况：病态矩阵。即使矩阵非奇异，如果条件数（最大奇异值与最小奇异值之比）非常大，微小的输入误差或舍入误差会被极度放大，导致解不可信。
// 一个简单的病态矩阵例子：希尔伯特矩阵 // H[i][j] = 1.0 / (i + j + 1)，当阶数增大时，条件数急剧增加。 std::vector<std::vector<double>> createHilbertMatrix(int n) { std::vector<std::vector<double>> H(n, std::vector<double>(n)); for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { H[i][j] = 1.0 / (i + j + 1.0); // 注意使用 1.0 避免整数除法 } } return H; }对于病态问题，高斯消元法（即使有选主元）可能给出错误结果。更专业的数值计算库会提供条件数估计，或者推荐使用迭代 refinement 等技术来改善解的质量。
4.2 浮点数比较与EPSILON的选择
这是数值计算中最经典的坑。永远不要用==直接比较浮点数。我们使用一个全局的EPS（epsilon）作为容忍度。EPS的选择至关重要：
- 太小（如
1e-15）：可能因为舍入误差将本应视为零的值误判为非零，导致算法失败。 - 太大（如
1e-6）：可能将实际很小的主元误判为零，错误地抛出奇异异常。
通常，对于双精度double，EPS取1e-10到1e-12是一个合理的范围。一个更稳健的策略是使用相对误差：fabs(a) < EPS * max_pivot，其中max_pivot是当前列中找到的最大主元值。这能自适应矩阵元素的尺度。
// 改进的选主元与奇异性判断（相对误差版） int pivot_row = k; double max_val = std::fabs(A[k][k]); for (int i = k + 1; i < n; ++i) { double val = std::fabs(A[i][k]); if (val > max_val) { max_val = val; pivot_row = i; } } // 使用相对容忍度判断 if (max_val < EPS * norm_of_A) { // norm_of_A 可以是矩阵 A 的某一范数（如无穷范数） throw std::runtime_error("Matrix is singular or nearly singular."); }4.3 原地操作与内存布局
我们的实现是原地修改矩阵A和向量b的。这节省了内存，但破坏了原始数据。如果调用者需要保留原始矩阵，需要在函数内部进行拷贝。
std::vector<double> gaussianElimination(const std::vector<std::vector<double>>& A_orig, const std::vector<double>& b_orig) { // 创建副本 auto A = A_orig; auto b = b_orig; // ... 后续消元求解过程 ... }此外，std::vector<std::vector<double>>这种“向量的向量”表示法，可能导致内存不连续，影响缓存效率。对于高性能计算，使用一维数组并按行优先或列优先顺序存储是更好的选择。
5. 高级话题与扩展应用
基础的高斯消元法远不止于求解Ax = b。
5.1 计算矩阵的行列式
高斯消元过程中，矩阵的行列式很容易得到。对于一个上三角矩阵U，其行列式等于主对角线上所有元素的乘积。在进行行交换时，行列式的符号会改变。因此，在消元过程中，我们只需要记录行交换的次数，并在最后乘上(-1)^(交换次数)。
double determinant(std::vector<std::vector<double>> A) { // 传值，不修改原矩阵 int n = A.size(); double det = 1.0; int swap_count = 0; for (int k = 0; k < n; ++k) { // 选主元... if (pivot_row != k) { std::swap(A[k], A[pivot_row]); swap_count++; } det *= A[k][k]; // 累乘主元 // 如果主元为0，行列式就是0 if (fabs(A[k][k]) < EPS) return 0.0; // 消元过程... for (int i = k + 1; i < n; ++i) { double factor = A[i][k] / A[k][k]; for (int j = k + 1; j < n; ++j) { // 注意：计算行列式时，可以只更新下三角，但这里为清晰展示 A[i][j] -= factor * A[k][j]; } // A[i][k] 被消为0，不影响行列式 } } // 根据行交换次数调整符号 if (swap_count % 2 == 1) { det = -det; } return det; }5.2 求解矩阵的逆
求一个方阵A的逆矩阵A⁻¹，等价于求解n个线性方程组：A * X = I，其中I是单位矩阵，X的每一列就是A⁻¹的对应列。我们可以通过高斯-若尔当消元法一次性完成。其思路是将增广矩阵[A | I]通过行变换化为[I | A⁻¹]。
std::vector<std::vector<double>> matrixInverse(const std::vector<std::vector<double>>& A) { int n = A.size(); // 创建增广矩阵 Aug = [A | I] std::vector<std::vector<double>> Aug(n, std::vector<double>(2 * n, 0.0)); for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { Aug[i][j] = A[i][j]; } Aug[i][n + i] = 1.0; // 单位矩阵部分 } // 高斯-若尔当消元 for (int k = 0; k < n; ++k) { // 选主元 (在左侧 n 列中选) int pivot_row = k; // ... 选主元逻辑 ... // 交换行 // 归一化主元行：让主元变为1 double pivot = Aug[k][k]; for (int j = k; j < 2 * n; ++j) { Aug[k][j] /= pivot; } // 消去其他行的第 k 列元素（不仅是下方，是全部其他行） for (int i = 0; i < n; ++i) { if (i != k) { double factor = Aug[i][k]; for (int j = k; j < 2 * n; ++j) { Aug[i][j] -= factor * Aug[k][j]; } } } } // 提取逆矩阵（增广矩阵的右半部分） std::vector<std::vector<double>> invA(n, std::vector<double>(n)); for (int i = 0; i < n; ++i) { for (int j = 0; j < n; ++j) { invA[i][j] = Aug[i][n + j]; } } return invA; }高斯-若尔当消元法将消元和回代合并为一步，直接得到简化行阶梯形（单位矩阵）。
5.3 模意义下的高斯消元（异或方程组）
在算法竞赛中，经常遇到系数和未知数取值仅为0或1的线性方程组，运算在模2（即异或）下进行。例如，解决一些开关灯、博弈或图论问题。这时，算法可以极大简化，因为只有加（异或）和乘（与）操作，且不需要处理浮点数精度。
// 使用 std::bitset 优化模2高斯消元 const int MAXN = 1000; // 根据问题规模调整 std::bitset<MAXN> matrix[MAXN]; // 每一行是一个bitset std::vector<int> xorGaussianElimination(int n, int m) { // n个未知数，m个方程 std::vector<int> solution(n, 0); std::vector<int> where(n, -1); // 记录每个未知数被哪个方程的主元锁定 for (int col = 0, row = 0; col < n && row < m; ++col) { // 寻找当前列中第row行及以下为1的行 int sel = row; for (; sel < m; ++sel) { if (matrix[sel][col]) break; } if (sel == m) continue; // 该列全为0，是自由变元 // 交换行 std::swap(matrix[sel], matrix[row]); // 记录主元位置 where[col] = row; // 用第row行消去其他所有行第col列的1 for (int i = 0; i < m; ++i) { if (i != row && matrix[i][col]) { matrix[i] ^= matrix[row]; } } ++row; } // 构造解 for (int i = 0; i < n; ++i) { if (where[i] != -1) { solution[i] = matrix[where[i]][n]; // 假设增广矩阵最后一列是常数项 } // else: solution[i] 是自由变元，可任意赋值（通常赋0） } // 检查无解情况：如果存在一行系数全0但常数项为1 for (int i = 0; i < m; ++i) { bool all_zero = true; for (int j = 0; j < n; ++j) { if (matrix[i][j]) { all_zero = false; break; } } if (all_zero && matrix[i][n]) { return {}; // 返回空向量表示无解 } } return solution; }std::bitset的位运算效率极高，可以将算法复杂度从 O(n³) 优化到 O(n³ / ω)，其中 ω 是机器字长（通常为32或64）。
6. 性能优化与工程实践
在真实项目中，我们很少直接使用这种朴素的实现。以下是一些优化思路：
6.1 内存访问优化
二维vector的逐列访问缓存不友好。我们可以将矩阵按行优先存储在一维数组中，并通过A[i * n + j]访问元素。在消元的内层循环中，对连续内存的访问能极大提升缓存命中率。
// 一维数组表示矩阵 std::vector<double> A_flat(n * n); // 访问 A[i][j] 变为 A_flat[i * n + j] // 在消元循环中，对主元行 `k` 的访问是连续的，对目标行 `i` 的访问也是连续的。 for (int j = k; j < n; ++j) { A_flat[i * n + j] -= factor * A_flat[k * n + j]; }6.2 循环展开与指令级并行
现代CPU有流水线和SIMD指令（如SSE, AVX）。编译器在优化级别高时（如-O3）会自动进行一些循环展开和向量化。但我们也可以手动提示，例如将最内层的循环按固定步长展开。
// 手动循环展开示例（假设n是4的倍数） int j = k; for (; j + 3 < n; j += 4) { A[i][j] -= factor * A[k][j]; A[i][j+1] -= factor * A[k][j+1]; A[i][j+2] -= factor * A[k][j+2]; A[i][j+3] -= factor * A[k][j+3]; } for (; j < n; ++j) { A[i][j] -= factor * A[k][j]; }6.3 分块算法
对于非常大的矩阵，可以利用分块高斯消元法将矩阵划分为小块，使得计算更多地发生在高速缓存中，减少内存带宽的压力。这是专业数值库（如LAPACK）采用的核心技术之一。其基本思想是将矩阵分块，例如：
- 将矩阵划分为
A11, A12; A21, A22四个块。 - 对
A11进行LU分解（或高斯消元）。 - 更新
A12和A21。 - 更新
A22（称为 Schur 补）。 - 递归地对更新后的
A22进行分解。
6.4 使用现有的数值线性代数库
在绝大多数生产环境中，重新造轮子是不明智的。成熟的库经过了几十年的优化和测试。
- Eigen: C++模板库，接口优雅，性能卓越，易于集成。
- Armadillo: 语法类似MATLAB，对科学计算友好。
- LAPACK(通过MKL, OpenBLAS等调用): 工业标准，性能最强，但C接口稍显复杂。
使用Eigen求解Ax=b只需几行代码：
#include <Eigen/Dense> Eigen::MatrixXd A(n, n); Eigen::VectorXd b(n); // ... 填充 A 和 b ... Eigen::VectorXd x = A.partialPivLu().solve(b); // 使用部分主元LU分解 // 或者 A.colPivHouseholderQr().solve(b) 使用更稳定的分解7. 测试与验证
编写完算法，必须进行全面的测试。
7.1 构造测试用例
- 随机可逆矩阵：生成随机矩阵，并确保其行列式远离零。
- 已知解的问题：构造矩阵
A和解向量x_true，计算b = A * x_true，然后用我们的算法解Ax=b，比较x和x_true。 - 特殊矩阵：
- 对角占优矩阵：确保稳定性。
- 希尔伯特矩阵：测试病态情况。
- 对称正定矩阵：可以用Cholesky分解更高效，但高斯消元也应能解。
- 奇异矩阵：测试算法是否能正确检测并报告错误。
7.2 验证结果
由于浮点误差，我们不能直接比较x == x_true。应该计算残差r = b - A * x的范数（如L2范数），以及误差e = x - x_true的范数。
double residualNorm(const std::vector<std::vector<double>>& A, const std::vector<double>& b, const std::vector<double>& x) { int n = A.size(); double norm = 0.0; for (int i = 0; i < n; ++i) { double sum = 0.0; for (int j = 0; j < n; ++j) { sum += A[i][j] * x[j]; } double r = b[i] - sum; norm += r * r; } return std::sqrt(norm); } // 在测试中 double residual = residualNorm(A_original, b_original, x_solution); if (residual > 1e-8) { // 设置一个合理的容忍度 std::cerr << "Warning: Large residual: " << residual << std::endl; }7.3 性能基准测试
对于不同规模的矩阵（如 n=100, 500, 1000），测试算法的运行时间，验证其 O(n³) 的增长趋势。可以使用<chrono>库进行计时。
#include <chrono> auto start = std::chrono::high_resolution_clock::now(); auto x = gaussianElimination(A, b); auto end = std::chrono::high_resolution_clock::now(); auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start); std::cout << "Time elapsed: " << duration.count() << " ms" << std::endl;8. 完整代码示例与总结
将以上所有部分整合，我们得到一个相对完整、带有部分选主元和基本错误处理的高斯消元法C++实现。这个实现侧重于清晰和教育意义，在追求极致性能的场合还需要进一步优化。
#include <iostream> #include <vector> #include <cmath> #include <algorithm> #include <stdexcept> #include <random> class GaussianEliminationSolver { public: static constexpr double EPS = 1e-12; static std::vector<double> solve(std::vector<std::vector<double>> A, std::vector<double> b) { int n = A.size(); if (n == 0) return {}; if (A[0].size() != n || b.size() != n) { throw std::invalid_argument("Invalid matrix or vector dimensions."); } std::vector<double> x(n, 0.0); std::vector<int> row_permutation(n); // 可选：记录行交换，用于更高级的用途 std::iota(row_permutation.begin(), row_permutation.end(), 0); // 前向消元 for (int k = 0; k < n; ++k) { // 部分选主元 int pivot_row = k; double max_val = std::fabs(A[k][k]); for (int i = k + 1; i < n; ++i) { double val = std::fabs(A[i][k]); if (val > max_val) { max_val = val; pivot_row = i; } } if (max_val < EPS) { throw std::runtime_error("Matrix is singular or ill-conditioned."); } if (pivot_row != k) { std::swap(A[k], A[pivot_row]); std::swap(b[k], b[pivot_row]); std::swap(row_permutation[k], row_permutation[pivot_row]); } // 消去下方行 for (int i = k + 1; i < n; ++i) { double factor = A[i][k] / A[k][k]; // 如果乘数很小，可以跳过以减少操作（但可能引入误差） if (std::fabs(factor) < EPS) continue; // 从第k列开始消去 for (int j = k; j < n; ++j) { A[i][j] -= factor * A[k][j]; } b[i] -= factor * b[k]; // 可选：显式置零，增加可读性，对结果无影响 // A[i][k] = 0.0; } } // 回代 for (int i = n - 1; i >= 0; --i) { double sum = 0.0; for (int j = i + 1; j < n; ++j) { sum += A[i][j] * x[j]; } x[i] = (b[i] - sum) / A[i][i]; } return x; } // 一个简单的测试函数 static void test() { std::mt19937 rng(42); std::uniform_real_distribution<double> dist(-10.0, 10.0); int n = 5; std::vector<std::vector<double>> A(n, std::vector<double>(n)); std::vector<double> x_true(n); std::vector<double> b(n); // 生成随机可逆矩阵和真解 for (int i = 0; i < n; ++i) { x_true[i] = dist(rng); for (int j = 0; j < n; ++j) { A[i][j] = dist(rng); } // 使矩阵对角占优，确保可逆 A[i][i] += n * 10.0; } // 计算 b = A * x_true for (int i = 0; i < n; ++i) { b[i] = 0.0; for (int j = 0; j < n; ++j) { b[i] += A[i][j] * x_true[j]; } } try { auto x_solved = solve(A, b); std::cout << "True solution x_true:\n"; for (double val : x_true) std::cout << val << " "; std::cout << "\nComputed solution x_solved:\n"; for (double val : x_solved) std::cout << val << " "; std::cout << std::endl; // 计算误差 double max_err = 0.0; for (int i = 0; i < n; ++i) { max_err = std::max(max_err, std::fabs(x_solved[i] - x_true[i])); } std::cout << "Maximum absolute error: " << max_err << std::endl; if (max_err < 1e-8) { std::cout << "Test PASSED!" << std::endl; } else { std::cout << "Test FAILED!" << std::endl; } } catch (const std::exception& e) { std::cerr << "Error: " << e.what() << std::endl; } } }; int main() { GaussianEliminationSolver::test(); return 0; }实现一个高斯消元法，就像搭建一座微型的数值计算桥梁。它连接了抽象的线性代数理论和具体的计算机运算。通过这个过程，我们不仅学会了如何解方程，更深入理解了浮点运算的陷阱、算法稳定性的重要性和代码优化的多种维度。下次当你面临一个需要求解线性系统的问题时，不妨先想想是否真的需要搬出庞大的线性代数库，也许这个几百行的手写实现，就是最简洁、最直接的解决方案。至少在理解原理和调试问题的时候，拥有一个自己完全掌控的“干净”实现，会让你心里踏实得多。