矩阵的LU分解，怎么修改啊？

redcandle

下面是自己写的一段，总是不对。求达人修改~
谢谢！！
function [L,U]=lu_self(A)
[n,m]=size(A);
L=eye(n,m);
U=A;
for k=1:n-1
    if U(k,k)==0
        stop;
    else      
        for i=k+1:n
            L(i,k)=U(i,k)/U(k,k);         
            for j=k+1:n
            U(i,j)=U(i,j)-L(i,k)*U(k,j);
            end           
        end
    end
end
end

qibbxxt

1.我建议你测试一个小例子，设置断点，检查每一步，这样对你自己也是提高
2.stop我很很在程序中用，是不是可以考虑用其他的关键字呢？

xuxu1457

你的U明显不对的，你程序中U的上三角上始终和A一样 没有动；可以细读下 help里面lu分解的lu子程序

1. function [L,U,P] = lu(A)
   
2. %LU  Lower-upper triangular factorization of GF array.
   
3. %   [L,U,P] = LU(X) performs Gaussian elimination on the square
   
4. %   GF matrix X and returns arguments as does the LU function
   
5. %   of MATLAB for double matrices.  Specifically:
   
6. %
   
7. %   [L,U] = LU(X) stores an upper triangular matrix in U and a
   
8. %   "psychologically lower triangular matrix" (i.e. a product
   
9. %   of lower triangular and permutation matrices) in L, so
   
10. %   that X = L*U.  X must be square.
    
11. %
    
12. %   [L,U,P] = LU(X) returns lower triangular matrix L, upper
    
13. %   triangular matrix U, and permutation matrix P so that
    
14. %   P*X = L*U.
    
15. %
    
16. %   See also LSOLVE, USOLVE, INV.
    
17. 
    
18. %    Copyright 1996-2002 The MathWorks, Inc.
    
19. %    $Revision: 1.5.4.1 $  $Date: 2002/11/13 16:53:08 $
    
20. 
    
21. global GF_TABLE_M GF_TABLE_PRIM_POLY GF_TABLE1 GF_TABLE2
    
22. 
    
23. n=length(A);
    
24. U=A;
    
25. if ~isequal(U.m,GF_TABLE_M) | ~isequal(U.prim_poly,GF_TABLE_PRIM_POLY)
    
26.     [GF_TABLE_M,GF_TABLE_PRIM_POLY,GF_TABLE1,GF_TABLE2] = gettables(U);
    
27. end
    
28. %  piv=A;
    
29. piv=zeros(1,n);
    
30. for k=1:n-1,
    
31.     %max_row = find(U(k:n,k)==max(U(k:n,k))) + k - 1;  % find maximums
    
32.     max_row = find(U.x(k:n,k)~=0) + k - 1;  % find non zero elements
    
33.     if isempty(max_row)
    
34.         error('matrix not invertible')
    
35.     end
    
36.     max_row = max_row(1);        % allow for repeated maximums
    
37.     if (max_row ~= k),
    
38.         U.x([k max_row],1:n) = U.x([max_row k],1:n);  % swap rows
    
39.     end;
    
40.     piv(k) = max_row;
    
41.     %U(k+1:n,k) = U(k+1:n,k)/U(k,k);
    
42.     Ukk_inv = gf_mex(U.x(k,k),U.x(k,k),U.m,'rdivide',...
    
43.         A.prim_poly,GF_TABLE1,GF_TABLE2);
    
44.     U.x(k+1:n,k) = gf_mex(U.x(k+1:n,k),Ukk_inv(ones(length(k+1:n),1)),U.m,'times',...
    
45.         A.prim_poly,GF_TABLE1,GF_TABLE2);
    
46.     % U(k+1:n,k+1:n) = U(k+1:n,k+1:n) - U(k+1:n,k)*U(k,k+1:n);
    
47.     temp = gf_mex(U.x(k+1:n,k),U.x(k,k+1:n),U.m,'mtimes',...
    
48.         A.prim_poly,GF_TABLE1,GF_TABLE2);
    
49.     U.x(k+1:n,k+1:n) = gf_mex(U.x(k+1:n,k+1:n),temp,U.m,'plus',...
    
50.         A.prim_poly,GF_TABLE1,GF_TABLE2);
    
51. end;
    
52. 
    
53. %    L=tril(U)+eye(U)-diag(diag(U));
    
54. L=A;
    
55. L.x=uint16(eye(n));
    
56. for i=2:n
    
57.     L.x(i,1:i-1)=U.x(i,1:i-1);
    
58. end;
    
59. U.x=triu(U.x);
    
60. 
    
61. % Create Permutation matrix:
    
62. P=A;
    
63. P.x=uint16(eye(n));    % start with identity
    
64. for k=1:n-1,
    
65.     P.x([k piv(k)],:)=P.x([piv(k) k],:);   % exchange row k with row p(k)
    
66. end;
    
67. 
    
68. % Adjust L according to the Permutation matrix
    
69. if(nargout==2)
    
70.     L = P'*L;
    
71. end

复制代码