# 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

Mentor
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
Code:
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
Code:
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
Code:
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
Code:
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