From e38f28a0b5ca83d212cd2881e16fa77d93944857 Mon Sep 17 00:00:00 2001
From: Gromiusz <kgorka22@gmail.com>
Date: Mon, 13 Jan 2025 17:32:54 +0100
Subject: [PATCH] fix: remove numpy

---
 code/threads_indep.py             | 24 ++++----
 code/threads_indep_numba_numpy.py | 91 +++++++++++++++++++++++++++++++
 2 files changed, 101 insertions(+), 14 deletions(-)
 create mode 100644 code/threads_indep_numba_numpy.py

diff --git a/code/threads_indep.py b/code/threads_indep.py
index deae5381..1a4fd1b2 100644
--- a/code/threads_indep.py
+++ b/code/threads_indep.py
@@ -1,16 +1,14 @@
 import gc
 import time
-import numpy as np
 from numba import njit, prange
 from time_measurement import time_measurement_longest, longest_threads_time_accumulator, tests_time
 import linear_algebra_utils as linAlg
 
-# Funkcje równoległe z Numba
+
 @njit(parallel=True)
 def numba_matrix_vector_multiply(A, input_x, Ax):
-    n, m = A.shape
-    for i in prange(n):
-        Ax[i] = np.dot(A[i], input_x)
+    for i in prange(len(A)):
+        Ax[i] = sum(A[i][j] * input_x[j] for j in range(len(input_x)))
 
 @njit(parallel=True)
 def numba_vector_vector_subtraction(b, Ax, residual):
@@ -51,32 +49,30 @@ def RichardsonMethodThreads(A, b, lambda_min, lambda_max, max_iterations, x0=Non
     gc.disable()
     start_time = time.perf_counter()
 
-    A = np.array(A)
-    b = np.array(b)
-    x0 = np.array(x0) if x0 is not None else np.zeros_like(b)
-    x = x0.copy()
+    x0 = x0 if x0 is not None else [0.0] * len(b)
+    x = x0[:]
 
     omega = 2 / (lambda_min + lambda_max)
     n = len(b)
 
     for iteration in range(max_iterations):
-        Ax = np.zeros_like(x)
+        Ax = [0.0] * n
         matrix_vector_multiply(A, x, Ax)
         longest_threads_time_accumulator.save_lap_and_reset()
 
-        residual = np.zeros_like(b)
+        residual = [0.0] * n
         vector_vector_subtraction(b, Ax, residual)
         longest_threads_time_accumulator.save_lap_and_reset()
 
-        change_vector = np.zeros_like(residual)
+        change_vector = [0.0] * n
         scalar_vector_multiply(omega, residual, change_vector)
         longest_threads_time_accumulator.save_lap_and_reset()
 
-        _x = np.zeros_like(x)
+        _x = [0.0] * n
         vector_vector_addition(x, change_vector, _x)
         longest_threads_time_accumulator.save_lap_and_reset()
 
-        x = _x.copy()
+        x = _x[:]
 
         if linAlg.SequentialLinearAlgebraUtils.vector_norm(residual) < tol:
             break
diff --git a/code/threads_indep_numba_numpy.py b/code/threads_indep_numba_numpy.py
new file mode 100644
index 00000000..deae5381
--- /dev/null
+++ b/code/threads_indep_numba_numpy.py
@@ -0,0 +1,91 @@
+import gc
+import time
+import numpy as np
+from numba import njit, prange
+from time_measurement import time_measurement_longest, longest_threads_time_accumulator, tests_time
+import linear_algebra_utils as linAlg
+
+# Funkcje równoległe z Numba
+@njit(parallel=True)
+def numba_matrix_vector_multiply(A, input_x, Ax):
+    n, m = A.shape
+    for i in prange(n):
+        Ax[i] = np.dot(A[i], input_x)
+
+@njit(parallel=True)
+def numba_vector_vector_subtraction(b, Ax, residual):
+    for i in prange(len(b)):
+        residual[i] = b[i] - Ax[i]
+
+@njit(parallel=True)
+def numba_scalar_vector_multiply(omega, vector, result):
+    for i in prange(len(vector)):
+        result[i] = omega * vector[i]
+
+@njit(parallel=True)
+def numba_vector_vector_addition(input_x, vector, output_x):
+    for i in prange(len(input_x)):
+        output_x[i] = input_x[i] + vector[i]
+
+# Funkcje z dekoratorem
+@time_measurement_longest(longest_threads_time_accumulator)
+def matrix_vector_multiply(A, input_x, Ax):
+    numba_matrix_vector_multiply(A, input_x, Ax)
+
+@time_measurement_longest(longest_threads_time_accumulator)
+def vector_vector_subtraction(b, Ax, residual):
+    numba_vector_vector_subtraction(b, Ax, residual)
+
+@time_measurement_longest(longest_threads_time_accumulator)
+def scalar_vector_multiply(omega, vector, result):
+    numba_scalar_vector_multiply(omega, vector, result)
+
+@time_measurement_longest(longest_threads_time_accumulator)
+def vector_vector_addition(input_x, vector, output_x):
+    numba_vector_vector_addition(input_x, vector, output_x)
+
+# Metoda Richardson z obsługą wątków
+def RichardsonMethodThreads(A, b, lambda_min, lambda_max, max_iterations, x0=None, tol=1e-5):
+    longest_threads_time_accumulator.hard_reset()
+
+    gc.disable()
+    start_time = time.perf_counter()
+
+    A = np.array(A)
+    b = np.array(b)
+    x0 = np.array(x0) if x0 is not None else np.zeros_like(b)
+    x = x0.copy()
+
+    omega = 2 / (lambda_min + lambda_max)
+    n = len(b)
+
+    for iteration in range(max_iterations):
+        Ax = np.zeros_like(x)
+        matrix_vector_multiply(A, x, Ax)
+        longest_threads_time_accumulator.save_lap_and_reset()
+
+        residual = np.zeros_like(b)
+        vector_vector_subtraction(b, Ax, residual)
+        longest_threads_time_accumulator.save_lap_and_reset()
+
+        change_vector = np.zeros_like(residual)
+        scalar_vector_multiply(omega, residual, change_vector)
+        longest_threads_time_accumulator.save_lap_and_reset()
+
+        _x = np.zeros_like(x)
+        vector_vector_addition(x, change_vector, _x)
+        longest_threads_time_accumulator.save_lap_and_reset()
+
+        x = _x.copy()
+
+        if linAlg.SequentialLinearAlgebraUtils.vector_norm(residual) < tol:
+            break
+
+    end_time = time.perf_counter()
+    gc.enable()
+    total_time = end_time - start_time
+    sequential_time = total_time - longest_threads_time_accumulator.total_time
+
+    print(f"Total: {total_time:.3e}s, Seq: {sequential_time:.3e}s, Parallel (threads): {longest_threads_time_accumulator.total_time:.3e}s, Tests time: {tests_time.total_time:.3e}s")
+
+    return x, 0