% PPivotCompare Partial Pivoting Comparison
%   [x0,x1,x2] = PPivotCompare( alpha, flag )
%     Performs a comparison between numerical solutions to
%               alpha*x + y = 1
%                     x + y = 2
%     x0 is the "exact solution" computed by an explicit formula
%     x1 is the solution found by LU factorization with no pivoting
%     x2 is the solution found by LU factorization with partial pivoting
%
%     with flag==1 the results of each solution, and their residuals, are displayed
%     for all other values of flag, there is no intermediate output

% Written by Douglas Meade
% Created on 11 Sep 2005
% Revised on 13 Sep 2005 --- replaced forward and backward with \

function [x0,x1,x2] = PPivotCompare(alpha,flag)

format long g

A = [ alpha 1; 1 1 ];
b = [ 1; 2 ];

x0 = [1/(1-alpha); (1-2*alpha)/(1-alpha)];

% L-U Factorization without Partial Pivoting

%[L,U] = lu( A )
M1 = [1 0;-A(2,1)/A(1,1) 1];
U1 = M1*A;
L1 = [1 0; A(2,1)/A(1,1) 1];

%y1 = forward( L1, b );
%x1 = backward( U1, y1 );
y1 = L1\b;
x1 = U1\y1;

% L-U Factorization with Partial Pivoting

%[P,L,U] = lu( A )
if abs(A(2,1))>abs(A(1,1)),
   P2 = [0 1;1 0];
else
   P2 = eye(2);
end
A2 = P2*A;
M2 = [1 0;-A2(2,1)/A2(1,1) 1];
U2 = M2*A2;
L2 = [1 0; A2(2,1)/A2(1,1) 1];

c = P2*b;
%y2 = forward( L2, c );
%x2 = backward( U2, y2 );
y2 = L2\c;
x2 = U2\y2;

if flag==1,
   A2
   fprintf('alpha=%g, epsilon=%g\n', alpha, eps )
   fprintf('Exact Solution:\n')
   disp([x0 A*x0-b])
   fprintf('GE w/No Pivoting:\n')
   disp([x1 A*x1-b])
   fprintf('Partial Pivoting:\n')
   disp([x2 A*x2-b])
end