% Script to test lutx.m for
% LU factorization of an n by n matrix A using 
% Gauss elimination with partial row pivoting.
% lutxt.m
% After permuting the rows of A according to
% the index vector p: A is factored as A = L*U
% Output:
% L is lower triangular with the main diagonal part = 1s.
% U is upper triangular 
% Pvec is a vector indexing the permuted rows of A.
% K. Ming Leung, 01/26/03
% K. Ming Leung, 02/21/04 (version 2), more general tests.

clear all; format short g; format compact;
% This is example 2.13 on p.69 of Heath:
AA = [1 2 2; 4 4 2; 4 6 4]; Bvec = [3; 6; 10]; % Exact solution: [-1; 3; -1]
%AA = [ 4 4 2; 4 6 4; 1 2 2]; Bvec = [ 6; 10; 3];
%A = [ 4 12 8 4; 1 7 18 9; 2 9 20 20; 3 11 15 14];  % Example on p.135 Fausett
%n = 1000; AA = rand(n,n);  Bvec = rand(n,1); % MATLAB's lu function will perform pivoting in general!AA = ones(2,2); Bvec = ones(2,1);
[n,n] = size(AA);
[LL,UU, Pvec] = lutx(AA);
%%display(LL); display(UU);
% Matlab's own LU factorization (with partial pivot)
%%[LMATLAB,UMATLAB] = lu(AA),

% solve the linear system using the LU factors
 Yvec = forwardSubstitution(LL,Bvec(Pvec),n);
 Xvec = backSubstitution(UU,Yvec,n),
 
 % check with Matlabs linear system solver
 Xmatlab = AA\Bvec,
 format long g;
 % display both solution vectors side-by-side
%% [Xvec Xmatlab],
 % display maximum absolute error
 max(abs(Xvec-Xmatlab)),
 