Matlab debugging for LU decomposition

[joehardy]

these are my mfiles for Lu decomposition with partial pivoting and substitution

1.LUdecomp.m
function [x]=LUdecomp(a,b,n,tol,er,x)
for i=1:n
o(i)=0;
s(i)=0;
end
er=0;
decompose(a,n,tol,o,s,er)
if er~=-1
substitute(a,o,n,b,x);
end
end

2.decompose.m
function [a]=decompose(a,n,tol,o,s,er)
for i=1:n
o(i)=i;
s(i)=abs(a(i,1));
for j=2:n
if abs(a(i,j))>s(i)
s(i)=abs(a(i,j));
end
end
end
for k=1:n-1
pivot(a,o,s,n,k)
if abs(a(o(k),k)/s(o(k)))<tol
er=-1;
abs(a(o(k),k)/s(o(k)));
break
end
for i=k+1:n
factor=a(o(i),k)/a(o(k),k);
a(o(i),k)=factor;
for j=k+1:n
a(o(i),j)=a(o(i),j)-factor*a(o(k),j);
end
end
end
if abs(a(o(k),k)/s(o(k)))<tol
er=-1;
abs(a(o(k),k)/s(o(k)));
end
end

3.pivot.m
function [o]=pivot(a,o,s,n,k)
p=k;
big=abs(a(o(k),k)/s(o(k)));
for i=k+1:n
dummy=abs(a(o(i),k)/s(o(i)));
if dummy>big
big=dummy;
p=i;
end
end
dummy=o(p);
o(p)=o(k);
o(k)=dummy;
end

4.substitute.m

function [x]=substitute(a,o,n,b,x)
for i=2:n
sum=b(o(i));
for j=1:i-1
sum=sum-a(o(i),j)*b(o(j));
end
b(o(i))=sum;
end
x(n)=b(o(n))/a(o(n),n);
for i=n-1:-1:1
sum=0;
for j=i+1:n
sum=sum+a(o(i),j)*x(j);
end
x(i)=(b(o(i))-sum)/a(o(i),i);
end
end
thank you

 

[Mark44]

What's your question? From the thread title, you seem to be looking for help in debugging this code. Can you narrow down which function is not working correctly?

Also, it's a good idea to surround your code with [ code] and [ /code] tags (omit the leading space.

these are my mfiles for Lu decomposition with partial pivoting and substitution

1.LUdecomp.m

function [x]=LUdecomp(a,b,n,tol,er,x)
for i=1:n
    o(i)=0;
    s(i)=0;
end
er=0;
decompose(a,n,tol,o,s,er)
if er~=-1
    substitute(a,o,n,b,x);
end
end

2.decompose.m

function [a]=decompose(a,n,tol,o,s,er)
    for i=1:n
        o(i)=i;
        s(i)=abs(a(i,1));
        for j=2:n
            if abs(a(i,j))>s(i)
            s(i)=abs(a(i,j)); 
            end
        end
    end
for k=1:n-1
    pivot(a,o,s,n,k)
    if abs(a(o(k),k)/s(o(k)))<tol
        er=-1;
        abs(a(o(k),k)/s(o(k)));
    break
    end
    for i=k+1:n
        factor=a(o(i),k)/a(o(k),k);
        a(o(i),k)=factor;
        for j=k+1:n
            a(o(i),j)=a(o(i),j)-factor*a(o(k),j);
        end
    end
end
if abs(a(o(k),k)/s(o(k)))<tol
        er=-1;
        abs(a(o(k),k)/s(o(k)));
end
end

3.pivot.m

function [o]=pivot(a,o,s,n,k)
p=k;
big=abs(a(o(k),k)/s(o(k)));
for i=k+1:n
    dummy=abs(a(o(i),k)/s(o(i)));
    if dummy>big
        big=dummy;
        p=i;
    end
end
dummy=o(p);
o(p)=o(k);
o(k)=dummy;
end

4.substitute.m

function [x]=substitute(a,o,n,b,x)
for i=2:n
    sum=b(o(i));
    for j=1:i-1
        sum=sum-a(o(i),j)*b(o(j));
    end
    b(o(i))=sum;
end
x(n)=b(o(n))/a(o(n),n);
for i=n-1:-1:1
    sum=0;
    for j=i+1:n
        sum=sum+a(o(i),j)*x(j);
    end
    x(i)=(b(o(i))-sum)/a(o(i),i);
end
end

thank you