Question

Need a MATLAB code for LU factorization(both partial and complete pivoting) of 10 random matrices of order 5x5.

Please do not comment like I don't have computer or I don't know matlab. If you please answer otherwise you can skip.

Answer 1

%%% Matlab code LU decomposition with complete pivoting

clc;
close all;
clear all;
tol=10^(-8);
A=rand(5,5);
[n m] = size(A);
% pivot vectors
p = (1:n)';
q = (1:m)';
% temp storage
rt = zeros(m,1); % row temp
ct = zeros(n,1); % col temp
t = 0; % scalar temp
for k = 1:min(n-1,m)
% determine pivot
[cv ri] = max(abs(A(k:n,k:m)));
[rv ci] = max(cv);
rp = ri(ci)+k-1;
cp = ci+k-1;
  
% swap row
t = p(k);
p(k) = p(rp);
p(rp) = t;
rt = A(k,:);
A(k,:) = A(rp,:);
A(rp,:) = rt;
  
% swap col
t = q(k);
q(k) = q(cp);
q(cp) = t;
ct = A(:,k);
A(:,k) = A(:,cp);
A(:,cp) = ct;
  
if abs(A(k,k)) >= tol
rows = (k+1):n;
cols = (k+1):m;
A(rows,k) = A(rows,k)/A(k,k);
A(rows,cols) = A(rows,cols)-A(rows,k)*A(k,cols);
else
% stop factorization if pivot is too small
break
end
end
% final column swap if m > n
if m > n
% determine col pivot
[cv ci] = max(abs(A(n,n:m)));
cp = ci+n-1;
  
% swap col
t = q(n);
q(n) = q(cp);
q(cp) = t;
ct = A(:,n);
A(:,n) = A(:,cp);
A(:,cp) = ct;
end
% produce L and U matrices
% these are sparse if L and U are sparse
l = min(n,m);
L = tril(A(1:n,1:l),-1) + speye(n,l);
U = triu(A(1:l,1:m));
disp('Matrix is : ');
A
disp('L matrix is :');
L
disp('U matrix is :');
U

OUTPUT:

Matrix is :

A =

0.9340 0.9172 0.5678 0.4733 0.3371
0.0578 0.7778 0.7244 0.1693 0.5494
0.5683 0.0823 0.3714 -0.0318 0.2326
0.8342 -0.2769 0.2891 0.2773 -0.1844
0.0812 0.3563 -0.0495 0.9003 0.0844

L matrix is :

L =

1.0000 0 0 0 0
0.0578 1.0000 0 0 0
0.5683 0.0823 1.0000 0 0
0.8342 -0.2769 0.2891 1.0000 0
0.0812 0.3563 -0.0495 0.9003 1.0000

U matrix is :

U =

0.9340 0.9172 0.5678 0.4733 0.3371
0 0.7778 0.7244 0.1693 0.5494
0 0 0.3714 -0.0318 0.2326
0 0 0 0.2773 -0.1844
0 0 0 0 0.0844

%%%% Partial pivoting

clc;
close all;
clear all;
A=rand(5,5);
[m, n] = size(A);
L=eye(n);
P=eye(n);
U=A;
for k=1:m-1
pivot=max(abs(U(k:m,k)));
for j=k:m
if(abs(U(j,k))==pivot)
ind=j;
break;
end
end
U([k,ind],k:m)=U([ind,k],k:m);
L([k,ind],1:k-1)=L([ind,k],1:k-1);
P([k,ind],:)=P([ind,k],:);
for j=k+1:m
L(j,k)=U(j,k)/U(k,k);
U(j,k:m)=U(j,k:m)-L(j,k)*U(k,k:m);
end
  
end

disp('Matrix is : ');
A
disp('L matrix is :');
L
disp('U matrix is :');
U

OUTPUT:

Matrix is :

A =

0.4173 0.4893 0.7803 0.1320 0.2348
0.0497 0.3377 0.3897 0.9421 0.3532
0.9027 0.9001 0.2417 0.9561 0.8212
0.9448 0.3692 0.4039 0.5752 0.0154
0.4909 0.1112 0.0965 0.0598 0.0430

L matrix is :

L =

1.0000 0 0 0 0
0.9555 1.0000 0 0 0
0.4417 0.5960 1.0000 0 0
0.0526 0.5817 0.6577 1.0000 0
0.5195 -0.1474 -0.1958 -0.2738 1.0000

U matrix is :

U =

0.9448 0.3692 0.4039 0.5752 0.0154
0 0.5472 -0.1442 0.4065 0.8065
0 0 0.6878 -0.3644 -0.2527
0 0 0.0000 0.9150 0.0495
0 0 0 0 0.1179

>>