LU

A = [2,1,-2;3,2,2;5,4,3]; b = [1;7;12];
A0 = A; b0 = b;
n = size(A,1);
I = eye(n);
k = 0;  %Initialise phase counter

k = k+1; % Reset/update phase counter
% Locate position of element largest modulus 
% in the pivot element
[~, r] = max(abs(A(k:n,k)));
% Create pivot matrix to interchange 
% rows r & k
P1=I; P1([k,r],:) = P1([r,k],:);
% Interchange rows r & k
A=P1*A;
% Create the Gaussian Transformation matrix
i = k+1:n; M1=I; M1(i,k)=-A(i,k)/A(k,k)
% Zero out elements in below pivot element
A = M1*A

k = k+1;
% Locate position of element largest modulus 
% in the pivot element
[~, r] = max(abs(A(k:n,k)));
r = r + k -1; % Correct position in pivot column
% Create pivot matrix to interchange 
% rows r & k
P2=I; P2([k,r],:) = P2([r,k],:); 
% Interchange rows r & k of A
A=P2*A;
% Interchange rows/columns r & k of M1
M1 = P2*M1*P2
% Create the Gaussian Transformation matrix
i = k+1:n; M2=I; M2(i,k)=-A(i,k)/A(k,k)
% Zero out elements in below pivot element
A = M2*A


P = P2*P1; % Permutation matrix
L = I/(M2*M1); % Lower triangular matrix
U = A; % Upper triangular matrix

[Lm, Um, Pm] = lu(A0);

P, Pm
L, Lm
U, Um

y = L\(P*b)


x = U\y


xm = A0\b0


A = [3,1,4,-1;2,-2,-1,2;5,7,14,8;1,3,2,4]; b = [7;1;20;-4]; 
A0 = A; b0 = b;
n = size(A,1);
I = eye(n); 
k = 0; %Initialise phase counter



k = k+1; % Reset/update phase counter
% Locate position of element largest modulus 
% in the pivot element
[~, r] = max(abs(A(k:n,k)));
r = r+k-1; %Adjust to correct location
% Create pivot matrix to interchange 
% rows r & k
P1=I; P1([k,r],:) = P1([r,k],:); 
% Interchange rows r & k of A
A=P1*A;
% Create the Gaussian Transformation matrix
i = k+1:n; M1=I; M1(i,k)=-A(i,k)/A(k,k)
% Zero out elements in below pivot element
A = M1*A


k = k+1; % Reset/update phase counter
% Locate position of element largest modulus 
% in the pivot element
[~, r] = max(abs(A(k:n,k)));
r = r + k -1; %Adjust to correct location
% Create pivot matrix to interchange 
% rows r & k
P2=I; P2([k,r],:) = P2([r,k],:); A=P2*A;
% Interchange rows/columns r & k of M1
M1 = P2*M1*P2
% Create the Gaussian Transformation matrix
i = k+1:n; M2=I; M2(i,k)=-A(i,k)/A(k,k)
% Zero out elements in below pivot element
A = M2*A



k = k+1;% Reset/update phase counter
% Locate position of element largest modulus 
% in the pivot element
[~, r] = max(abs(A(k:n,k)));
r = r + k -1;%Adjust to correct location
% Create pivot matrix to interchange 
% rows r & k
P3=I; P3([k,r],:) = P3([r,k],:); A=P3*A;
% Interchange rows/columns r & k of M1, M2
M1 = P3*M1*P3, M2 = P3*M2*P3
i = k+1:n; M3=I; M3(i,k)=-A(i,k)/A(k,k)
% Zero out elements in below pivot element
A = M3*A


P = P3*P2*P1; % Permutation matrix
L = I/(M3*M2*M1); % Lower triangular matrix
U = A; % Upper triangular matrix

disp('Permution matrix, P')
disp(P)
disp('Lower triangular matrix, L')
disp(L)
disp('Upper triangular matrix, U')
disp(U)


A = [2,4,3,2;3,6,5,2;2,5,2,-3;4,5,14,14];
n = size(A,1); I = eye(n); 
k = 0; % Initialise phase counter
L = zeros(n); % Initialise Lower triangular matrix



k = k+1;
[~, r] = max(abs(A(k:n,k)));
r = r+k-1;
P1=I; P1([k,r],:) = P1([r,k],:); A=P1*A;
i = k+1:n; M1=I; M1(i,k)=-A(i,k)/A(k,k)
% Update lower triangular matrix at phase 1
L(i,k)=A(i,k)/A(k,k)
A = M1*A



k = k+1;
[~, r] = max(abs(A(k:n,k)));
r = r + k -1;
P2=I; P2([k,r],:) = P2([r,k],:); A=P2*A;
%M1 = P2*M1*P2
i = k+1:n; M2=I; M2(i,k)=-A(i,k)/A(k,k)
% Update lower triangular matrix at phase 2
L = P2*L; L(i,k)=A(i,k)/A(k,k)
A = M2*A



k = k+1;
[~, r] = max(abs(A(k:n,k)));
r = r + k -1;
P3=I; P3([k,r],:) = P3([r,k],:); A=P3*A;
%M1 = P3*M1*P3, M2 = P3*M2*P3
i = k+1:n; M3=I; M3(i,k)=-A(i,k)/A(k,k)
% Update lower triangular matrix at phase 2
L = P3*L; L(i,k)=A(i,k)/A(k,k)
A = M3*A


L = I + L;
U = A;
P = P3*P2*P1;

disp('Permution matrix, P')
disp(P)
disp('Lower triangular matrix, L')
disp(L)
disp('Upper triangular matrix, U')
disp(U)


A = [2,4,3,2;3,6,5,2;2,5,2,-3;4,5,14,14];
U = A; n = size(A,1); L = zeros(n);
I = eye(n);  P = I;

for k = 1:n-1
    [~, r] = max(abs(U(k:n,k)));
    r = r + k -1;
    % Permutation matrix to swap rows r & k
    Pk=I; Pk([k,r],:) = Pk([r,k],:); 
    % Swap rows r & k
    U=Pk*U
    % Create kth Gaussian Transformation matrix
    i = k+1:n; m = U(i,k)/U(k,k);
    Mk=I; Mk(i,k)=-m;
    % Update elements of L
    L = Pk*L; L(i,k)= m %U(i,k)/U(k,k)
    % Update permutation matrix
    P = Pk*P;
    % Zero out elements below pivot element
    U = Mk*U
end
L = I + L;
disp('Permution matrix, P')
disp(P)
disp('Lower triangular matrix, L')
disp(L)
disp('Upper triangular matrix, U')
disp(U)


[Lm, Um, Pm] = lu(A0);




P, Pm
L, Lm
U, Um


y = L\(P*b)


x = U\y


xm = A0\b0

%A = [2,4,3,2;3,6,5,2;2,5,2,-3;4,5,14,14];
A = [pi,-exp(1),sqrt(2),-sqrt(3);pi^2,exp(1),-exp(1)^2,3/7; sqrt(5),-sqrt(6),1,-sqrt(2);pi^3,exp(1)^2,-sqrt(7),1/9]
%b = [sqrt(11); 0; pi; sqrt(2)];
U = A; n = size(A,1); L = zeros(n);
I = eye(n);  P = I;

for k = 1:n-1
    % Create kth Gaussian Transformation matrix
    i = k+1:n; m = U(i,k)/U(k,k);
    Mk=I; Mk(i,k)=-m;
    % Update elements of L
    L(i,k)= m %U(i,k)/U(k,k)
    % Zero out elements below pivot element
    U = Mk*U
end
L = I + L;
disp('Permution matrix, P')
disp(P)
disp('Lower triangular matrix, L')
disp(L)
disp('Upper triangular matrix, U')
disp(U)
disp('Solution')
disp((U\(L\(P*b)))')