%% MCHI, cv1, file 1 (cv1_1) - Intro to MATLAB
%
% Author:       Martin Isoz
% Organisation: ICT Prague
% Date:         10. 09. 2014
%
% 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

%% Variables
% note: MATLAB = MATrix LABoratory, everything is a matrix

% Scalars (matrix 1x1)
a = 4;                                                                      %declare and define at a same time
b = 5;                                                                      %semicolon to surpress output

c = a+b;
disp(['writing output by disp, c = ' num2str(c)]);
fprintf(1,'writing output by fprintf, c = %d\n',c);

%% Linear equation solving
% Ax = b

A = floor(9*rand(5,5))+1;                                                   %create matrix 5x5 from numbers 1-10
A = A + diag(sum(A,2));                                                     %make it diagonally dominant
b = round(10*rand(5,1));                                                    %create matrix 5x1 from numbers 0-10

tic
x = A\b;
toc
disp('A = ');disp(A)
disp('b = ');disp(b)
disp('x = ');disp(x)
disp('Ax = ');disp(A*x);

% Gauss elimination - without pivotig
[m,n] = size(A);                                                            %size of matrix A
T     = [A b];                                                              %temporary matrix
disp('T0 = ');disp(T);
tic
for i = 2:m                                                                 %forward run
    T(i-1,i-1:end) = T(i-1,i-1:end)./T(i-1,i-1);
    T(i:m,:) = (T(i:m,:) - T(i:m,i-1)*T(i-1,:));
%     disp(['T' num2str(i-1) ' = ']);disp(T);
end
T(end,m:end) = T(end,m:end)./T(end,m);
for i = m-1:-1:1                                                            %backward run
   T(1:i,:) = T(1:i,:) - T(1:i,i+1)*T(i+1,:);
%    disp(['T' num2str(i) ' = ']);disp(T);
end
toc
disp('T_Final = ');disp(T);

% Jacobi method - iterative lin. eq. solver
maxiter = 20;                                                               %maximal number of iterations
eps     = 1e-2;                                                             %preset difference
delta   = eps+1;
xOld    = zeros(m,1);
xNew    = xOld;
LU      = diag(diag(A))-A;
iterMat = zeros(m,maxiter);

tic
k = 1;
while k <= maxiter && delta > eps
%     xNew = (LU*xOld + b)./diag(A);
    for i = 1:numel(xOld)
        xNew(i) = (b(i) - A(i,1:i-1)*xOld(1:i-1) - A(i,i+1:end)*xOld(i+1:end))/A(i,i);
    end
%     for i = 1:numel(xOld)
%         xNew(i) = (b(i) - A(i,1:i-1)*xNew(1:i-1) - A(i,i+1:end)*xOld(i+1:end))/A(i,i);
%     end
    delta = norm(A*xNew-b);
    xOld = xNew;
    iterMat(:,k) = xOld;
    k = k+1;
end
toc
disp('Jacobi iteration process')
disp(iterMat)