%% Class 3: Gaussian elimination WITHOUT pivoting
% You need to download the m-file lunp.m and put it in the same folder as your other m-files.
%% Solving a linear system using Gaussian elimination WITHOUT pivoting
% We want to solve the linear system Ax=b

clearvars                         % Always clear variables at the beginning
format compact                    % omit blank lines between results
format short g                    % show each entry with 5 digits accuracy
A = [4 3 3; 8 7 9; -8 -1 3]
b = [1;4;2]

% We use Gaussian elimination without pivoting:
% elimination in column 1: For j=2,3,\ldots,n: subtract \ell_{j1}  times row 1 from row j
% elimination in column 2: For j=3,4,\ldots,n: subtract \ell_{j2}  times row 2 from row j
% ...
% This transforms the matrix A to an upper triangular matrix U.
% The lower triangular matrix L has 1's on the diagonal and the multipliers \ell_{jk} below the diagonal.

[L,U] = lunp(A)

% For a given right hand side vector b:
% Solve Ly=b by forward substitution: (Matlab's \ command checks whether the matrix is triangular and then performs substitution in the correct order)

y = L\b

% Solve Ux=y by back substitution to obtain the solution of the linear system Ax=b:

x = U\y

%% Problems with this algorithm:
%% (1) Algorithm can fail even for nonsingular matrix A when zero pivot is encountered:
% We want to solve the linear system \left\lbrack \begin{array}{cc}
% 0 & 1\\
% 2 & 0
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_1 \\
% x_2
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% 3\\
% 4
% \end{array}\right\rbrack. The algorithm stops since it encounters the pivot u_{11}=0.

A = [0 1; 2 0]
[L,U] = lunp(A)

%% (2) Algorithm can give very bad results due to tiny pivots:
% We want to solve the linear system \left\lbrack \begin{array}{ccc}
% \ldotp 1 & \ldotp 3 & 0\\
% \ldotp 3 & \ldotp 9 & 1\\
% 0 & 1 & 1
% \end{array}\right\rbrack \left\lbrack \begin{array}{c}
% x_1 \\
% x_2 \\
% x_3
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% 1\\
% 1\\
% 1
% \end{array}\right\rbrack.

A = [.1 .3 0; .3 .9 1; 0 1 1]
b = [1;1;1]

% The exact solution is

xex = A\b

% We now use Gaussian elimination without pivoting:

[L,U] = lunp(A)

% Note that we obtain a tiny pivot element u_{22}\approx 2.22\cdot 10^{-16}. In exact arithmetic we would have obtained a zero pivot.

y = L\b
xhat = U\y

% There is no division by zero, but Matlab alerts us that there may be loss of accuracy.
% Note that the computed solution \hat x is completely wrong: 100% relative error for x, and 58% relative error for \hat b=A\hat x.

relerr_x = norm(xhat-xex,Inf)/norm(xex,Inf)  % relative error of xhat
bhat = A*xhat;
relerr_b = norm(bhat-b,Inf)/norm(b)          % backward error: relative error of bhat = A*xhat