clc
clear

% simple matrix for which there is no LUdecomp
A = [0,2;3,1]

% LUdecomp(A)  % you can check that this will cause an error

[P,L,U] = PALUdecomp(A)

disp('** computed P^{-1}.L.U')
disp(P'*L*U)
disp('** check matrix A')
disp(A)

% slightly less trivial example where the error appears in step k=2
% the problem arises because B_{22} = B_{21} B_{12} / B_{11}
% so in the usual LU decomp then A^(1)_{22} is zero
% B = [ 1,2,3 ; 2,4,2 ; 1,3,2 ]

% LUdecomp(B)  % you can check that this will cause an error

% [P,L,U] = PALUdecomp(B)
% 
% 
% disp('** computed P^{-1}.L.U')
% disp(P'*L*U)
% disp('** check matrix B')
% disp(B)


% check the V2 decomp as well...
% 
% disp('** V2')
% [P2,L,U] = PALUdecompV2(B)

% check that LUdecomp of P.B gives the same answer
% 
% [L,U] = LUdecomp(P*B)

% bigger matrix...
%
% nn=800
% A = rand(nn);
% b = zeros(nn);
% 
% disp('** test PALU decomp')
% [P,L,U] = PALUdecomp(A);
% 
% disp('** top left block of L')
% disp( L( 1:5 , 1:5 ) )
% disp('** top left block of U')
% disp( U( 1:5 , 1:5 ) )
% checkLU = L*U;
% disp('** top left block of LU')
% disp( checkLU( 1:5 , 1:5 ) )
% disp('** top left block of PA')
% PA = P*A;
% disp( PA( 1:5 , 1:5 ) )

% % estimate "error" as sum of squares of elements of (LU-PA)
% estErr = trace( (L*U-P*A)*(L*U-P*A)' );
% disp(['** est error (roundoff) ',num2str(estErr)])
% disp(['** rms error per element ',num2str(sqrt(estErr)/nn/nn)])

% 
% disp('** test PALU decompV2')
% [P,L,U] = PALUdecompV2(A);
% 
% disp('** top left block of L')
% disp( L( 1:5 , 1:5 ) )
% disp('** top left block of U')
% disp( U( 1:5 , 1:5 ) )

% disp('** test LU decomp')
% [L,U] = LUdecomp(A);
% disp('** top left block of L')
% disp( L( 1:5 , 1:5 ) )
% disp('** top left block of U')
% disp( U( 1:5 , 1:5 ) )
% checkLU = L*U ;
% disp('** top left block of LU')
% disp( checkLU( 1:5 , 1:5 ) )

% % estimate "error" as sum of squares of elements of (LU-A)
% estErr = trace( (L*U-A)*(L*U-A)' );
% disp(['** est error (roundoff) ',num2str(estErr)])
% disp(['** rms error per element ',num2str(sqrt(estErr)/nn/nn)])


 
