LU分解(2)

function [P, L, U] = zgaxpyplu(A)

%compute LU decomposition based on Gaxpy operation with pivoting

%aimed at improve the stability of LU decomposition

[m, n] = size(A);

if m ~= n

    error(‘current support square matrix only‘)

end

P = eye(n);

L = eye(n);

U = zeros(n);

for k = 1:n

    v = zeros(n, 1);

    if k == 1

        v(k:end) = A(k:end, 1);

    else

        U(1:k-1, k) = L(1:k-1, 1:k-1) \ A(1:k-1, k);

        v(k:n) = A(k:n, k) - L(k:n, 1:k-1)*U(1:k-1, k);

    end

    %find the largest element in v(k:end)

    [max_value, max_index] = max(v(k:end));

    max_index = max_index + k - 1;

    if max_index ~= k

        %exchange the max_index-th row and k-th row

        A([k max_index], k+1:n) = A([max_index k], k+1:n);

        %exchange the max_index-th row and k-th row in L

        if k > 1

            L([k max_index], 1:k-1) = L([max_index k], 1:k-1);

        end

        P([k max_index], :) = P([max_index k], :);

        v([k max_index]) = v([max_index k]);

    end

    if (abs(v(k)) > 1e-6) && (k < n)

        v(k+1:end) = v(k+1:end) / v(k);

        L(k+1:end, k) = v(k+1:end);

    end

    U(k, k) = v(k);

end

运行这段代码1500次，每次都使用随机生成的矩阵和Matlab的LU分解做对比，结果如下：

mean of my lu : 2.02547e-15
variance of my lu : 9.06349e-29
mean of matlab lu : 3.65765e-16
variance of matlab lu : 1.98304e-32

将每个算法和同类型的算法相比，使用gaxpy的算法都比之前的方法要好。但是现在有一个问题是使用pivot的算法居然比不使用pivot的算法没有足够显著的改进，当然这也证明了使用基于gaxpy运算可以提高LU分解的稳定性。

总结：LU分解是用于求解线性方程组最基本的矩阵分解方法，对于LU分解可以使用高斯消元法，也可以使用基于gaxpy的运算方法。但是从目前的结果来看，使用gaxpy运算可以极大的提高LU分解的稳定性。但是作为成熟的方法，使用lapack中对应的方法是最好的选择。