% TestLinearSolve.m
% This file will be used to test the programs you wrote.  
% You must write 
%    (1)  myLUP which returns L, U, and P. 
%    (2)  myLUSolve which returns a solution to Ax = b. 

clear
clc

b = [0 -10 -39 -16]';




A1 = [6  -2  2    4; 
      12 -8  4   10; 
      3  -13 3    3; 
     -6  4   2 -18];
 

 [L,U,P] = myLUP(A1);
 disp('In the LUP factorization of A1')
 disp('P is')
 disp(P)
 disp('and L*U - P*A1 is')
 disp(L*U - P*A1)

 
 % Now I want to solve Ax = b or PAx = Pb or LUx = Pb  so 
 x =  myLUSolve(L,U,P*b);
 disp('the solution to A1 x = b is')
 disp(['x = [', num2str(x'), ']'])
 disp('and the residual is ')
 disp(['r = [', num2str((A1*x - b)'), ' ]'])
 
 
A2 = [6  -2  2    4; 
      12 -4  4   10; 
      3  -13 3    3; 
     -6  4   2 -18];
 
[L,U,P] = myLUP(A2);
disp('In the LUP factorization of A2, P is')
disp(P)
disp('and L*U - P*A2 is')
disp(L*U - P*A2)