diff --git a/benchmarks/gaussian_elimination_and_lu_decomposition_benchmark.py b/benchmarks/gaussian_elimination_and_lu_decomposition_benchmark.py
new file mode 100644
index 0000000..c5ad086
--- /dev/null
+++ b/benchmarks/gaussian_elimination_and_lu_decomposition_benchmark.py
@@ -0,0 +1,184 @@
+from __future__ import annotations
+
+import argparse
+import os
+import random
+import time
+from collections.abc import Callable
+
+import matplotlib.pyplot as plt
+
+from matrix_calculus.lu_decomposition import LUDecomposition
+from matrix_calculus.matrix_processor import MatrixProcessor
+
+
+def generate_matrix(n: int, seed: int = 0) -> list[list[float]]:
+    random.seed(seed)
+
+    A = [[random.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
+
+    for i in range(n):
+        A[i][i] += n  # stabilność
+
+    return A
+
+
+def measure_time(func: Callable, *args) -> float:
+    start = time.perf_counter()
+    func(*args)
+    end = time.perf_counter()
+    return end - start
+
+
+def benchmark_gaussian_elimination(out_dir: str, n_max: int):
+    sizes = list(range(10, n_max + 1, 10))
+
+    times_no_pivot = []
+    times_pivot = []
+
+    for n in sizes:
+        M = generate_matrix(n, seed=42)
+
+        t1 = measure_time(MatrixProcessor.gaussian_elimination, M, False)
+        t2 = measure_time(MatrixProcessor.gaussian_elimination, M, True)
+
+        times_no_pivot.append(t1)
+        times_pivot.append(t2)
+
+        print(f"[Gauss] n={n} done")
+
+    path = os.path.join(out_dir, "gauss_benchmark.png")
+
+    plt.figure()
+    plt.plot(sizes, times_no_pivot, label="bez pivotingu")
+    plt.plot(sizes, times_pivot, label="z pivotingiem")
+    plt.xlabel("Rozmiar macierzy (n)")
+    plt.ylabel("Czas [s]")
+    plt.title("Eliminacja Gaussa – porównanie")
+    plt.legend()
+    plt.grid()
+    plt.savefig(path)
+
+    print(f"[SAVE] {path}")
+
+
+def benchmark_lu_decomposition(out_dir: str, n_max: int):
+    sizes = list(range(10, n_max + 1, 10))
+
+    times_no_pivot = []
+    times_pivot = []
+
+    for n in sizes:
+        M = generate_matrix(n, seed=123)
+
+        t1 = measure_time(LUDecomposition.decompose, M, False)
+        t2 = measure_time(LUDecomposition.decompose, M, True)
+
+        times_no_pivot.append(t1)
+        times_pivot.append(t2)
+
+        print(f"[LU] n={n} done")
+
+    path = os.path.join(out_dir, "lu_benchmark.png")
+
+    plt.figure()
+    plt.plot(sizes, times_no_pivot, label="bez pivotingu")
+    plt.plot(sizes, times_pivot, label="z pivotingiem")
+    plt.xlabel("Rozmiar macierzy (n)")
+    plt.ylabel("Czas [s]")
+    plt.title("LU faktoryzacja – porównanie")
+    plt.legend()
+    plt.grid()
+    plt.savefig(path)
+
+    print(f"[SAVE] {path}")
+
+
+def print_example_results():
+    M = [
+        [2.0, 1.0, 1.0],
+        [4.0, -6.0, 0.0],
+        [-2.0, 7.0, 2.0],
+    ]
+
+    print("\n====================")
+    print("MACIERZ WEJŚCIOWA")
+    print("====================")
+    for row in M:
+        print(row)
+
+    print("\n====================")
+    print("GAUSS BEZ PIVOTINGU")
+    print("====================")
+    result = MatrixProcessor.gaussian_elimination(M, pivoting=False)
+    for row in result:
+        print(row)
+
+    print("\n====================")
+    print("GAUSS Z PIVOTINGIEM (bez normalizacji przekątnej)")
+    print("====================")
+    result = MatrixProcessor.gaussian_elimination(
+        M,
+        pivoting=True,
+        normalize_diagonal=False,
+    )
+    for row in result:
+        print(row)
+
+    print("\n====================")
+    print("LU BEZ PIVOTINGU")
+    print("====================")
+    L, U = LUDecomposition.decompose(M, pivoting=False)
+    print("L:")
+    for row in L:
+        print(row)
+    print("U:")
+    for row in U:
+        print(row)
+
+    print("\n====================")
+    print("LU Z PIVOTINGIEM")
+    print("====================")
+    P, L, U = LUDecomposition.decompose(M, pivoting=True)
+    print("P:")
+    for row in P:
+        print(row)
+    print("L:")
+    for row in L:
+        print(row)
+    print("U:")
+    for row in U:
+        print(row)
+
+
+def benchmark_gaussian_elimination_and_lu_decomposition(out_dir: str, n_max: int):
+    os.makedirs(out_dir, exist_ok=True)
+
+    benchmark_gaussian_elimination(out_dir, n_max)
+    benchmark_lu_decomposition(out_dir, n_max)
+    print_example_results()
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        "--output",
+        type=str,
+        default="results",
+        help="folder na wykresy",
+    )
+
+    parser.add_argument(
+        "--n",
+        type=int,
+        default=150,
+        help="maksymalny rozmiar macierzy",
+    )
+
+    args = parser.parse_args()
+
+    benchmark_gaussian_elimination_and_lu_decomposition(
+        out_dir=args.output,
+        n_max=args.n,
+    )
diff --git a/matrix_calculus/linalg/__init__.py b/matrix_calculus/linalg/__init__.py
new file mode 100644
index 0000000..e24b365
--- /dev/null
+++ b/matrix_calculus/linalg/__init__.py
@@ -0,0 +1,6 @@
+from __future__ import annotations
+
+from .matrix_properties import det
+from .matrix_properties import rank
+
+__all__ = ["det", "rank"]
diff --git a/matrix_calculus/linalg/matrix_properties.py b/matrix_calculus/linalg/matrix_properties.py
new file mode 100644
index 0000000..3056a74
--- /dev/null
+++ b/matrix_calculus/linalg/matrix_properties.py
@@ -0,0 +1,32 @@
+from __future__ import annotations
+
+from ..matrix_processor import MatrixProcessor
+
+
+def det(matrix: list[list[float]]) -> float:
+    """Calculate the determinant of a square matrix."""
+    processed_matrix, stats = MatrixProcessor.gaussian_elimination(
+        matrix,
+        stats=True,
+        normalize_diagonal=False,
+    )
+
+    if not stats["reversible"]:
+        return 0.0
+
+    determinant = (-1) ** stats["rows_swaps"]
+    for i in range(len(processed_matrix)):
+        determinant *= processed_matrix[i][i]
+
+    return determinant
+
+
+def rank(matrix: list[list[float]]) -> float:
+    """Calculate the rank of a matrix."""
+    _, stats = MatrixProcessor.gaussian_elimination(
+        matrix,
+        stats=True,
+        normalize_diagonal=False,
+    )
+
+    return stats["rank"]
diff --git a/matrix_calculus/matrix_processor.py b/matrix_calculus/matrix_processor.py
index b3f9652..17ddabe 100644
--- a/matrix_calculus/matrix_processor.py
+++ b/matrix_calculus/matrix_processor.py
@@ -1,12 +1,61 @@
 from __future__ import annotations
 
+from typing import Any
+from typing import Literal
+from typing import overload
+
 
 class MatrixProcessor:
+    @overload
+    @staticmethod
+    def gaussian_elimination(
+        matrix: list[list[float]],
+        pivoting: bool = ...,
+        stats: Literal[True] = ...,
+        normalize_diagonal: bool = ...,
+    ) -> tuple[list[list[float]], dict[str, Any]]:
+        pass
+
+    @overload
     @staticmethod
-    def gaussian_elimination(matrix: list[list[float]]) -> list[list[float]]:
+    def gaussian_elimination(
+        matrix: list[list[float]],
+        pivoting: bool = ...,
+        stats: Literal[False] = ...,
+        normalize_diagonal: bool = ...,
+    ) -> list[list[float]]:
+        pass
+
+    @staticmethod
+    def gaussian_elimination(
+        matrix: list[list[float]],
+        pivoting: bool = True,
+        stats: bool = False,
+        normalize_diagonal: bool = True,
+    ) -> list[list[float]] | tuple[list[list[float]], dict]:
+        if pivoting:
+            result = MatrixProcessor.gaussian_elimination_with_pivoting(
+                matrix,
+                normalize_diagonal,
+            )
+        else:
+            result = MatrixProcessor.gaussian_elimination_without_pivoting(
+                matrix,
+                normalize_diagonal,
+            )
+
+        return result[0] if not stats else result
+
+    @staticmethod
+    def gaussian_elimination_without_pivoting(
+        matrix: list[list[float]],
+        normalize_diagonal: bool = True,
+    ) -> tuple[list[list[float]], dict]:
         M = [row[:] for row in matrix]
         n = len(M)
         eps = 1e-12
+        rows_swaps = 0
+        rank = n
 
         def swap_rows(row1, row2):
             for col in range(n):
@@ -22,25 +71,37 @@ def swap_rows(row1, row2):
                         diagonal_element_value = M[row][diagonal_element]
 
                         swap_rows(diagonal_element, row)
+                        rows_swaps += 1
 
             if abs(diagonal_element_value) < eps:
-                raise ValueError("Matrix is singular")
+                rank -= 1
+                continue
 
-            for col in range(diagonal_element, n):
-                M[diagonal_element][col] /= diagonal_element_value
+            if normalize_diagonal:
+                for col in range(diagonal_element, n):
+                    M[diagonal_element][col] /= diagonal_element_value
 
             for row in range(diagonal_element + 1, n):
-                factor = M[row][diagonal_element]
+                factor = (
+                    M[row][diagonal_element]
+                    if normalize_diagonal
+                    else M[row][diagonal_element] / diagonal_element_value
+                )
                 for col in range(diagonal_element, n):
                     M[row][col] -= factor * M[diagonal_element][col]
 
-        return M
+        return M, {"rows_swaps": rows_swaps, "rank": rank, "reversible": rank == n}
 
     @staticmethod
-    def gaussian_elimination_pivoting(matrix: list[list[float]]) -> list[list[float]]:
+    def gaussian_elimination_with_pivoting(
+        matrix: list[list[float]],
+        normalize_diagonal: bool = True,
+    ) -> tuple[list[list[float]], dict]:
         M = [row[:] for row in matrix]
         n = len(M)
         eps = 1e-12
+        rows_swaps = 0
+        rank = n
 
         def swap_rows(row1, row2):
             for col in range(n):
@@ -57,14 +118,24 @@ def swap_rows(row1, row2):
 
             if highest_pivot_row != diagonal_element:
                 swap_rows(diagonal_element, highest_pivot_row)
+                rows_swaps += 1
 
             pivot_value = M[diagonal_element][diagonal_element]
             if abs(pivot_value) < eps:
-                raise ValueError("Matrix is singular")
+                rank -= 1
+                continue
+
+            if normalize_diagonal:
+                for col in range(diagonal_element, n):
+                    M[diagonal_element][col] /= pivot_value
 
             for row in range(diagonal_element + 1, n):
-                factor = M[row][diagonal_element] / pivot_value
+                factor = (
+                    M[row][diagonal_element]
+                    if normalize_diagonal
+                    else M[row][diagonal_element] / pivot_value
+                )
                 for col in range(diagonal_element, n):
                     M[row][col] -= factor * M[diagonal_element][col]
 
-        return M
+        return M, {"rows_swaps": rows_swaps, "rank": rank, "reversible": rank == n}
diff --git a/reports/imgs/gauss_benchmark.png b/reports/imgs/gauss_benchmark.png
new file mode 100644
index 0000000..d5f0e9c
Binary files /dev/null and b/reports/imgs/gauss_benchmark.png differ
diff --git a/reports/imgs/lu_benchmark.png b/reports/imgs/lu_benchmark.png
new file mode 100644
index 0000000..6e940da
Binary files /dev/null and b/reports/imgs/lu_benchmark.png differ
diff --git a/reports/matrix_gaussian_elimination_and_lu_decomposition_report.md b/reports/matrix_gaussian_elimination_and_lu_decomposition_report.md
new file mode 100644
index 0000000..15cfe0b
--- /dev/null
+++ b/reports/matrix_gaussian_elimination_and_lu_decomposition_report.md
@@ -0,0 +1,220 @@
+# Eliminacja Gaussa i LU faktoryzacja
+
+*Autorzy: Maja Byrecka, Michał Kowalczyk*
+
+## Eliminacja Gaussa bez pivotingu
+
+### Pseudokod
+
+```pseudo
+Wejście: macierz A ∈ ℝ^{n×n}
+
+Dla k = 0 do n-1:
+
+    pivot = A[k][k]
+
+    Jeśli |pivot| < ε:
+        Spróbuj znaleźć wiersz i > k taki, że |A[i][k]| > ε
+        Jeśli taki wiersz istnieje:
+            zamień wiersz k z i
+            pivot = A[k][k]
+        W przeciwnym razie:
+            zgłoś błąd (macierz osobliwa)
+
+    # normalizacja wiersza pivotowego
+    Dla j = k do n-1:
+        A[k][j] = A[k][j] / pivot
+
+    # eliminacja poniżej pivota
+    Dla i = k+1 do n-1:
+        factor = A[i][k]
+        Dla j = k do n-1:
+            A[i][j] = A[i][j] - factor * A[k][j]
+
+Wyjście: macierz w postaci trójkątnej górnej
+```
+
+### Złożoność
+
+- Czasowa: \( O(n^3) \)
+- Pamięciowa: \( O(n^2) \)
+  (tworzona jest kopia macierzy; możliwy wariant in-place)
+
+## Eliminacja Gaussa z pivotingiem
+
+### Pseudokod
+
+```pseudo
+Wejście: macierz A ∈ ℝ^{n×n}
+
+Dla k = 0 do n-1:
+
+    # wybór najlepszego pivota (partial pivoting)
+    max_row = k
+    max_value = |A[k][k]|
+
+    Dla i = k+1 do n-1:
+        Jeśli |A[i][k]| > max_value:
+            max_value = |A[i][k]|
+            max_row = i
+
+    Jeśli max_row ≠ k:
+        zamień wiersz k z max_row
+
+    pivot = A[k][k]
+
+    Jeśli |pivot| < ε:
+        zgłoś błąd (macierz osobliwa)
+
+    # eliminacja poniżej pivota
+    Dla i = k+1 do n-1:
+        factor = A[i][k] / pivot
+        Dla j = k do n-1:
+            A[i][j] = A[i][j] - factor * A[k][j]
+
+Wyjście: macierz w postaci trójkątnej górnej
+```
+
+### Przykładowe wyniki
+
+Dla macierzy wejściowej:
+`M = [[2.0, 1.0, 1.0], [4.0, -6.0, 0.0], [-2.0, 7.0, 2.0]]`
+
+- Eliminacja Gaussa bez pivotingu:
+`[[1.0, 0.5, 0.5], [0.0, 1.0, 0.25], [0.0, 0.0, 1.0]]`
+
+- Eliminacja Gasussa z pivotingiem (bez normalizacji przekątnej):
+`[[4.0, -6.0, 0.0], [0.0, 4.0, 1.0], [0.0, 0.0, 1.0]]`
+
+### Złożoność
+
+- Czasowa: \( O(n^3) \)
+- Pamięciowa: \( O(n^2) \)
+
+---
+
+## LU faktoryzacja bez pivotingu
+
+### Pseudokod
+
+```pseudo
+Wejście: macierz A ∈ ℝ^{n×n}
+
+Inicjalizuj:
+    L = macierz zerowa n×n
+    U = macierz zerowa n×n
+
+Dla i = 0 do n-1:
+
+    # wyznaczanie wiersza U
+    Dla k = i do n-1:
+        sum = 0
+        Dla j = 0 do i-1:
+            sum = sum + L[i][j] * U[j][k]
+        U[i][k] = A[i][k] - sum
+
+    L[i][i] = 1
+
+    # wyznaczanie kolumny L
+    Dla k = i+1 do n-1:
+        sum = 0
+        Dla j = 0 do i-1:
+            sum = sum + L[k][j] * U[j][i]
+
+        Jeśli U[i][i] == 0:
+            zgłoś błąd (dzielenie przez zero)
+
+        L[k][i] = (A[k][i] - sum) / U[i][i]
+
+Wyjście: macierze L (dolnotrójkątna) i U (górnotrójkątna)
+```
+
+### Złożoność
+
+- Czasowa: \( O(n^3) \)
+- Pamięciowa: \( O(n^2) \)
+
+---
+
+## LU faktoryzacja z pivotingiem
+
+### Pseudokod
+
+```pseudo
+Wejście: macierz A ∈ ℝ^{n×n}
+
+Inicjalizuj:
+    U = kopia A
+    L = macierz zerowa n×n
+    P = macierz jednostkowa n×n
+
+Dla k = 0 do n-1:
+
+    # wybór pivota (partial pivoting)
+    pivot = indeks i ≥ k maksymalizujący |U[i][k]|
+
+    Jeśli |U[pivot][k]| < ε:
+        zgłoś błąd (macierz osobliwa)
+
+    # zamiana wierszy
+    zamień wiersze k i pivot w U
+    zamień wiersze k i pivot w P
+
+    Jeśli k > 0:
+        zamień odpowiednie elementy w L (kolumny < k)
+
+    # eliminacja
+    Dla j = k+1 do n-1:
+        L[j][k] = U[j][k] / U[k][k]
+        Dla col = 0 do n-1:
+            U[j][col] = U[j][col] - L[j][k] * U[k][col]
+
+Dla i = 0 do n-1:
+    L[i][i] = 1
+
+Wyjście: macierze P, L, U takie, że PA = LU
+```
+
+### Przykładowe wyniki
+
+Dla macierzy wejściowej:
+`M = [[2.0, 1.0, 1.0], [4.0, -6.0, 0.0], [-2.0, 7.0, 2.0]]`
+
+- LU faktoryzacja bez pivotingu:
+`L = [[1.0, 0.0, 0.0], [2.0, 1.0, 0.0], [-1.0, -1.0, 1.0]]`
+`U = [[2.0, 1.0, 1.0], [0.0, -8.0, -2.0], [0.0, 0.0, 1.0]]`
+
+- LU faktoryzacja z pivotingiem:
+`P = [[0. 1. 0.], [1. 0. 0.], [0. 0. 1.]]`
+`L = [[1. 0. 0.], [0.5 1. 0.], [-0.5  1. 1.]]`
+`U = [[ 4. -6.  0.], [0. 4. 1.], [0. 0. 1.]]`
+
+### Złożoność
+
+- Czasowa: \( O(n^3) \)
+- Pamięciowa: \( O(n^2) \)
+
+## Porównanie czasów wykonywania algorytmów
+
+### Eliminacja Gaussa
+
+Obie wersje algorytmu mają taką samą złożoność obliczeniową, więc zgodnie z oczekiwaniami wykresy czasów ich wykonywania się pokrywają.
+
+<img src="imgs/gauss_benchmark.png" style="max-width:100%; height:auto;" />
+
+### LU faktoryzacja
+
+Tak samo jak w przypadku eliminacji Gaussa, obie wersje algorytmu faktoryzacji LU mają taką samą złożoność obliczeniową. Powodem rozbieżności jest różnica implementacji między nimi - wersja z pivotingiem wykonuje się znacznie szybciej ze względu na użycie biblioteki numerycznej przyspieszającej operacje na tablicach.
+
+<img src="imgs/lu_benchmark.png" style="max-width:100%; height:auto;" />
+
+
+
+---
+
+Wykonanie:
+
+* Język: Python3
+* Podział zadań:
+    * Faktoryzacja LU z pivotingiem i bez pivotingu - Maja Byrecka
+    * Eliminacja Gaussa z pivotingiem i bez pivotingu - Michał Kowalczyk
diff --git a/reports/resources/gaussian_elimination_and_lu_decomposition/code.txt b/reports/resources/gaussian_elimination_and_lu_decomposition/code.txt
new file mode 100644
index 0000000..4b327a7
--- /dev/null
+++ b/reports/resources/gaussian_elimination_and_lu_decomposition/code.txt
@@ -0,0 +1,380 @@
+
+========== lu_decomposition_test.py ==========
+from __future__ import annotations
+
+import numpy as np
+import pytest  # noqa: F401
+
+from matrix_calculus.lu_decomposition import LUDecomposition
+from matrix_calculus.matrix_multiplier import MatrixMultiplier
+
+
+def test_lu_no_pivot():
+    matrix = [[1, 2, 3], [5, 6, 7], [9, 10, 11]]
+    L, U = LUDecomposition.decompose(matrix=matrix, pivoting=False)
+    assert L == [[1, 0, 0], [5, 1, 0], [9, 2, 1]]
+    assert U == [[1, 2, 3], [0, -4, -8], [0, 0, 0]]
+
+
+def test_lu_with_pivot_singular():
+    # Test singular matrix (determinant equal to 0)
+    with pytest.raises(ZeroDivisionError):
+        matrix = [[1, 2, 3], [5, 6, 7], [9, 10, 11]]
+        LUDecomposition.decompose(matrix=matrix, pivoting=True)
+
+
+def test_lu_with_pivot_correct():
+    # Test correct decomposition
+    matrix = [[1, 2, 3], [4, 5, 6], [8, 10, 10]]
+    P, L, U = LUDecomposition.decompose(matrix=matrix, pivoting=True)
+    L == [[1, 0, 0], [4, 1, 0], [8, 2, 1]]
+    assert np.allclose(
+        MatrixMultiplier.traditional_multiplication(P, matrix),
+        MatrixMultiplier.traditional_multiplication(L, U),
+    )
+
+
+def test_lu_example():
+    date = 8 + 24
+    np.random.seed(date)
+    matrix = np.random.rand(date, date) * 10 // 1
+    # Make sure that generated matrix is not singular
+    matrix += np.eye(date) * 100
+
+    # LU without pivoting
+    L, U = LUDecomposition.decompose(matrix, pivoting=False)
+    assert np.allclose(
+        MatrixMultiplier.traditional_multiplication(L, U),
+        matrix,
+    )
+
+    # LU with pivoting
+    P, L, U = LUDecomposition.decompose(matrix, pivoting=True)
+    assert np.allclose(
+        MatrixMultiplier.traditional_multiplication(L, U),
+        MatrixMultiplier.traditional_multiplication(P, matrix),
+    )
+
+
+========== matrix_multiplier.py ==========
+from __future__ import annotations
+
+from matrix_calculus.utils import _next_power_of_2
+from matrix_calculus.utils import _pad_matrix
+from matrix_calculus.utils import _strassen
+from matrix_calculus.utils import _unpad_matrix
+
+
+class MatrixMultiplier:
+
+    @staticmethod
+    def traditional_multiplication(
+        matrix_1: list[list[float]],
+        matrix_2: list[list[float]],
+    ) -> list[list[float]]:
+
+        if len(matrix_1[0]) != len(matrix_2):
+            raise ValueError(
+                "Invalid matrices sizes. Given matrices cannot be multiplied.",
+            )
+
+        n_rows = len(matrix_1)
+        n_cols = len(matrix_2[0])
+
+        res = [[0.0 for _ in range(n_cols)] for _ in range(n_rows)]
+
+        for row in range(n_rows):
+            for col in range(n_cols):
+                for k in range(len(matrix_1[0])):
+                    res[row][col] += matrix_1[row][k] * matrix_2[k][col]
+
+        return res
+
+    @staticmethod
+    def strassen_matrix_multiplication(
+        matrix_1: list[list[float]],
+        matrix_2: list[list[float]],
+    ) -> list[list[float]]:
+
+        rows_1 = len(matrix_1)
+        cols_1 = len(matrix_1[0])
+        rows_2 = len(matrix_2)
+        cols_2 = len(matrix_2[0])
+
+        if not (rows_1 == cols_1 and rows_2 == cols_2 and rows_1 == rows_2):
+            raise ValueError(
+                "Both matrices must be squared matrices of the same size",
+            )
+
+        n = len(matrix_1)
+
+        if n == 1:
+            return [[matrix_1[0][0] * matrix_2[0][0]]]
+
+        m = _next_power_of_2(n)
+
+        if m == n:
+            return _strassen(matrix_1, matrix_2, n)
+
+        A_pad = _pad_matrix(matrix_1, m)
+        B_pad = _pad_matrix(matrix_2, m)
+
+        C_pad = _strassen(A_pad, B_pad, m)
+
+        return _unpad_matrix(C_pad, n)
+
+
+========== utils.py ==========
+from __future__ import annotations
+
+
+def _zeroes_matrix(n):
+    return [[0 for _ in range(n)] for _ in range(n)]
+
+
+# sub_matrix is a tuple of the form (matrix, row_start, col_start)
+
+
+def _add_sub_matrices(
+    sub_matrix_1,
+    sub_matrix_2,
+    n,
+    result,
+    result_matrix_offset=(0, 0),
+):
+    m1, r1, c1 = sub_matrix_1
+    m2, r2, c2 = sub_matrix_2
+    result_offset_row, result_offset_col = result_matrix_offset
+    for i in range(n):
+        for j in range(n):
+            result[result_offset_row + i][result_offset_col + j] = (
+                m1[r1 + i][c1 + j] + m2[r2 + i][c2 + j]
+            )
+    return result
+
+
+def _subtract_sub_matrices(
+    sub_matrix_1,
+    sub_matrix_2,
+    n,
+    result,
+    result_matrix_offset=(0, 0),
+):
+    m1, r1, c1 = sub_matrix_1
+    m2, r2, c2 = sub_matrix_2
+    result_offset_row, result_offset_col = result_matrix_offset
+    for i in range(n):
+        for j in range(n):
+            result[result_offset_row + i][result_offset_col + j] = (
+                m1[r1 + i][c1 + j] - m2[r2 + i][c2 + j]
+            )
+    return result
+
+
+def _set_sub_matrix(sub_matrix, other_matrix, n):
+    m, r, c = sub_matrix
+    for i in range(n):
+        for j in range(n):
+            m[i + r][j + c] = other_matrix[i][j]
+    return m
+
+
+def _sub_matrix_to_matrix(sub_matrix, n, result):
+    m, r, c = sub_matrix
+    for i in range(n):
+        for j in range(n):
+            result[i][j] = m[r + i][c + j]
+    return result
+
+
+def _add_matrices(*matrices, result=None):
+    assert result is not None
+    for i in range(len(result)):
+        for j in range(len(result[0])):
+            result[i][j] = matrices[0][i][j]
+            for k in range(1, len(matrices)):
+                result[i][j] += matrices[k][i][j]
+    return result
+
+
+def _subtract_matrix_from(target_matrix, matrix_to_subtract):
+    for i in range(len(target_matrix)):
+        for j in range(len(target_matrix[0])):
+            target_matrix[i][j] -= matrix_to_subtract[i][j]
+    return target_matrix
+
+
+def _strassen(matrix_1, matrix_2, n):
+    if n == 2:
+        c11 = matrix_1[0][0] * matrix_2[0][0] + matrix_1[0][1] * matrix_2[1][0]
+        c12 = matrix_1[0][0] * matrix_2[0][1] + matrix_1[0][1] * matrix_2[1][1]
+        c21 = matrix_1[1][0] * matrix_2[0][0] + matrix_1[1][1] * matrix_2[1][0]
+        c22 = matrix_1[1][0] * matrix_2[0][1] + matrix_1[1][1] * matrix_2[1][1]
+        return [[c11, c12], [c21, c22]]
+
+    else:
+        sub_matrix_dim = n >> 1
+
+        buffer_1 = _zeroes_matrix(sub_matrix_dim)
+        buffer_2 = _zeroes_matrix(sub_matrix_dim)
+
+        a11 = (matrix_1, 0, 0)
+        a12 = (matrix_1, 0, sub_matrix_dim)
+        a21 = (matrix_1, sub_matrix_dim, 0)
+        a22 = (matrix_1, sub_matrix_dim, sub_matrix_dim)
+
+        b11 = (matrix_2, 0, 0)
+        b12 = (matrix_2, 0, sub_matrix_dim)
+        b21 = (matrix_2, sub_matrix_dim, 0)
+        b22 = (matrix_2, sub_matrix_dim, sub_matrix_dim)
+
+        m1 = _strassen(
+            _add_sub_matrices(a11, a22, sub_matrix_dim, buffer_1),
+            _add_sub_matrices(b11, b22, sub_matrix_dim, buffer_2),
+            sub_matrix_dim,
+        )
+        m2 = _strassen(
+            _add_sub_matrices(a21, a22, sub_matrix_dim, buffer_1),
+            _sub_matrix_to_matrix(b11, sub_matrix_dim, buffer_2),
+            sub_matrix_dim,
+        )
+        m3 = _strassen(
+            _sub_matrix_to_matrix(a11, sub_matrix_dim, buffer_1),
+            _subtract_sub_matrices(b12, b22, sub_matrix_dim, buffer_2),
+            sub_matrix_dim,
+        )
+        m4 = _strassen(
+            _sub_matrix_to_matrix(a22, sub_matrix_dim, buffer_1),
+            _subtract_sub_matrices(b21, b11, sub_matrix_dim, buffer_2),
+            sub_matrix_dim,
+        )
+        m5 = _strassen(
+            _add_sub_matrices(a11, a12, sub_matrix_dim, buffer_1),
+            _sub_matrix_to_matrix(b22, sub_matrix_dim, buffer_2),
+            sub_matrix_dim,
+        )
+        m6 = _strassen(
+            _subtract_sub_matrices(a21, a11, sub_matrix_dim, buffer_1),
+            _add_sub_matrices(b11, b12, sub_matrix_dim, buffer_2),
+            sub_matrix_dim,
+        )
+        m7 = _strassen(
+            _subtract_sub_matrices(a12, a22, sub_matrix_dim, buffer_1),
+            _add_sub_matrices(b21, b22, sub_matrix_dim, buffer_2),
+            sub_matrix_dim,
+        )
+
+        c = _zeroes_matrix(n)
+
+        c11 = _add_matrices(m1, m4, m7, result=buffer_1)
+        _subtract_matrix_from(c11, m5)
+        _set_sub_matrix((c, 0, 0), c11, sub_matrix_dim)
+
+        c12 = _add_matrices(m3, m5, result=buffer_1)
+        _set_sub_matrix((c, 0, sub_matrix_dim), c12, sub_matrix_dim)
+
+        c21 = _add_matrices(m2, m4, result=buffer_1)
+        _set_sub_matrix((c, sub_matrix_dim, 0), c21, sub_matrix_dim)
+
+        c22 = _add_matrices(m1, m3, m6, result=buffer_1)
+        _subtract_matrix_from(c22, m2)
+        _set_sub_matrix(
+            (c, sub_matrix_dim, sub_matrix_dim),
+            c22,
+            sub_matrix_dim,
+        )
+
+        return c
+
+
+def _next_power_of_2(n):
+    return 1 if n == 0 else 2 ** ((n - 1).bit_length())
+
+
+def _pad_matrix(matrix, new_size):
+    old_size = len(matrix)
+    padded = [[0] * new_size for _ in range(new_size)]
+
+    for i in range(old_size):
+        for j in range(old_size):
+            padded[i][j] = matrix[i][j]
+
+    return padded
+
+
+def _unpad_matrix(matrix, size):
+    return [row[:size] for row in matrix[:size]]
+
+
+========== lu_decomposition.py ==========
+from __future__ import annotations
+
+import numpy as np
+
+
+class LUDecomposition:
+    @staticmethod
+    def decompose(
+        matrix: list[list[float]],
+        pivoting=False,
+    ) -> tuple:
+        if pivoting is True:
+            return LUDecomposition.decompose_with_pivot(matrix=matrix)
+        else:
+            return LUDecomposition.decompose_no_pivot(matrix=matrix)
+
+    @staticmethod
+    def decompose_no_pivot(
+        matrix: list[list[float]],
+    ) -> tuple[list[list[float]], list[list[float]]]:
+        n = len(matrix)
+        L = [[0.0] * n for _ in range(n)]
+        U = [[0.0] * n for _ in range(n)]
+
+        for i in range(n):
+            for k in range(i, n):
+                sum_u = sum(L[i][j] * U[j][k] for j in range(i))
+                U[i][k] = matrix[i][k] - sum_u
+            L[i][i] = 1.0
+            for k in range(i + 1, n):
+                sum_l = sum(L[k][j] * U[j][i] for j in range(i))
+                if U[i][i] == 0:
+                    raise ZeroDivisionError(
+                        "U diagonal element equals zero. Cannot decompose without pivoting.",
+                    )
+                L[k][i] = (matrix[k][i] - sum_l) / U[i][i]
+
+        return L, U
+
+    @staticmethod
+    def decompose_with_pivot(
+        matrix: list[list[float]],
+    ) -> tuple[list[list[float]], list[list[float]], list[list[float]]]:
+        A = np.array(matrix).copy().astype(float)
+        n = A.shape[0]
+
+        L = np.zeros((n, n))
+        U = A.copy()
+        P = np.eye(n)
+
+        for k in range(n):
+            pivot = np.argmax(np.abs(U[k:n, k])) + k
+            if U[pivot, k] < 1e-5:
+                raise ZeroDivisionError(
+                    "Matrix is singular. Cannot use LU decomposition with pivoting.",
+                )
+
+            # Switch rows
+            U[[k, pivot], :] = U[[pivot, k], :]
+            P[[k, pivot], :] = P[[pivot, k], :]
+            if k > 0:
+                L[[k, pivot], :k] = L[[pivot, k], :k]
+
+            for j in range(k + 1, n):
+                L[j, k] = U[j, k] / U[k, k]
+                U[j, :] = U[j, :] - L[j, k] * U[k, :]
+
+        for i in range(n):
+            L[i, i] = 1.0
+
+        return P, L, U
diff --git a/reports/resources/gaussian_elimination_and_lu_decomposition/report.pdf b/reports/resources/gaussian_elimination_and_lu_decomposition/report.pdf
new file mode 100644
index 0000000..273d571
Binary files /dev/null and b/reports/resources/gaussian_elimination_and_lu_decomposition/report.pdf differ
diff --git a/tests/gaussian_elimination_test.py b/tests/gaussian_elimination_test.py
index c78e090..99dce8e 100644
--- a/tests/gaussian_elimination_test.py
+++ b/tests/gaussian_elimination_test.py
@@ -6,44 +6,68 @@
 
 
 def test_gaussian_elimination():
+    print("\nTesting Gaussian elimination without pivoting...")
+
     # Test 1x1 matrix
     matrix = [[2]]
-    res = MatrixProcessor.gaussian_elimination(matrix)
+    res = MatrixProcessor.gaussian_elimination(matrix, pivoting=False)
     assert res == [[1.0]]
 
     # Test 2x2 matrix
     matrix = [[2, 4], [1, 3]]
-    res = MatrixProcessor.gaussian_elimination(matrix)
+    res = MatrixProcessor.gaussian_elimination(matrix, pivoting=False)
     assert res == [[1.0, 2.0], [0.0, 1.0]]
 
     # Test 3x3 singular matrix
-    with pytest.raises(ValueError):
-        matrix = [[2, 4, 6], [1, 3, 5], [0, 2, 4]]
-        MatrixProcessor.gaussian_elimination(matrix)
+    matrix = [[2, 4, 6], [1, 3, 5], [0, 2, 4]]
+    res, stats = MatrixProcessor.gaussian_elimination(
+        matrix,
+        pivoting=False,
+        stats=True,
+    )
+    print(res)
+    print(stats)
 
     # Test 3x3 matrix
     matrix = [[2, 4, 6], [1, 2, 5], [0, 2, 4]]
-    res = MatrixProcessor.gaussian_elimination(matrix)
+    res = MatrixProcessor.gaussian_elimination(matrix, pivoting=False)
     assert res == [[1.0, 2.0, 3.0], [0.0, 1.0, 2.0], [0.0, 0.0, 1.0]]
 
 
 def test_gaussian_elimination_pivoting():
+    print("\nTesting Gaussian elimination with pivoting...")
+
     # Test 1x1 matrix
     matrix = [[2]]
-    res = MatrixProcessor.gaussian_elimination_pivoting(matrix)
+    res = MatrixProcessor.gaussian_elimination(matrix, normalize_diagonal=False)
     assert res == [[2.0]]
 
     # Test 2x2 matrix
     matrix = [[2, 4], [1, 3]]
-    res = MatrixProcessor.gaussian_elimination_pivoting(matrix)
+    res = MatrixProcessor.gaussian_elimination(matrix, normalize_diagonal=False)
     assert res == [[2.0, 4.0], [0.0, 1.0]]
 
+    # Test 2x2 singular matrix
+    matrix = [[2, 4], [4, 8]]
+    res, stats = MatrixProcessor.gaussian_elimination(
+        matrix,
+        normalize_diagonal=False,
+        stats=True,
+    )
+    print(res)
+    print(stats)
+
     # Test 3x3 singular matrix
-    with pytest.raises(ValueError):
-        matrix = [[2, 4, 6], [1, 3, 5], [0, 2, 4]]
-        MatrixProcessor.gaussian_elimination_pivoting(matrix)
+    matrix = [[2, 4, 6], [1, 3, 5], [0, 2, 4]]
+    res, stats = MatrixProcessor.gaussian_elimination(
+        matrix,
+        stats=True,
+        normalize_diagonal=False,
+    )
+    print(res)
+    print(stats)
 
     # Test 3x3 matrix
     matrix = [[2, 4, 6], [1, 2, 5], [0, 2, 4]]
-    res = MatrixProcessor.gaussian_elimination_pivoting(matrix)
+    res = MatrixProcessor.gaussian_elimination(matrix, normalize_diagonal=False)
     assert res == [[2, 4, 6], [0, 2, 4], [0, 0, 2]]
diff --git a/tests/linalg_matrix_properties_test.py b/tests/linalg_matrix_properties_test.py
new file mode 100644
index 0000000..3defc6c
--- /dev/null
+++ b/tests/linalg_matrix_properties_test.py
@@ -0,0 +1,54 @@
+from __future__ import annotations
+
+import numpy as np
+import pytest  # noqa: F401
+
+import matrix_calculus.linalg as linalg
+
+
+def test_determinant():
+    # Test 1x1 matrix
+    matrix = [[5]]
+    assert linalg.det(matrix) == 5.0
+
+    # Test 2x2 matrix
+    matrix = [[1, 2], [3, 4]]
+    assert linalg.det(matrix) == -2.0
+
+    # Test 3x3 matrix
+    matrix = [[6, 1, 1], [4, -2, 5], [2, 8, 7]]
+    assert linalg.det(matrix) == np.linalg.det(matrix)
+
+    # Test 2x2 singular matrix
+    matrix = [[2, 4], [1, 2]]
+    assert linalg.det(matrix) == 0.0
+
+    # Test 3x3 singular matrix
+    matrix = [[2, 4, 6], [1, 2, 3], [1, 7, 5]]
+    assert linalg.det(matrix) == 0.0
+
+
+def test_rank():
+    # Test 1x1 matrix
+    matrix = [[5]]
+    assert linalg.rank(matrix) == 1
+
+    # Test 2x2 full rank matrix
+    matrix = [[1, 2], [3, 4]]
+    assert linalg.rank(matrix) == 2
+
+    # Test 3x3 full rank matrix
+    matrix = [[6, 1, 1], [4, -2, 5], [2, 8, 7]]
+    assert linalg.rank(matrix) == np.linalg.matrix_rank(matrix)
+
+    # Test 2x2 singular matrix
+    matrix = [[2, 4], [1, 2]]
+    assert linalg.rank(matrix) == 1
+
+    # Test 3x3 singular matrix
+    matrix = [[2, 4, 6], [1, 2, 3], [1, 7, 5]]
+    assert linalg.rank(matrix) == 2
+
+    # Test 3x3 singular matrix
+    matrix = [[2, 4, 6], [1, 2, 3], [4, 8, 12]]
+    assert linalg.rank(matrix) == 1