核心概述
针对同一矩阵 A 的重复求解场景，先进行 LU/LUP 分解，将后续每次求解 Ax=b 的代价降为 O(n^2)，并在实数域支持部分主元、在模域以乘逆元替代除法。
原文引述
Description: LU decomposition (with optional partial pivoting) for repeated solves Ax=b. After O(n^3) factorization, each solve costs O(n^2). Supports modular inverses for mod p. Time: Factor O(n^3); each solve O(n^2) Status: tested ——摘自本节点现有英文注释
-
函数/接口（据实现）
- struct LU { vector<vector
> L, U; vector P; }; - LU decompLU(const vector<vector
>& A) - 返回 L、U 与置换向量 P（若需要部分主元）。
- vector
solveLU(const LU& lu, const vector & b) - 使用预分解的 LU 快速求解 Ax=b。
- 可选：solveMultipleLU(lu, B) 同时解多个右端向量（矩阵 B）。
- struct LU { vector<vector
-
LU 分解（无主元或部分主元）
- 目标：A = P·L·U，其中 P 为置换矩阵，L 为单位下三角，U 为上三角。
- 算法（带部分主元）：
- 初始化 L=I，U=A，P 为单位置换。
- 对 k=0..n−1：
- 选取主元行 p（|U[p][k]| 最大）；若 p≠k，交换 U 的行 k 与 p，记录 P[k]=p。
- 对 i=k+1..n−1：
- L[i][k] = U[i][k] / U[k][k]（模域乘逆元）。
- 对 j=k..n−1：U[i][j] = U[i][j] − L[i][k]·U[k][j]。
- 若不需要主元（模域且已知可逆），可跳过主元选择。
-
求解步骤（前代 + 回代）
- 置换右端：b’ = P·b（按 P 向量重排）。
- 前代：解 L·y = b’：
- y[0] = b’[0]；
- y[i] = b’[i] − Σ_{j=0}^{i−1} L[i][j]·y[j]。
- 回代：解 U·x = y：
- x[n−1] = y[n−1] / U[n−1][n−1]；
- x[i] = (y[i] − Σ_{j=i+1}^{n-1} U[i][j]·x[j]) / U[i][i]。
- 模域下除法改为乘逆元。
-
数值与实现注意
- 实数域：部分主元提升数值稳定；若 U[k][k] 极小，可能矩阵奇异。
- 模域：若模为素数，所有非零元素可逆；非素数模需检查 gcd。
- 原地分解：可将 L 的严格下三角与 U 存放在同一矩阵中（U 的对角线不存，因 L 对角线为 1）。
- 多右端向量：可批量前代与回代，减少重复。
-
输出语义
- decompLU 返回 LU 结构；solveLU 返回解向量 x。
-
复杂度
- 分解 O(n^3)；每次求解 O(n^2)。
-
使用与注意
- 适合同一矩阵 A、多个 b 的场景（如线性规划、矩阵求逆）。
- 若仅需一次求解，直接高斯消元更简洁（见 114-数值算法-线性方程求解）。
- 模 2 系统见 116-数值算法-GF2线性求解（位运算加速）。
- 三对角系统见 117-数值算法-三对角求解（O(n)）。
- 求逆矩阵：LU 分解后解 n 个单位向量可得逆（见 107-数值算法-矩阵求逆）。