%% MCHI, cv1, file 1 (cv1_1) - Intro to MATLAB
%
% Author:       Martin Isoz
% Organisation: ICT Prague
% Date:         20. 09. 2017
%
% License: This code is published under MIT License, please do not abuse
% it.
%
%
% MATLAB is obect extension of C++ (well, was at least) => 
%  => EVERYTHING is an OBJECT

%% Clear environment
close all                                                                   %closes all figures and visible windows (w/o main, editor...)
clear all                                                                   %clears all the variables in base workspace (base vs caller workspace)
clc                                                                         %clears the command window outputs

%% Set the script variables
systemSize = 5;
numOfRHS   = 1;

%% LU decomposition (without pivoting)
% A = LU; Ax = b -> LUx = b -> Ly = b -> Ux = y
% suitable for repeated solution of a system with the same system matrix
% but different right hand sides

A = floor(9*rand(systemSize,systemSize))+1;                                 %create matrix systemSizexsystemSize from numbers 1-10
A = A + diag(sum(A,2));                                                     %make it diagonally dominant

% LU decomposition of A (a rebranded Gaussian elimination)
[m,n] = size(A);                                                            %size of matrix A
T     = A;                                                                  %temporary matrix
L     = eye(m,n);                                                           %prepare the output matrix
% disp('T0 = ');disp(T);
tic
for i = 2:m+1                                                               %forward run of Gaussian elimination
    Ln              = eye(m,n);                                             %initialize the matrix for storing current L
    Ln(i:end,i-1) = -T(i:end,i-1)./T(i-1,i-1);
    L = L-tril(Ln,-1);                                                      %I do not want to repeatedly add the diagonal of Ln
    T = Ln*T;
%     disp(['L' num2str(i-1) ' = ']);disp(Ln);
%     disp(['T' num2str(i-1) ' = ']);disp(T);
end
U = T;

% Note: L = sum(-tril(Ln,-1),n=2,m+1) + eye(size(A)) corresponds to the
%       product of inverses of Ln (can be proved)

% display the results
disp('A  = ');disp(A);
disp('U  = ');disp(U);
disp('L  = ');disp(L);

% control output
disp('A - L*U');disp(A-L*U);

%% Solution of a SLAE using the LU decomposition of the system matrix
% solve the system Ax = b -> LUx = b for different b

for cycleCount = 1:numOfRHS

    % prepare the current RHS (right hand side)
    b = round(10*rand(systemSize,1));                                       %create matrix systemSizex1 from numbers 0-10
    
    % solution of the system Ly = b by forward substitution
    T = [L b];
    for i = 2:m                                                             %forward run
        T(i:m,:) = (T(i:m,:) - T(i:m,i-1)*T(i-1,:));
    end
    y = T(:,m+1);                                                           %save the intermediate solution
    disp('L*y-b = ');disp(L*y-b);
    
    % solution of the system Ux = y by backward substitution
    T = [U y];
    T(end,m:end) = T(end,m:end)./T(end,m);
    for i = m-1:-1:1                                                        %backward run
        T(i,i:end) = T(i,i:end)./T(i,i);
        T(1:i,:) = T(1:i,:) - T(1:i,i+1)*T(i+1,:);
    end
    x = T(:,m+1);                                                           %save the final solution
    
    disp('U*x-y = ');disp(U*x-y);
    disp('A*x-b = ');disp(A*x-b);
end