带主元置换的自定义LU分解大矩阵运行结果异常问题求助

安全阿岳

2026-9-30

代码错误排查与修复方案

核心错误原因

行主元交换逻辑仅同步置换了上三角矩阵U和置换矩阵P，未同步置换下三角矩阵L中已经计算完成的前j-1列元素。
2x2矩阵运算时外层循环j最大值为1，前j-1=0列无已计算内容，不会触发错误；3x3矩阵仅在j=2时需要置换L的前1列，若测试用例在该轮不需要行交换也不会暴露问题，矩阵规模增大后行交换概率提升，错误就会稳定复现。
存在冗余的维度计算逻辑：代码开头先执行n = size(A)得到二维维度向量，后续又重新执行[n, m] = size(A)获取行列数，虽然不会直接导致错误，但属于不严谨的写法，容易引发后续维度匹配问题。
可选优化：原代码中L矩阵系数计算行末尾未加分号，运行时会打印大量中间结果，建议补充分号关闭不必要的输出。

修复后完整代码

function [L, U, P] = luFactor(A)
[n, m] = size(A);
% 先校验矩阵是否为方阵
if n ~= m          
    error('Check dimensions of A')
end
U = A;          % U = 上三角矩阵
L = eye(n);     % L = 下三角矩阵
P = eye(n);     % P = 初始置换矩阵

for j = 1:n-1
    % 找当前列主元所在行
    [~, ind] = max(abs(U(j:n, j)));   
    pivot_row = ind + j - 1;
    % 置换U矩阵对应行
    [U(j, :), U(pivot_row, :)] = deal(U(pivot_row, :), U(j, :));
    % 置换P矩阵对应行
    [P(j, :), P(pivot_row, :)] = deal(P(pivot_row, :), P(j, :));
    % 新增：置换L矩阵已计算的前j-1列对应行
    if j > 1
        [L(j, 1:j-1), L(pivot_row, 1:j-1)] = deal(L(pivot_row, 1:j-1), L(j, 1:j-1));
    end

    % 消元计算
    for i = j:n-1
        L(i+1, j) = U(i+1, j) / U(j, j);
        U(i+1, :) = U(i+1, :) - L(i+1, j) * U(j, :);
    end
end
end