refac

.gitignore

@@ -1,2 +1,5 @@
+# LaTeX
+build
+
 # MATLAB
 *.asv

01_Direct-Methods-for-Solving-Linear-Systems/Code/main.m

@@ -23,8 +23,6 @@ b = P*b;
 y = Alg3(L, b);     % Forward substitution
 x = Q*Alg4(U, y)	% Backward substitution
 
-
-
 %% Problem 4
 
 A = [0.835, 0.667;

01_Direct-Methods-for-Solving-Linear-Systems/README.md

@@ -0,0 +1,10 @@
+## Lab 1 - Direct Methods for Solving Linear Systems
+
+
+
+### Bibliography
+
+References in the code comments are represented by `[number]`:
+
+  - 1
+  - 2 – Golub, van Loan, Matrix Computations, 3rd edition

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg1.m

@@ -0,0 +1,26 @@
+function A = Alg1(A)
+% Algorithm 1: Outer Product Gaussian Elimination
+% Performs a gaussian eliminaion on a square matrix A.
+
+[m, n] = size(A);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+for k = 1:n-1
+    if det(A(1:k, 1:k)) < eps
+        error('Matrix is not nonsingular!')
+    end
+end
+
+% The following algorithm is based on the Algrotihm 3.2.1 from [2].
+
+for k = 1 : m-1
+    rows = k + 1 : m;
+    A(rows, k) = A(rows, k)/A(k, k);
+    A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows);
+end
+ 
+end
+

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg10.m

@@ -0,0 +1,24 @@
+function A = Alg10(A)
+% Algorithm 10: Cholesky (Banachiewicz) factorization.
+% Matrix A has to be symmetric and positive-definite:
+% all d(i, i) > 0 for A = LDL^T.
+% We are using Outer product Version (Golub, Load, Alg. 4.2.2), which
+% computes lower triangular G, such that A = G*G^T
+
+[m, n] = size(A);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+for k = 1:m
+   A(k,k) = sqrt(A(k,k));
+   A(k+1:n, k) = A(k+1:n, k)/A(k,k);
+   for j = k+1:n
+       A(j:n, j) = A(j:n, j) - A(j:n, k)*A(j, k);
+   end
+end
+
+A = tril(A);
+
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg11.m

@@ -0,0 +1,22 @@
+function [Q R] = Alg11(A)
+% Algorithm 11: QR factorization via Householder transformation.
+
+[m, n] = size(A);
+
+if ~ (m>=n)
+   error("m has to be greater or equal n!") 
+end
+
+for j = 1:n
+    [v, beta] = householder(A(j:m, j))
+    A(j:m, j:n) = (eye(m-j+1)-beta*(v*v.'))*A(j:m, j:n)
+    if j < m
+        A(j+1:m, j) = v(2:m-j+1)
+    end
+end
+
+R = triu(A(1:n, 1:n))
+
+Q = tril(A, -1)
+
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg12.m

@@ -0,0 +1,21 @@
+function [Q, R] = Alg12(A)
+% Algorithm 11 - The  QR  algorithm  by  the  Householder  transformation
+% [Q, R] = Alg12(A)
+% Based on [GVL96]: Algorithm 5.2.4
+[m, n] = size(A);
+
+if ~ (m >= n)
+   error("m has to be greater or equal n!") 
+end
+
+Q = eye(m, m);
+for j = 1:n
+    for i = m:-1:j+1
+        [c, s] = Givens(A(i-1, j), A(i, j));
+        A(i-1:i, j:n) = [c s; -s c]'*A(i-1:i, j:n);
+        Q(:, i-1:i) = Q(:, i-1:i)*[c s ; -s c];
+    end
+end
+R = A;
+  
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg13.m

@@ -0,0 +1,19 @@
+function [Q, R] = Alg13(A)
+% Algorithm 13: The QR algorithm by the Gram-Schmidt ortogonalization
+% [Q, R] = Alg13(A)
+% Based on [GVL96]: Section 5.2.8
+
+[m,n] = size(A);
+
+if ~ (m >= n)
+   error("m has to be greater or equal n!") 
+end
+
+R = zeros(n, n);
+Q = zeros(m, n);
+for k = 1:n
+    R(k, k) = norm(A(:, k));
+    Q(:, k) = A(:, k)/R(k, k);
+    R(k, k+1:n) = Q(:, k)'*A(:, k+1:n);
+    A(:, k+1:n) = A(:, k+1:n) - Q(:, k)*R(k, k+1:n);
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg2.m

@@ -0,0 +1,47 @@
+function [P, Q, L, U] = Alg2(A)
+% Algorithm 2: Gaussian Elimination with Complete Pivoting.
+% [P, Q, L, U] = Alg2_gaussian_elimination_with_complete_pivoting(A)
+% computes the complete pivoting factorization PAQ = LU.
+
+[m, n] = size(A);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+% p and q are permutation vectors - respectively rows and columns
+p = 1:m;
+q = 1:m;
+
+% The following algorithm is based on the Algrotihm 3.4.2 from [2].
+
+for k = 1 : m-1
+    i = k:m;
+    j = k:m;
+    % Find the maximum entry to be the next pivot
+    [max_val, rows_of_max_in_col] = max(abs(A(i, j)));
+    [~, max_col] = max(max_val);
+    max_row = rows_of_max_in_col(max_col);
+    % Assign value of mu and lambda in respect to the main matrix A
+    [mi, lm] = deal(max_row+k-1, max_col+k-1);
+    % Interchange the rows and columns of matrix A...
+    A([k mi], 1:m) = deal(A([mi k], 1:m));
+    A(1:m, [k lm]) = deal(A(1:m, [lm k]));
+    % ...and respective permutation vectors entries.
+    p([k, mi]) = p([mi, k]);
+    q([k, lm]) = q([lm, k]);
+    % Perform Gaussian elimination with the greatest pivot
+    if A(k, k) ~= 0
+        rows = k+1 : m;
+        A(rows, k) = A(rows, k)/A(k, k);
+        A(rows, rows) = A(rows, rows) - A(rows, k) * A(k, rows);
+    end
+end
+
+I = eye(m);
+U = triu(A);
+L = tril(A, -1) + I;
+P = I(p, :);
+Q = I(:, q);
+
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg3.m

@@ -0,0 +1,29 @@
+function b = Al3(L, b)
+% Algorithm 3: Forward Substitution
+% b = Alg3_forward_substitution(L, b) overwrites b with the solution to
+% Lx = b.
+
+[m, n] = size(L);
+
+if m ~= n
+    error('Matrix is not square!')
+end
+
+if length(b) ~= m
+    error('Vector b has wrong length!')
+end
+
+if det(L) < eps
+    error("Matrix is not nonsingular!")
+end
+
+% The following algorithm is based on the Algrotihm 3.1.1 from [2].
+
+% b(m, :) so that matrices are also accepted
+b(1, :) = b(1, :)/L(1,1);
+for i = 2:m
+    b(i, :) = (b(i, :) - L(i, 1:i-1)*b(1:i-1, :))/L(i, i);
+end
+
+end
+

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg4.m

@@ -0,0 +1,33 @@
+function b = Alg4(U,b)
+% Argorithm 4: Back Substitution 
+% b = Alg4_back_substitution(U,b) returns vetor b with solution to the
+% Ux = b.
+
+[m, n] = size(U);
+
+if U ~= triu(U)
+    error('Matrix U is not upper triangular!')
+end
+
+if m ~= n
+    error('Matrix is not squared!')
+end
+
+if length(b) ~= m
+    error('Vector b has wrong length!')
+end
+
+% if det(U) < eps
+%     error('Matrix is not nonsingular!')
+% end
+
+
+% The following algorithm is based on the Algrotihm 3.1.2 from [2].
+
+% b(m, :) so that matrices are also accepted
+b(m, :) = b(m, :)/U(m, m);
+for i = m-1:-1:1
+    b(i, :) = (b(i, :) - U(i, i+1 : m)*b(i+1 : m, :))/U(i, i);
+end
+
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg5.m

@@ -0,0 +1,21 @@
+function A = Alg5(A)
+% Algorithm 5: Gauss-Jordan Elimination
+% A = Alg5_gauss_jordan_elimination(A) performs Gauss-Jordan elimination
+% on an augmented matrix A.
+
+[m, ~] = size(A);
+
+for k = 1 : m
+    
+    row = A(k, :);
+    row = row/row(k);
+    A(k, :) = row;
+    
+    for i = 1 : m
+        if i ~= k && A(i, k) ~= 0
+            A(i, :) = A(i, :)-(A(i, k))*row;
+        end
+    end
+end
+
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg6_RREF.m

@@ -0,0 +1,46 @@
+function A = Alg6_RREF(A)
+% Algorithm 6: Reduced Row Echelon Form (RREF)
+% A = Alg6_RREF(A) returns RREF of matrix A. 
+
+[m, n] = size(A);
+
+j = 0;
+
+for k = 1 : m
+    j = j + 1;
+    if j > n
+        break
+    end
+    % We want the left-most coefficient to be 1 (pivot)
+    row = A(k, :);
+    if row(k) == 0
+        j = j + 1;
+    end
+    row = row/row(j);
+    A(k, :) = row;
+    
+    for i = 1 : m
+        if i ~= k
+            A(i, :) = A(i, :)-(A(i, j))*row;
+        end
+    end
+    
+    for i = k + 1 : m
+       A(i:end, k+1:end); % Partial matrix (in which we are looking for non-zero pivots)
+       A(i:end, k+1); % Left-most column
+       if ~any(A(i:end, k+1)) % If the left-most column has only zeros check the next one
+           k = k + 1;
+       end
+       A(i:end, k+1:end);
+       if A(i, k+1) == 0
+           non_zero_row = find(A(i:end,k+1), 1);
+           if isempty(non_zero_row)
+               continue
+           end
+           A([i, i+non_zero_row-1], :) = deal(A([i+non_zero_row-1, i], :));
+       end
+    end
+end
+
+end
+

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg7.m

@@ -0,0 +1,35 @@
+function [L, U] = Alg7(A)
+% LU Factorization without pivoting
+% [L, U] = Alg7(A) decomposes matrix A using Crout's algorithm into
+% U - upper triangular matrix and L - unit lower triangular matrix
+% such that A = LU.
+
+[m, n] = size(A);
+
+% This solution is mathematically correct, however computationaly
+% inefficient.
+%
+% A = Alg1_outer_product_gaussian_elimination(A);
+% 
+% U = triu(A);
+% L = tril(A, -1) + eye(m);
+
+% Instead, we should use the Crout's algorithm
+
+L = zeros(m, n);
+U = eye(m, n);
+
+for k = 1:n
+    
+    for j = k : n
+        U(k, j) = A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j));
+    end
+    
+    for i = k:n
+        L(i, k) = (A(i, k) - dot(L(i,1:k-1), U(1:k-1, k))) / U(k,k);
+    end
+    
+    
+end
+
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/algorithms/Alg8.m

@@ -0,0 +1,34 @@
+function [L, U, P] = Alg8(A)
+% LU Factorization with partial pivoting
+% [L, U, P] = Alg8(A) decomposes matrix A using Crout's algorithm into
+% U - upper triangular matrix and L - unit lower triangular matrix
+% using partial pivoting, interchanging A's rows, such that AP = LU.
+
+[m, n] = size(A);
+
+L = zeros(m, n);
+U = eye(m, n);
+P = eye(m, n);
+
+for k = 1:m
+    if A(k, k) == 0
+        non_zero_col = find(A(k, :), 1);
+        A(:, [k non_zero_col]) = deal(A(:, [non_zero_col k]));
+        P(:, [k non_zero_col]) = deal(P(:, [non_zero_col k]));
+    end
+end
+
+for k = 1:n
+    
+    for j = k : n
+        U(k, j) = A(k, j) - dot(L(k, 1:k-1), U(1:k-1, j));
+    end
+    
+    for i = k:n
+        L(i, k) = (A(i, k) - dot(L(i,1:k-1), U(1:k-1, k))) / U(k,k);
+    end
+    
+    
+end
+
+end

01_Direct-Methods-for-Solving-Linear-Systems/Report/bibliography.bib

@@ -0,0 +1,13 @@
+@misc{Zdunek,
+    author = {Zdunek, Rafał},
+    title = {Lecture notes in {Numerical Algorithms}},
+    year = {2021}
+}
+
+@Book{GoluVanl96,
+  Title                    = {Matrix Computations},
+  Author                   = {Golub, Gene H. and Van Loan, Charles F.},
+  Publisher                = {The Johns Hopkins University Press},
+  Year                     = {1996},
+  Edition                  = {Third}
+}

01_Direct-Methods-for-Solving-Linear-Systems/Report/main.tex

@@ -0,0 +1,108 @@
+\documentclass[a4paper]{article}
+\usepackage{WUST-AAE-NMaO-Report}
+\usepackage{natbib}
+
+\usepackage{amsmath}
+% \usepackage{amssymb}
+% \usepackage[polish,main=english]{babel}
+% \usepackage{titlesec}
+% \usepackage{vmargin}
+% \usepackage{graphicx}
+% \graphicspath{{images/}}
+% \usepackage{fancyhdr}
+% \usepackage{systeme}
+% \usepackage{xcolor}
+% \usepackage{indentfirst}
+% \usepackage{xspace}
+% \newcommand{\MATLAB}{\textsc{Matlab}\xspace}
+% \usepackage{listings}
+% \usepackage{matlab-prettifier}
+% \usepackage{hyperref}
+% \usepackage{multicol}
+
+
+% \newcommand{\sectionbreak}{\clearpage}
+% \setmarginsrb{3 cm}{2.5 cm}{3 cm}{2.5 cm}{1 cm}{1.5 cm}{1 cm}{1.5 cm}
+
+% Augmented matrix
+% \newenvironment{amatrix}[2]{%
+  % \left[\begin{array}{@{}*{#1}{c} | c *{#2}{c} @{}}
+% }{%
+  % \end{array}\right]
+% }
+
+% Bold letters in matrices
+% \newcommand{\matr}[1]{\mathbf{#1}}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\title{Numerical Algorithms – Report 1: Direct Methods for Solving Linear Systems}
+\author{Sergiusz Warga 230757}
+\date{\today}
+\reporttutor{dr hab. inż. Rafał Zdunek}
+
+\begin{document}
+
+\maketitle
+\tableofcontents
+\pagebreak
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Problems calculations}
+\input{problems/Problem_1}
+% \input{problems/Problem2}\\
+% \input{problems/Problem3}\\
+% \input{problems/Problem4}\\
+% \input{problems/Problem5}\\
+% \input{problems/Problem6}\\
+% \input{problems/Problem7}\\
+% \input{problems/Problem8}\\
+% \input{problems/Problem9}\\
+% \input{problems/Problem10}\\
+% \input{problems/Problem11}\\
+% \input{problems/Problem12}\\
+% \input{problems/Problem13}\\
+% \input{problems/Problem14}\\
+% \input{problems/Problem15}
+
+% \section{Algorithms}
+% \subsection{Algorithm 1 -- Gaussian elimination}\label{algorithm:1}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg1.m}
+% \subsection{Algorithm 2 -- Gaussian  elimination  with  Complete  Pivoting}\label{algorithm:2}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg2.m}
+% \subsection{Algorithm 3 -- Forward Substitution}\label{algorithm:3}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg3.m}
+% \subsection{Algorithm 4 -- Back Substitution}\label{algorithm:4}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg4.m}
+% \subsection{Algorithm 5 -- The Gauss-Jordan elimination algorithm}\label{algorithm:5}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg5.m}
+% \subsection{Algorithm 6 -- The RREF algorithm}\label{algorithm:6}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg6_RREF.m}
+% \subsection{Algorithm 7 -- The LU  factorization without pivoting}\label{algorithm:7}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg7.m}
+% \subsection{Algorithm 8 -- The LU  factorization with partial pivoting}\label{algorithm:8}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg8.m}
+% % \subsection{Algorithm 9 -- A family of the Cholesky  factorization  algorithms}\label{algorithm:9}
+% % \lstinputlisting[style=Matlab-editor]{algorithms/Alg9.m}
+% \stepcounter{subsection}
+% \subsection{Algorithm 10 -- The Cholesky factorization}\label{algorithm:10}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg10.m}
+% \subsection{Algorithm 11 -- The  QR  algorithm  by  the  Householder  transformation}\label{algorithm:11}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg11.m}
+% \subsection{Algorithm 12 -- The QR algorithm by the Givens rotations}\label{algorithm:12}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg12.m}
+% \subsection{Algorithm 13 -- The QR algorithm by the Gram-Schmidt orthogonalization}\label{algorithm:13}
+% \lstinputlisting[style=Matlab-editor]{algorithms/Alg13.m}
+
+%%%%%%%%%%%%%%%%%%%
+%% BIBLIOGRAPHY %%%
+%%%%%%%%%%%%%%%%%%%
+
+\clearpage
+
+\nocite{Zdunek, GoluVanl96}
+\bibliographystyle{alpha}
+\bibliography{bibliography}
+
+\end{document}

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem1.m

@@ -0,0 +1,13 @@
+A = [2, -1, 0, 0;
+    -1, 2, -1, 0;
+    0, -1, 2, -1;
+    0, 0, -1, 2];
+b = [0;0;0;5];
+B = Alg5([A b])
+
+B =
+
+    1.0000         0         0         0    1.0000
+         0    1.0000         0         0    2.0000
+         0         0    1.0000         0    3.0000
+         0         0         0    1.0000    4.0000

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem10.tex

@@ -0,0 +1,58 @@
+\subsection{Problem 10}
+Transform  the  following  matrix  to  the  RREF,  determine $rank(\mathbf{A})$ and  identify  the columns corresponding to the basic and free variables. 
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    2 & 4 & 4 & 6 & 2 \\
+    3 & 6 & 6 & 9 & 6 \\
+    1 & 2 & 4 & 5 & 3 
+    \end{bmatrix}
+\end{equation*}
+
+Check whether the symmetric matrix $\mathbf{A}$ is positive-definite. If so, apply the Cholesky factorization. Then, compute its inverse.
+
+\subsubsection*{Solution}
+
+\begin{equation*}
+    \begin{split}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    2 & 4 & 4 & 6 & 2 \\
+    3 & 6 & 6 & 9 & 6 \\
+    1 & 2 & 4 & 5 & 3 
+    \end{bmatrix}
+    \xrightarrow{\substack{R_2 - 2R_1 \\ R_3 - 3R_1 \\ R_4 - R_1}}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    0 & 0 & 0 & 0 & 0 \\
+    0 & 0 & 0 & 0 & 3 \\
+    0 & 0 & 2 & 2 & 2
+    \end{bmatrix}
+    \xrightarrow{pivot}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    0 & 0 & 2 & 2 & 2 \\
+    0 & 0 & 0 & 0 & 3 \\
+    0 & 0 & 0 & 0 & 0
+    \end{bmatrix}
+    \\
+    \xrightarrow{\substack{R_2 / 2 \\ R_3 / 3}}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 1 \\
+    0 & 0 & 1 & 1 & 1 \\
+    0 & 0 & 0 & 0 & 1 \\
+    0 & 0 & 0 & 0 & 0
+    \end{bmatrix}
+    \xrightarrow{\substack{R_1 - R_3 \\ R_2 - R_3}}
+    \begin{bmatrix}
+    1 & 2 & 2 & 3 & 0 \\
+    0 & 0 & 1 & 1 & 0 \\
+    0 & 0 & 0 & 0 & 1 \\
+    0 & 0 & 0 & 0 & 0
+    \end{bmatrix}
+    \end{split}
+\end{equation*}
+
+Rank of the given matrix is equal to 3.
+

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem11.tex

@@ -0,0 +1,17 @@
+\subsection{Problem 11}
+
+Compute the LU factorization for
+
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+    \phantom{-}1 & \phantom{-}3 & 3 & 2 \\
+    \phantom{-}2 & \phantom{-}6 & 9 & 5 \\
+    -1 & -3 & 3 & 0
+    \end{bmatrix}
+\end{equation*}
+
+
+
+
+Determine a set of basic variables and a set of free variables, and find a homogeneous solution to $\mathbf{Ax = 0}$. What is the rank of $\mathbf{A}$? 

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem12.tex

@@ -0,0 +1,14 @@
+\subsection{Problem 12}
+
+Compute the QR factorization of the matrix:
+
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+    0 & -1 & -3  \\
+    0 & \phantom{-}0 & -2 \\
+    0 & -2 & -1 
+    \end{bmatrix}
+\end{equation*}
+
+How many flops (multiplications/divisions, additions/subtractions) are needed to perform the QR factorization with the Householder transformations and Givens rotations?

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem13.tex

@@ -0,0 +1,3 @@
+\subsection{Problem 13}
+
+Let $\mathbf{Ax = b}$ , where $A = I_N \bigotimes C^TC$, the symbol $\bigotimes$ denotes the Kronecker product, $I_N \in R^{N \times N}$ is an identity matrix, $C \in R^{M \times M}$ is a random matrix with a uniform distribution, $M = 100$, and $N=50$, and $x \square N(0, I_{MN})$. Find the direct method that solves the above system  of  linear  equations  with  the  lowest  computational  cost.  Estimate  the  cost  with  a  roughly calculated number of flops and with the elapsed time.

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem14.tex

@@ -0,0 +1,4 @@
+\subsection{Problem 14}
+
+ Solve $\mathbf{Ax=b}$, where $A \in R ^ {10 \times 10}$ is  the  Hilbert  matrix, and,  $x \square N(0, I_{10})$ with the Gaussian elimination in Matlab. Then make a small change in an entry of $\mathbf{A}$ or $\mathbf{b}$, and compare the solutions.
+  

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem15.tex

@@ -0,0 +1,59 @@
+\subsection{Problem 15}
+
+Find the currents in the circuit, assuming all resistances are 10 Ohm.
+\begin{figure}[h]
+    \centering
+    \includegraphics[width=0.5\textwidth]{images/electro.png}
+    \label{fig:electro}
+\end{figure}
+
+\subsubsection*{Solution}
+
+\begin{equation*}
+    \left\{\begin{matrix}
+I_{L1}(R_1 + R_2) - I_{L3}(R_2) = -E_2\\ 
+I_{L2}(R_3 + R_4) - I_{L3}(R_3) = E_2\\ 
+-I_{L1}(R_2) - I_{L2}(R_3) + I_{L3}(R_2 + R_3) = -E_1
+\end{matrix}\right.
+\end{equation*}
+
+After substitution we have:
+
+\begin{equation*}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{3}{1}
+        20 & 0 & -10 & -10 \\
+        0 & 20 & -10 & 10 \\
+        -10 & -10 & 20 & -20
+    \end{amatrix}
+    \xrightarrow[magic]{}
+    \begin{amatrix}{3}{1}
+        1 & 0 & 0 & -1.5 \\
+        0 & 1 & 0 & -0.5 \\
+        0 & 0 & 1 & -2
+    \end{amatrix}
+\end{equation*}
+
+\begin{equation*}
+\left\{\begin{matrix}
+I_1 = -I_{L3} \\
+I_2 = -I_{L1} + I_{L2} \\
+I_3 = I_{L1} - I_{L3} \\
+I_4 = -I_{L1} \\
+I_5 = I_{L2} - I_{L3} \\
+I_6 = -I_{L2}
+\end{matrix}\right.
+\end{equation*}
+
+\begin{equation*}
+\left\{\begin{matrix}
+I_1 = 2 \\
+I_2 = 1 \\
+I_3 = 0.5 \\
+I_4 = 1.5 \\
+I_5 = 1.5 \\
+I_6 = 0.5
+\end{matrix}\right.
+\end{equation*}

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem2.m

@@ -0,0 +1,15 @@
+A = [1, 1, 1;
+     1, 1, 2;
+     1, 2, 2];
+b = [1;2;1];
+[P, Q, L, U] = Alg2(A);
+% PAQ = LU
+% Ly = b and Ux = y
+y = Alg3(L, P*b);     % Forward substitution
+x = Q*Alg4(U, y)      % Backward substitution
+
+x =
+
+     1
+    -1
+     1

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem2.tex

@@ -0,0 +1,52 @@
+\subsection{Problem 2}\label{problem:2}
+Solve the following system of linear equations using the Gaussian elimination with pivoting:
+\begin{equation*}
+    \systeme{x_1+x_2+x_3=1,x_1+x_2+2x_3=2,x_1+2x_2+2x_3=1}
+\end{equation*}
+Explain why the Gaussian elimination without pivoting does not work.
+\subsubsection*{Mathematics}
+As was stated at the end of the theoretical introduction to the previous problem, Gaussian Elimination algorithms require the diagonal entry to be different than zero (to not normalize the row dividing by zero). While in manual computations such a division would be noticed, in \MATLAB their result is \texttt{Inf} and algorithm proceeds with replacing subsequent coefficients with this value, resulting in gibberish output. With pivoting this behaviour can be avoided. The basic idea is to find a maximum absolute value in the zeroed column below the pivot and interchange their rows. This is called a partial pivoting. Searching in the whole sub-matrix spanned between pivot and $A_{n,m}$, therefore interchanging not only rows but also columns, is called a complete pivoting.
+\subsubsection*{Solution}
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{3}{1}
+        1 & 1 & 1 & 1\\
+        1 & 1 & 2 & 2\\
+        1 & 2 & 2 & 1
+    \end{amatrix} \xrightarrow[\substack{R_2-R_1\\R_3-R_1}]{}
+    \begin{amatrix}{3}{1}
+        1 & 1 & 1 & 1\\
+        0 & 0 & 1 & 1\\
+        0 & 1 & 1 & 0
+    \end{amatrix} \xrightarrow[R_2\leftrightarrow R_3]{}
+    \begin{amatrix}{3}{1}
+        1 & 1 & 1 & 1\\
+        0 & 1 & 1 & 0\\
+        0 & 0 & 1 & 1
+    \end{amatrix}\rightarrow\\
+    \rightarrow
+    \systeme{x_1+x_2+x_3=1,x_2+x_3=0,x_3=1}
+    \rightarrow
+    \systeme{x_1+x_2+x_3=\phantom{-}1,x_2=-1,x_3=\phantom{-}1}
+    \rightarrow
+    \begin{cases}
+    x_1=\phantom{-}1\\
+    x_2=-1\\
+    x_3=\phantom{-}1
+    \end{cases}
+\end{split}
+\end{equation*}
+
+To verify manual solution the Algorithm~\ref{algorithm:2} was used. As is explained in detail in~\cite[section 3.4.2]{GoluVanl96} it transforms matrix $\matr{P}\matr{A}\matr{U}$ into a form from which lower unit triangular and upper triangular matrices may be extracted, so that:
+\begin{align*}
+   \matr{P}\matr{A}\matr{U} &= \matr{L}\matr{U} & \matr{A}\matr{x}&=\matr{b}
+\end{align*}
+therefore
+\begin{align*}
+    \matr{A}&=\matr{P^T}\matr{L}\matr{U}\matr{Q^T} & \matr{P^T}\matr{L}\matr{U}\matr{Q^T}\matr{x} &= \matr{b}\\
+    & & \matr{L}\matr{U}\matr{Q^T}\matr{x} &= \matr{P}\matr{b}
+\end{align*}
+\lstinputlisting[style=Matlab-editor]{problems/Problem2.m}

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem3.tex

@@ -0,0 +1,140 @@
+\subsection{Problem 3}
+Solve the system:
+\begin{equation*}
+    \systeme{0.0001x_1+x_2=1,x_1+x_2=2}
+\end{equation*}
+with the Gaussian elimination with and without pivoting at the round-off error limited to 3 significant digits. Compute the condition number of the system matrix.
+\subsubsection*{Mathematics}
+Not only pivoting helps with $a_{ii}=0$ problem but, by choosing the biggest entry, ensures the stability of the LU factorization. As explained in detail in~~\cite[section 3.3]{GoluVanl96} for the factorization errors to be small, the diagonal entries of $\matr{U}$ matrix has to be small. Please note that in both solutions and condition number computations precision was limited to 3 significant digits.
+\subsubsection*{Solution without pivoting}
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{2}{1}
+        10^{-4} & 1 & 1\\
+        1 & 1 & 2
+    \end{amatrix}
+    \xrightarrow[\substack{R_2-10^{4}R_1}]{}
+    \begin{amatrix}{2}{1}
+        10^{-4} & 1 & 1\\
+        0 & -10^{4} & -10^{4}
+    \end{amatrix}
+\end{split}
+\end{equation*}
+gives us
+\begin{equation*}
+    \begin{cases}
+    x_1=0\\
+    x_2=1\\
+    \end{cases}
+\end{equation*}
+which is obviously contradictory to $R_2$.
+
+\subsubsection*{Solution with pivoting}
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{2}{1}
+        10^{-4} & 1 & 1\\
+        1 & 1 & 2
+    \end{amatrix}
+    \xrightarrow[R_1\leftrightarrow R_2]{}
+    \begin{amatrix}{2}{1}
+        1 & 1 & 2 \\
+        10^{-4} & 1 & 1
+    \end{amatrix}
+    \xrightarrow[\substack{R_2-10^{-4}R_1}]{}
+    \begin{amatrix}{2}{1}
+        1 & 1 & 2 \\
+        0 & 1 & 1
+    \end{amatrix}
+\end{split}
+\end{equation*}
+gives us
+\begin{equation*}
+    \begin{cases}
+    x_1=1\\
+    x_2=1
+    \end{cases}
+\end{equation*}
+\subsubsection*{Condition number}
+Condition number should describe the stability of a computed matrix. The bigger the condition number $\kappa$ the less stable the (coefficient) matrix is. It is computed as: 
+\begin{equation}
+    \kappa(\mathbf{A})=||\mathbf{A}||_2\cdot||\mathbf{A}^{-1}||_2
+\end{equation}
+where 
+\begin{equation}\label{eqn:2-norm}
+\begin{split}
+    ||\mathbf{A}||_2 = \sqrt{\lambda_{max}A^TA}
+\end{split}
+\end{equation}
+$\matr{A}^{-1}$ may be computed via Gauss-Jordan algorithm with pivoting:
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{I}
+    \end{amatrix} = 
+    \begin{amatrix}{2}{2}
+        \frac{1}{1000} & 1 & 1 & 0\\
+        1 & 1 & 0 & 1
+    \end{amatrix}
+    \xrightarrow[R_1\leftrightarrow R_2]{}
+    \begin{amatrix}{2}{2}
+        1 & 1 & 0 & 1\\
+        \frac{1}{1000} & 1 & 1 & 0
+    \end{amatrix}
+    \xrightarrow[\substack{R_2-\frac{1}{1000}R_2}]{}\\
+    \xrightarrow[\substack{R_2-\frac{1}{1000}R_2}]{}
+    \begin{amatrix}{2}{2}
+        1 & 1 & 0 & 1\\
+        0 & 1 & 1 & 0
+    \end{amatrix}
+    \xrightarrow[\substack{R_1-R_2}]{}
+    \begin{amatrix}{2}{2}
+        1 & 0 & -1 & 1\\
+        0 & 1 & 1 & 0
+    \end{amatrix} = 
+    \begin{amatrix}{1}{1}
+        \matr{I} & \matr{A^{-1}}
+    \end{amatrix}
+\end{split}
+\end{equation*}
+Now, using the equation~\ref{eqn:2-norm}
+
+\parbox{0.5\textwidth}{
+\begin{gather*}
+    \matr{A}^T\matr{A} = 
+    \begin{bmatrix}
+        1 & 1 \\
+        1 & 2
+    \end{bmatrix} \\
+    \det\begin{pmatrix}
+        1-\lambda & 1 \\
+        1 & 2-\lambda
+    \end{pmatrix} = 0 \\
+    \lambda=\frac{3\pm\sqrt{5}}{2}
+    \end{gather*}
+}
+\parbox{0.5\textwidth}{
+\begin{gather*}
+     \matr{A}^{-1^T}\matr{A}^{-1} = 
+    \begin{bmatrix}
+        2 & -1 \\
+        -1 & 1
+    \end{bmatrix} \\
+    \det\begin{pmatrix}
+        2-\lambda & -1 \\
+        -1 & 1-\lambda
+    \end{pmatrix} = 0 \\
+    \lambda=\frac{3\pm\sqrt{5}}{2}
+\end{gather*}
+}
+we finally obtain
+\begin{equation*}
+    \kappa(\mathbf{A})=\frac{3\pm\sqrt{5}}{2}\approx2.618
+\end{equation*}
+which is a very low number indicating, that the matrix $\matr{A}$ is computationally stable.

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem4.m

@@ -0,0 +1,23 @@
+A = [0.835, 0.667;
+     0.333, 0.266];
+b = [0.168; 0.067];
+bp = [0.168; 0.066];
+B = Alg5([A b])
+Bp = Alg5([A bp])
+kappa = cond(A)
+
+B =
+
+    1.0000         0    1.0000
+         0    1.0000   -1.0000
+
+
+Bp =
+
+    1.0000         0 -666.0000
+         0    1.0000  834.0000
+
+
+kappa =
+
+   1.3238e+06

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem4.tex

@@ -0,0 +1,37 @@
+\subsection{Problem 4}
+Solve the system of linear equations:
+\begin{equation*}
+    \systeme{0.835x_1 + 0.667x_2 = 0.168,0.333x_1 + 0.266x_2 = 0.067}
+\end{equation*}
+
+Then slightly perturb $b_2$ from $0.067$ to $0.066$ and compute the solution to the perturbed system. Explain the change in the solution by computing the condition number of the system matrix.
+\subsubsection*{Solution}
+\begin{equation*}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{2}{1}
+        0.835 & 0.667 & 0.168\\
+        0.333 & 0.266 & 0.067
+    \end{amatrix} \xrightarrow[R_2-0.333R_1]{\frac{R_1}{0.835}}
+    \begin{amatrix}{2}{1}
+        1 & 0.7988 & 0.2012\\
+        0 & 0 & 0
+    \end{amatrix}
+\end{equation*}
+
+which yields infinite many solutions. Perturbing $b_2$ we obtain:
+\begin{equation*}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b'}
+    \end{amatrix} \xrightarrow[R_2-0.333R_1]{\frac{R_1}{0.835}}
+    \begin{amatrix}{2}{1}
+        1 & 0.7988 & 0.2012\\
+        0 & 0 & -0.001
+    \end{amatrix}
+\end{equation*}
+
+which yields none solutions. Performing the same calculations with Algorithm~\ref{algorithm:5}:
+\lstinputlisting[style=Matlab-editor]{problems/Problem4.m}
+
+With enough precision \MATLAB presented two valid solutions. Both of them differed by a factor of hundreds, which may be explained by a huge coefficient number (way over a million).

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem5.m

@@ -0,0 +1,17 @@
+A = [2, 1, 2;
+     1, 2, 3;
+     4, 1, 2];
+[P, Q, L, U] = Alg2(A);
+I = P*eye(3);
+% Ly = b and Ux = y
+y = Alg3(L, I);     % Forward substitution
+x = Q*Alg4(U, y);   % Backward substitution
+inv(A) - x
+
+ans =
+
+   1.0e-15 *
+
+   -0.2220    0.0555    0.1110
+    0.8882   -0.4441   -0.6661
+         0    0.1110         0

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem5.tex

@@ -0,0 +1,17 @@
+\subsection{Problem 5}
+Find $\mathbf{A^{-1}}$ to:
+
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+        2 & 1 & 2 \\
+        1 & 2 & 3 \\
+        4 & 1 & 2 
+    \end{bmatrix}
+\end{equation*}
+
+solving the system $\mathbf{AX=I_{3}}$.
+
+\subsubsection*{Solution}
+\lstinputlisting[style=Matlab-editor]{problems/Problem5.m}
+which judges well our algorithms.

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem6.m

@@ -0,0 +1,26 @@
+A = [1 2 3 4;
+    -1 1 2 1;
+     0 2 1 3;
+     0 0 1 1];
+b = [1; 1; 1; 1];
+% PA = LU
+[L, U, P] = Alg8(A);
+det(A) - prod(diag(P.'*U))
+% Ly = b and Ux = y
+y = Alg3(L, P*b);   % Forward substitution
+x = Alg4(U, y);   % Backward substitution
+x - A\b
+
+ans =
+
+     0
+
+
+ans =
+
+   1.0e-15 *
+
+    0.4441
+         0
+         0
+         0

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem6.tex

@@ -0,0 +1,26 @@
+\subsection{Problem 6}
+
+Apply the LU factorization to the matrix
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+        \phantom{-}1 & 2 & 3 & 4 \\
+        -1 & 1 & 2 & 1 \\
+        \phantom{-}0 & 2 & 1 & 3 \\
+        \phantom{-}0 & 0 & 1 & 1
+    \end{bmatrix}
+\end{equation*}
+Then calculate $\det(\matr{A})$ using the matrix $\matr{U}$. Finally solve $\matr{A}\matr{x}=\matr{b}$ for $\matr{b}=[1\dots1]^T$.
+\subsubsection*{Mathematics}
+The LU Factorization decomposes matrix $\matr{A}$ into an upper triangular matrix $\matr{U}$ and unit lower triangular matrix $\matr{L}$, so that
+\begin{equation*}
+    \matr{A} = \matr{L}\matr{U}
+\end{equation*}
+Using the property of a determinant and a fact, that $\det(\matr{L}) = 1$ one can calculate $\det(\matr{A})$ as
+\begin{equation*}
+    \det(\matr{A}) = \det(\matr{L}\matr{U}) = \det(\matr{L}) \cdot \det(\matr{U}) = \det(\matr{U}) = u_{11}\cdots u_{nn}
+\end{equation*}
+\subsubsection*{Solution}
+A glance at the matrix $\matr{A}$ tells us that it is not strictly diagonally dominant, so a pivoting algorithm should be used.
+%\footnote{$\matr{A}\in\mathbb{R}^{n\times n}$ is \textit{strictly diagonally dominant} if $|a_{ii}|>\sum_{j=1 j i}^n|a_{ij}|$}
+\lstinputlisting[style=Matlab-editor]{problems/Problem6.m}

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem7.m

@@ -0,0 +1,24 @@
+A = hilb(5);
+[L, U] = Alg7(A);
+cond(A)
+L*U - A
+prod(diag(U)) - det(A)
+
+ans =
+
+   4.7661e+05
+   
+ans =
+
+   1.0e-16 *
+
+         0         0         0         0         0
+         0         0         0         0         0
+         0         0         0         0         0
+         0         0         0         0   -0.1388
+         0         0         0   -0.1388         0
+
+
+ans =
+
+   1.6620e-23

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem7.tex

@@ -0,0 +1,21 @@
+\subsection{Problem 7}
+
+Let $\mathbf{A} = \left[a_{ij}\right] \in\mathbb{R}^{N\times N}$ with, $a_{ij} = \frac{1}{i + j - 1}$ (Hilbert  matrix). For  $N=5$, perform  the LU factorization of the matrix $\mathbf{A}$. Then, compute det($\mathbf{A}$).
+\subsubsection*{Mathematics}
+Hilbert matrices are an example of ill-conditioned matrices, which -- with their high $\kappa$ -- makes the numerical computations highly unstable. For our  $\mathbf{H} = \left[h_{ij}\right] \in\mathbb{R}^{5\times 5}$:
+
+\begin{equation*}
+    \matr{H} = 
+    \begin{bmatrix}
+               1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \\
+     \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\
+     \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} \\
+     \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} \\
+     \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} 
+    \end{bmatrix}
+\end{equation*}
+
+$\kappa(\matr{H})\approx4.7\cdot10^5$ which is a huge value similar to this in the Problem 4.
+\subsubsection*{Solution}
+\lstinputlisting[style=Matlab-editor]{problems/Problem7.m}
+The above results may look promising, but a determinant calculated both by \MATLAB and our algorithms is of the order of $4\cdot10^{-12}$, which for numerical computations is highly not useful.

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem8.m

@@ -0,0 +1,32 @@
+A = [ 1 -1  0  0;
+     -1  2 -1  0;
+      0 -1  2 -1;
+      0  0 -1  2];
+if all(eig(A) > 0)
+    C = Alg10(A)
+    C*C.' - A
+    inv(C)
+end
+
+C =
+
+     1     0     0     0
+    -1     1     0     0
+     0    -1     1     0
+     0     0    -1     1
+
+
+ans =
+
+     0     0     0     0
+     0     0     0     0
+     0     0     0     0
+     0     0     0     0
+
+
+ans =
+
+     1     0     0     0
+     1     1     0     0
+     1     1     1     0
+     1     1     1     1

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem8.tex

@@ -0,0 +1,16 @@
+\subsection{Problem 8}
+Let
+\begin{equation*}
+    \matr{A} = 
+    \begin{bmatrix}
+        \phantom{-}1 & -1 & \phantom{-}0 & \phantom{-}0 \\
+        -1 & \phantom{-}2 & -1 & \phantom{-}0 \\
+        \phantom{-}0 & -1 & \phantom{-}2 & -1 \\
+        \phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}2
+    \end{bmatrix}
+\end{equation*}
+
+Check whether the symmetric matrix $\mathbf{A}$ is positive-definite. If so, apply the Cholesky factorization. Then, compute its inverse.
+
+\lstinputlisting[style=Matlab-editor]{problems/Problem8.m}
+

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem9.tex

@@ -0,0 +1,3 @@
+\subsection{Problem 9}
+Let  $\mathbf{A}$  be  the  Pascal  matrix  of  order  100  (\textit{pascal}  function  in Matlab).  Check  whether the matrix is positive definite. If so, apply the Cholesky factorization and give interpretation of the obtained factor.
+

01_Direct-Methods-for-Solving-Linear-Systems/Report/problems/Problem_1.tex

@@ -0,0 +1,65 @@
+\subsection{Problem 1}
+
+Solve the system and find the pivots when:
+\begin{equation*}
+    \systeme{2u-v=0,-u+2v-w=0,-v+2w-z=0,-w+2z=5}
+\end{equation*}
+
+\subsubsection*{Mathematics}
+The solution to a linear system like the one presented above may be obtained via Gauss-Jordan Elimination. An augmented matrix
+$\matr{B} = \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+\end{amatrix}$, where $\matr{A}$ is a coefficient matrix, and $\matr{b}$ is a column vector of constant terms, is assembled. Then, the algorithm implemented in \ref{algorithm:5} is performed. Note, that none of the $a_{ii}$ can be equal to zero for this algorithm to work properly\footnote{Explanation of this property is presented in \ref{problem:2}}.
+
+\subsubsection*{Solution}
+
+Manual calculations were performed. In the following transformations a row normalization is presented above the arrow and then a zeroing is below it.
+
+\begin{equation*}
+\begin{split}
+    \begin{amatrix}{1}{1}
+        \matr{A} & \matr{b}
+    \end{amatrix} = 
+    \begin{amatrix}{4}{1}
+        \phantom{-}2 & -1 & \phantom{-}0 & \phantom{-}0 & 0\\
+        -1 & \phantom{-}2 & -1 & \phantom{-}0 & 0\\
+        \phantom{-}0 & -1 & \phantom{-}2 & -1 & 0\\
+        \phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}2 & 5
+    \end{amatrix} \xrightarrow[R_2 + R_1]{\frac{1}{2}R_1} 
+    \begin{amatrix}{4}{1}
+        \phantom{-}1 & -\frac{1}{2} & \phantom{-}0 & \phantom{-}0 & 0\\
+        \phantom{-}0 & \phantom{-}\frac{3}{2} & -1 & \phantom{-}0 & 0\\
+        \phantom{-}0 & -1 & \phantom{-}2 & -1 & 0\\
+        \phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}2 & 5
+    \end{amatrix} \xrightarrow[\substack{R_1+\frac{1}{2}R_2\\R_3+R_2}]{\frac{2}{3}R_2} \\
+    \xrightarrow[R_3+\frac{2}{3}R_2]{\frac{2}{3}R_2}
+    \begin{amatrix}{4}{1}
+        \phantom{-}1 & \phantom{-}0 & -\frac{1}{3} & \phantom{-}0 & 0\\
+        \phantom{-}0 & \phantom{-}1 & -\frac{2}{3} & \phantom{-}0 & 0\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}\frac{4}{3} & -1 & 0\\
+        \phantom{-}0 & \phantom{-}0 & -1 & \phantom{-}2 & 5
+    \end{amatrix} \xrightarrow[\substack{R_1+\frac{1}{3}R_3\\R_2+\frac{2}{3}R_3\\R_4+R_3}]{\frac{3}{4}R_3}
+    \begin{amatrix}{4}{1}
+        \phantom{-}1 & \phantom{-}0 & \phantom{-}0 & -\frac{1}{4} & 0\\
+        \phantom{-}0 & \phantom{-}1 & \phantom{-}0 & -\frac{1}{2} & 0\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}1 & -\frac{3}{4} & 0\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}0 & -\frac{4}{3} & 5
+    \end{amatrix} \xrightarrow[\substack{R_1+\frac{1}{4}R_4\\R_2+\frac{1}{2}R_4\\R_3+\frac{3}{4}R_4}]{-\frac{3}{4}R_4}\\
+    \xrightarrow[\substack{R_1+\frac{1}{4}R_4\\R_2+\frac{1}{2}R_4\\R_3+\frac{3}{4}R_4}]{-\frac{3}{4}R_4}
+    \begin{amatrix}{4}{1}
+        \phantom{-}1 & \phantom{-}0 & \phantom{-}0 & \phantom{-}0 & 1\\
+        \phantom{-}0 & \phantom{-}1 & \phantom{-}0 & \phantom{-}0 & 2\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}1 & \phantom{-}0 & 3\\
+        \phantom{-}0 & \phantom{-}0 & \phantom{-}0 & \phantom{-}1 & 4
+    \end{amatrix} \rightarrow
+    \begin{cases}
+    u=1\\
+    v=2\\
+    w=3\\
+    z=4
+    \end{cases}
+\end{split}
+\end{equation*}
+
+Above result was verified with the Algorithm~\ref{algorithm:5}:
+\lstinputlisting[style=Matlab-editor]{problems/Problem1.m}

02_Iterative-Linear-Solvers/Report/main.tex

@@ -0,0 +1,15 @@
+\documentclass[a4paper]{article}
+\usepackage{WUST-AAE-NMaO-Report}
+
+\title{Numerical Algorithms – Report 1: Iterative Linear Solvers}
+\author{Sergiusz Warga 230757}
+\date{\today}
+\reporttutor{dr hab. inż. Rafał Zdunek}
+
+\begin{document}
+
+\maketitle
+\tableofcontents
+\pagebreak
+
+\end{document}

Direct Methods for Solving Linear Systems/README.md

@@ -1,8 +0,0 @@
-# Direct Methods for Solving Linear Systems
-
-## Bibliography
-
-References in the code comments are represented by `[number]`:
-
-  - 1
-  - 2 – Golub, van Loan, Matrix Computations, 3rd edition

Dockerfile

@@ -0,0 +1,6 @@
+FROM alpine:3.16.0
+
+RUN apk update && \
+    apk add --no-cache \
+    texlive-full=20220403.62885-r2 \
+    make=4.3-r0

LICENSE

@@ -0,0 +1,13 @@
+            DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
+                    Version 2, December 2004
+
+ Copyright (C) 2004 Sam Hocevar <sam@hocevar.net>
+
+ Everyone is permitted to copy and distribute verbatim or modified
+ copies of this license document, and changing it is allowed as long
+ as the name is changed.
+
+            DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
+   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
+
+  0. You just DO WHAT THE FUCK YOU WANT TO.

Makefile

@@ -0,0 +1,63 @@
+BUILDDIR = `pwd`/build
+DOCKER_IMAGE = latex-builder
+LATEX_OPTIONS = -pdf
+#  -interaction=nonstopmode
+PACKAGE_NAME = WUST-AAE-NMaO-Report
+
+SOURCE_DIRS = $(shell ls -d */ | grep -v "^build")
+SOURCE_TEXS = $(SOURCE_DIRS:=Report/main.tex)
+PDFBUILDS = $(addprefix $(BUILDDIR)/, $(SOURCE_DIRS:/=.pdf))
+
+
+.DEFAULT_GOAL := docker
+
+
+.PHONY : clean
+clean :
+	-rm -r $(BUILDDIR) $(PACKAGE_NAME).sty
+
+
+.PHONY : all
+all : $(PDFBUILDS)
+
+
+.PHONY : sty
+sty : $(PACKAGE_NAME).sty
+$(PACKAGE_NAME).sty : $(PACKAGE_NAME).ins $(PACKAGE_NAME).dtx
+	mkdir -p $(BUILDDIR)
+	yes | latex -output-directory=$(BUILDDIR) $(PACKAGE_NAME).ins
+	cp $(BUILDDIR)/$@ $@
+
+
+$(PDFBUILDS) : $(BUILDDIR)/%.pdf: $(PACKAGE_NAME).sty %/Report/main.tex
+	mkdir -p $(BUILDDIR)
+	cp logo-pwr-2016.pdf $(BUILDDIR)/
+	latexmk $(LATEX_OPTIONS) \
+		-cd \
+		-jobname=$* \
+		-output-directory=$(BUILDDIR) \
+		$*/Report/main
+
+
+.PHONY : hadolint
+hadolint : Dockerfile
+	docker run --rm --interactive hadolint/hadolint < Dockerfile
+
+
+.PHONY : image
+image : Dockerfile
+	docker build \
+		--tag $(DOCKER_IMAGE) \
+		.
+
+
+.PHONY : docker
+docker : image ## Compile the project via the latex-builder docker image
+	docker run \
+		--rm \
+		--interactive \
+		--workdir /data \
+		--volume `pwd`:/data \
+		--name=$(DOCKER_IMAGE) \
+		$(DOCKER_IMAGE) \
+		sh -c "make all"

README.md

@@ -0,0 +1 @@
+# Numerical Methods and Optimization

WUST-AAE-NMaO-Report.dtx

@@ -0,0 +1,145 @@
+% \iffalse meta-comment
+%
+% DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
+%
+% Author: Sergiusz Warga <sergiusz.warga@gmail.com>
+%
+% \fi
+%
+% \iffalse
+%<*driver>
+\ProvidesFile{WUST-AAE-NMaO-Report.dtx}
+%</driver>
+%<package>\NeedsTeXFormat{LaTeX2e}[2005/12/01]
+%<package>\ProvidesPackage{WUST-AAE-NMaO-Report}
+%<*package>
+  [2023/01/06 v0.1.0 WUST AAE NMaO Report package]
+%</package>
+%
+%<*driver>
+\documentclass{ltxdoc}
+\usepackage[
+  hidelinks,
+  bookmarks,
+  bookmarksopen,
+  bookmarksopenlevel=2,
+  pdftitle={WUST AAE NMaO Report LaTeX package},
+  pdfauthor={Sergiusz Warga}
+]{hyperref}
+\EnableCrossrefs
+\CodelineIndex
+\RecordChanges
+\setlength\hfuzz{15pt}      % don’t show so many
+\hbadness=7000              % over- and underfull box warnings
+\begin{document}
+  \DocInput{WUST-AAE-NMaO-Report.dtx}
+\end{document}
+%</driver>
+% \fi
+%
+% \CheckSum{0}
+%
+% \changes{v0.1.0}{2023/01/06}{Initial version}
+%
+% \GetFileInfo{WUST-AAE-NMaO-Report.dtx}
+%
+% \DoNotIndex{\@date, \@title}
+% \DoNotIndex{\\}
+% \DoNotIndex{\begin}
+% \DoNotIndex{\DeclareRobustCommand, \dimexpr, \doublerulesep}
+% \DoNotIndex{\end}
+% \DoNotIndex{\fbox, \fboxrule, \fboxsep}
+% \DoNotIndex{\gdef}
+% \DoNotIndex{\hoffset, \Huge}
+% \DoNotIndex{\Large}
+% \DoNotIndex{\newlength, \newpage, \null}
+% \DoNotIndex{\oddsidemargin}
+% \DoNotIndex{\paperwidth}
+% \DoNotIndex{\renewcommand}
+% \DoNotIndex{\setlength}
+% \DoNotIndex{\textwidth}
+% \DoNotIndex{\uppercase}
+% \DoNotIndex{\vfill}
+%
+% \title{WUST AAE NMaO Report package\thanks
+%   {This document corresponds to \texttt{WUST-AAE-NMaO-Report}~\fileversion,
+%       dated \filedate.}}\author{Sergiusz Warga\texttt{sergiusz.warga@gmail.com}}
+%
+% \maketitle
+%
+% \StopEventually{\PrintChanges \PrintIndex}
+%
+% \section{Implementation}
+%
+% \subsection{Packages}
+%    \begin{macrocode}
+\RequirePackage{graphicx}
+\RequirePackage{xspace}
+%    \end{macrocode}
+%
+% \subsection{Dimentions}
+%
+% Text area.
+%
+%    \begin{macrocode}
+\setlength{\textwidth}{140mm}
+%    \end{macrocode}
+%
+% Center the text area.
+%
+%    \begin{macrocode}
+\setlength{\hoffset}{
+  \dimexpr
+  (\paperwidth-\textwidth)/2-\oddsidemargin-1in
+}
+
+\DeclareRobustCommand*{\reportgroup}[1]{\gdef\@reportgroup{#1}}
+\DeclareRobustCommand*{\reporttutor}[1]{\gdef\@reporttutor{#1}}
+\DeclareRobustCommand*{\reportclassdate}[1]{\gdef\@reportclassdate{#1}}
+\DeclareRobustCommand*{\reportlabassistant}[1]{\gdef\@reportlabassistant{#1}}
+
+\reportgroup{}
+\reporttutor{}
+\reportclassdate{}
+\reportlabassistant{}
+
+%    \end{macrocode}
+% \DescribeMacro{\maketitle}
+%    \begin{macrocode}
+\renewcommand{\maketitle}{
+\begin{titlepage}
+  \thispagestyle{empty}
+  \begin{center}
+    \includegraphics{logo-pwr-2016.pdf}
+    \vspace{\baselineskip}
+    \rule{\linewidth}{0.2 mm}\\[0.4 cm]
+    {\huge \bfseries \@title}\\
+    \rule{\linewidth}{0.2 mm}\\[1.5 cm]
+  \vfill
+  \large
+  \textsc{\@author}\\
+  \vspace{\baselineskip}
+  \textsc{Tutor: \@reporttutor}\\
+  \textsc{Laboratory group: Monday 13:15}\\
+
+  \end{center}
+\end{titlepage}
+}
+%    \end{macrocode}
+% MATLAB command.
+%    \begin{macrocode}
+\newcommand{\MATLAB}{\textsc{Matlab}\xspace}
+%    \end{macrocode}
+% Bold letters in matrices
+%    \begin{macrocode}
+\newcommand{\matr}[1]{\mathbf{#1}}
+%    \end{macrocode}
+% Augmented matrix
+%    \begin{macrocode}
+\newenvironment{amatrix}[2]{%
+  \left[\begin{array}{@{}*{#1}{c} | c *{#2}{c} @{}}
+}{%
+  \end{array}\right]
+}
+%    \end{macrocode}
+% \Finale

WUST-AAE-NMaO-Report.ins

@@ -0,0 +1,38 @@
+%%
+%% DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
+%%
+%% Author: Sergiusz Warga <sergiusz.warga@gmail.com>
+%%
+\input docstrip.tex
+\keepsilent % Turn off DocStrip activity reports
+
+\usedir{tex/latex/WUST-AAE-NMaO-Report}
+
+\preamble
+
+This is a generated file.
+
+To what the fuck you want with it.
+
+\endpreamble
+
+\generate{
+\file{WUST-AAE-NMaO-Report.sty}
+{\from{WUST-AAE-NMaO-Report.dtx}{package}}
+}
+
+\obeyspaces
+\Msg{****************************************************}
+\Msg{*                                                  *}
+\Msg{* To finish the installation you have to move the  *}
+\Msg{* following file into a directory searched by TeX: *}
+\Msg{*                                                  *}
+\Msg{*     WUST-AAE-NMaO-Report.sty                     *}
+\Msg{*                                                  *}
+\Msg{* To produce the documentation run the file        *}
+\Msg{* WUST-AAE-NMaO-Report.dtx through LaTeX.          *}
+\Msg{*                                                  *}
+\Msg{* Happy TeXing!                                    *}
+\Msg{*                                                  *}
+\Msg{****************************************************}
+\endbatchfile

WUST-AAE-NMaO-Report.sty

@@ -0,0 +1,62 @@
+%%
+%% This is file `WUST-AAE-NMaO-Report.sty',
+%% generated with the docstrip utility.
+%%
+%% The original source files were:
+%%
+%% WUST-AAE-NMaO-Report.dtx  (with options: `package')
+%% 
+%% This is a generated file.
+%% 
+%% To what the fuck you want with it.
+%% 
+\NeedsTeXFormat{LaTeX2e}[2005/12/01]
+\ProvidesPackage{WUST-AAE-NMaO-Report}
+  [2023/01/06 v0.1.0 WUST AAE NMaO Report package]
+\RequirePackage{graphicx}
+\RequirePackage{xspace}
+\setlength{\textwidth}{140mm}
+\setlength{\hoffset}{
+  \dimexpr
+  (\paperwidth-\textwidth)/2-\oddsidemargin-1in
+}
+
+\DeclareRobustCommand*{\reportgroup}[1]{\gdef\@reportgroup{#1}}
+\DeclareRobustCommand*{\reporttutor}[1]{\gdef\@reporttutor{#1}}
+\DeclareRobustCommand*{\reportclassdate}[1]{\gdef\@reportclassdate{#1}}
+\DeclareRobustCommand*{\reportlabassistant}[1]{\gdef\@reportlabassistant{#1}}
+
+\reportgroup{}
+\reporttutor{}
+\reportclassdate{}
+\reportlabassistant{}
+
+\renewcommand{\maketitle}{
+\begin{titlepage}
+  \thispagestyle{empty}
+  \begin{center}
+    \includegraphics{logo-pwr-2016.pdf}
+    \vspace{\baselineskip}
+    \rule{\linewidth}{0.2 mm}\\[0.4 cm]
+    {\huge \bfseries \@title}\\
+    \rule{\linewidth}{0.2 mm}\\[1.5 cm]
+  \vfill
+  \large
+  \textsc{\@author}\\
+  \vspace{\baselineskip}
+  \textsc{Tutor: \@reporttutor}\\
+  \textsc{Laboratory group: Monday 13:15}\\
+
+  \end{center}
+\end{titlepage}
+}
+\newcommand{\MATLAB}{\textsc{Matlab}\xspace}
+\newcommand{\matr}[1]{\mathbf{#1}}
+\newenvironment{amatrix}[2]{%
+  \left[\begin{array}{@{}*{#1}{c} | c *{#2}{c} @{}}
+}{%
+  \end{array}\right]
+}
+\endinput
+%%
+%% End of file `WUST-AAE-NMaO-Report.sty'.