www.pudn.com > mylu.rar > mylu.m

function x=mylu(A,b) 
%Input  -A is an N*N matrix 
%       -b is N*1 matrix 
%Output -x is an N*1 matrix containing the solution to Ax=b 
[N,N]=size(A);x=zeros(N,1);y=zeros(N,1); 
U0=A;L0=eye(N);L0t=zeros(N);P0=eye(N);R=zeros(N,N); 
for p=1:N-1 
    %Find the pivot row for colum p 
    [max1,j]=max(abs(U0(p:N,p))); 
    %Interchange row p and j 
    R(p,:)=1:N;P1=eye(N); 
    c=P1(p,:);P1(p,:)=P1(j+p-1,:);P1(j+p-1,:)=c; 
    d=R(p,p);R(p,p)=R(p,j+p-1);R(p,j+p-1)=d; 
    U0=P1*U0;P0=P1*P0; 
    if A(p,p)==0  
        'A is singular.No unique solution'  
        break 
    end 
    L1=zeros(N); 
    for k=p+1:N   
        mult=U0(k,p)/U0(p,p);U0(k,p:N)=U0(k,p:N)-mult*U0(p,p:N);L1(k,p)=-mult; 
    end 
    L0t=L0t+L1; 
end 
L0t=L0t+eye(N);R(N,:)=1:N; 
for p=1:N-1 
    L1=eye(N);L1(:,p)=L0t(:,p); 
    for k=p:N-1 P2=eye(N); 
        for n=1:N P1(:,n)=P2(:,R(k+1,n));end  
        L1=P1*L1*(eye(N)/P1); 
    end 
    L0=L0*(eye(N)/L1); 
end 
%Solve for y 
b=P0*b;y(1)=b(1);for k=2:N y(k)=b(k)-L0(k,1:k-1)*y(1:k-1);end 
%Solve for x 
x(N)=y(N)/U0(N,N);for k=N-1:-1:1 x(k)=(y(k)-U0(k,k+1:N)*x(k+1:N))/U0(k,k);end