Merge pull request #7 from kuhyx/ola

two more matrices to test, processes implementation

code/eigenvalue_methods.py

@@ -1,7 +1,10 @@
+import numpy as np 
 class EigenvalueMethods:
     @staticmethod
     def power_method(LinAlgType, A, max_iter, tol=1e-6):
-        n = len(A)
+        if isinstance(A, list): #słabe, szkoda czasu, trzeba przypilnować, żeby od razu każda macierz była tego samego typu
+            A = np.array(A)
+        n = A.shape[0]
         x = [1] * n
         lambda_old = 0
         
@@ -17,9 +20,17 @@ class EigenvalueMethods:
 
     @staticmethod
     def inverse_power_method(LinAlgType, A, max_iter, tol=1e-6):
-        n = len(A)
+        import scipy
+        if scipy.sparse.issparse(A):
+            A = A.toarray()  # Convert sparse matrix to dense array
+
+        if isinstance(A, list):
+            A = np.array(A)  # Convert list to NumPy array if needed
+        n =  A.shape[0]
         I = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
+
         A_inv = [LinAlgType.gaussian_elimination(A.tolist(), I_col) for I_col in I]
+
         A_inv = list(map(list, zip(*A_inv)))
 
         return 1 / EigenvalueMethods.power_method(LinAlgType, A_inv, max_iter, tol)

code/linear_algebra_utils.py

@@ -1,8 +1,11 @@
 import math
+import itertools
+import operator
+from multiprocessing import Pool
 from abc import ABC, abstractmethod
 from concurrent.futures import ThreadPoolExecutor
 from functools import partial
-from time_measurement import time_measurement, threads_time_accumulator
+from time_measurement import time_measurement, time_accumulator
 
 class LinearAlgebraUtils(ABC):
     @staticmethod
@@ -72,7 +75,7 @@ class SequentialLinearAlgebraUtils(ABC):
     
     @staticmethod
     def vector_norm(v):
-        return sum(x*x for x in v)**0.5
+        return math.sqrt(sum(x*x for x in v))
 
     @staticmethod
     def vector_scalar_divide(x, scalar):
@@ -171,7 +174,7 @@ class ThreadsLinearAlgebraUtils(ABC):
 
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def dot_product(v1, v2):
         chunks = ThreadsLinearAlgebraUtils.divide_vectors_to_chunks(v1, v2)
         with ThreadPoolExecutor(max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS) as executor:
@@ -179,7 +182,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return sum(results)
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def matrix_vector_multiply(A, x):
         chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(A)
         with ThreadPoolExecutor(max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS) as executor:
@@ -188,7 +191,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return [item for sublist in results for item in sublist]
     
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def vector_norm(v):
         chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(v)
 
@@ -201,7 +204,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return total_sum**0.5
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def vector_scalar_divide(x, scalar):
         chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(x)
 
@@ -210,7 +213,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return [item for sublist in results for item in sublist]
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def matrix_scalar_multiply(A, w):
         chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(A)
         with ThreadPoolExecutor(max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS) as executor:
@@ -218,7 +221,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return [item for sublist in results for item in sublist]
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def vector_vector_subtraction(v1, v2):
         chunks = ThreadsLinearAlgebraUtils.divide_vectors_to_chunks(v1, v2)
         with ThreadPoolExecutor(max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS) as executor:
@@ -227,7 +230,7 @@ class ThreadsLinearAlgebraUtils(ABC):
 
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def vector_vector_addition(v1, v2):
         chunks = ThreadsLinearAlgebraUtils.divide_vectors_to_chunks(v1, v2)
         with ThreadPoolExecutor(max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS) as executor:
@@ -235,7 +238,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return [item for sublist in results for item in sublist]
     
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def scalar_vector_multiply(omega, vector):
         chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(vector)
         with ThreadPoolExecutor(max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS) as executor:
@@ -244,7 +247,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return [item for sublist in results for item in sublist]
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def matrix_norm(A):
         chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(A)
 
@@ -258,13 +261,13 @@ class ThreadsLinearAlgebraUtils(ABC):
         return math.sqrt(total_sum)
     
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def divide_matrixes_to_chunks(A, B):
         chunk_size = len(A) // ThreadsLinearAlgebraUtils.NUM_THREADS
         return [(A[i:i + chunk_size], B[i:i + chunk_size]) for i in range(0, len(A), chunk_size)]
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def matrix_matrix_subtraction(A, B):
 
         def subtract_chunk(pair):
@@ -277,7 +280,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return [row for chunk in results for row in chunk]
 
     @staticmethod
-    @time_measurement(threads_time_accumulator)
+    @time_measurement(time_accumulator)
     def gaussian_elimination(A, b):
         n = len(A)
         M = [row[:] for row in A]
@@ -310,3 +313,161 @@ class ThreadsLinearAlgebraUtils(ABC):
                 M[k][-1] -= M[k][i] * x[i]
 
         return x
+
+@time_measurement(time_accumulator)
+def process_row(params):
+    A, k, i = params
+    factor = A[i][k] / A[k][k]
+    return [A[i][j] - factor * A[k][j] for j in range(len(A[0]))]
+
+@time_measurement(time_accumulator)
+def divide_by_scalar(pair):
+    xi, scalar = pair
+    return xi / scalar
+
+@time_measurement(time_accumulator)
+def multiply_by_scalar(pair):
+    element, scalar = pair
+    return element * scalar
+
+class ProcessLinearAlgebraUtils:
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def dot_product(v1, v2):
+        with Pool() as pool:
+            result = pool.starmap(ProcessLinearAlgebraUtils.multiply_elements, zip(v1, v2))
+        return sum(result)
+    
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def multiply_elements(x, y):
+        return x * y
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def matrix_vector_multiply_row(params):
+        row, vector = params
+        return SequentialLinearAlgebraUtils.dot_product(row, vector)
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def matrix_vector_multiply(A, x):
+        with Pool() as pool:
+            result = pool.map(ProcessLinearAlgebraUtils.matrix_vector_multiply_row, [(row, x) for row in A])
+        return list(result)
+    
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def vector_norm(v):
+        with Pool() as pool:
+            squared = pool.map(ProcessLinearAlgebraUtils.square, v)
+        return math.sqrt(sum(squared))
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def square(x):
+        return x * x
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def vector_scalar_divide(x, scalar):
+        with Pool() as pool:
+            result = pool.map(divide_by_scalar, [(xi, scalar) for xi in x])
+        return list(result)
+    
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def divide_vector_by_scalar(x, scalar):
+        with Pool() as pool:
+            result = pool.map(ProcessLinearAlgebraUtils.vector_scalar_divide, [(xi, scalar) for xi in x])
+        return list(result)
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def matrix_scalar_multiply_row(params):
+        row, w = params
+        return [w * element for element in row]
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def matrix_scalar_multiply(A, w):
+        with Pool() as pool:
+            result = pool.map(ProcessLinearAlgebraUtils.matrix_scalar_multiply_row, [(row, w) for row in A])
+        return result
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def vector_vector_operation(params):
+        v1, v2, op = params
+        return op(v1, v2)
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def vector_vector_subtraction(v1, v2):
+        with Pool() as pool:
+            result = pool.map(ProcessLinearAlgebraUtils.vector_vector_operation, zip(v1, v2, itertools.repeat(operator.sub)))
+        return list(result)
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def vector_vector_addition(v1, v2):
+        with Pool() as pool:
+            result = pool.map(ProcessLinearAlgebraUtils.vector_vector_operation, zip(v1, v2, itertools.repeat(operator.add)))
+        return list(result)
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def scalar_vector_multiply(omega, vector):
+        with Pool() as pool:
+            result = pool.map(multiply_by_scalar, [(element, omega) for element in vector])
+        return list(result)
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def matrix_norm(A):
+        with Pool() as pool:
+            row_sums = pool.map(lambda row: sum(x ** 2 for x in row), A)
+        return math.sqrt(sum(row_sums))
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def matrix_matrix_subtraction(A, B):
+        def subtract_rows(row_pair):
+            return [a - b for a, b in zip(*row_pair)]
+
+        with Pool() as pool:
+            result = pool.starmap(subtract_rows, zip(A, B))
+        return result
+
+    @staticmethod
+    @time_measurement(time_accumulator)
+    def gaussian_elimination(A, b):
+        try:
+            n = len(A)
+            A = [list(row) + [b_i] for row, b_i in zip(A, b)] 
+
+            for k in range(n):
+                # Pivoting
+                max_index = max(range(k, n), key=lambda x: abs(A[x][k]))
+                if A[max_index][k] == 0:
+                    raise ValueError("Matrix is singular and cannot be solved.")
+                A[k], A[max_index] = A[max_index], A[k]
+
+                # Parallel row processing
+                with Pool() as pool:
+                    results = pool.map(process_row, [(A, k, i) for i in range(k + 1, n)])
+
+                # Update remaining rows in matrix
+                for i in range(k + 1, n):
+                    A[i] = results[i - (k + 1)]
+
+            # Back substitution
+            x = [0] * n
+            for i in range(n - 1, -1, -1):
+                sum_ax = sum(A[i][j] * x[j] for j in range(i + 1, n))
+                x[i] = (A[i][-1] - sum_ax) / A[i][i]
+
+            return x
+        except Exception as e:
+            print(f"Error during Gaussian elimination: {e}")
+            return None

code/matrix_generator.py

@@ -1,4 +1,5 @@
 import numpy as np
+import scipy.io
 
 class MatrixGenerator:
     @staticmethod
@@ -17,10 +18,39 @@ class MatrixGenerator:
         return spd_matrix
 
     @staticmethod
-    def generate_random_matrix_and_vector(size):
-        A = MatrixGenerator.generate_spd_matrix(size)
-        b = np.random.uniform(-1, 1, size)
-        return A, b
-
     def generate_identity_matrix(size):
         return np.eye(size)
+    
+    @staticmethod
+    def generate_alternating_vector(size):
+        return np.tile([1, 2], int(np.ceil(size / 2)))[:size]
+    
+    @staticmethod
+    def get_matrix_from_file(file_path, problem):
+        mat_contents = scipy.io.loadmat(file_path)
+        problem_record = mat_contents['Problem'][0][0]
+        A = np.array(problem_record[problem])
+        if scipy.sparse.issparse(A):
+            A = A.toarray()
+        return A
+
+    @staticmethod
+    def generate_matrix_and_vector(type, size=None):
+        if type == 'spd':
+            if size is None:
+                raise ValueError("Size must be provided for SPD matrix generation.")
+            matrix = MatrixGenerator.generate_spd_matrix(size)
+            vector = np.random.uniform(-1, 1, size)
+        elif type == 'nemeth12':
+            matrix = -1 * MatrixGenerator.get_matrix_from_file("nemeth12.mat", 1)
+            size = matrix.shape[0]
+            vector = MatrixGenerator.generate_alternating_vector(size)
+        elif type == 'poli3':
+            matrix = MatrixGenerator.get_matrix_from_file("poli3.mat", 2)
+            size = matrix.shape[0]
+            vector = MatrixGenerator.generate_alternating_vector(size)
+        else:
+            raise ValueError("Invalid type specified. Choose 'spd', 'nemeth12', or 'poli3'.")
+        
+        return matrix, vector
+

code/processing_type.py

@@ -3,3 +3,4 @@ from enum import Enum, auto
 class ProcessingType(Enum):
     SEQUENTIAL = auto()
     THREADS = auto()
+    PROCESSES = auto()

code/richardson_method.py

@@ -1,14 +1,13 @@
-from linear_algebra_utils import LinearAlgebraUtils
-from linear_algebra_utils import SequentialLinearAlgebraUtils
-from linear_algebra_utils import ThreadsLinearAlgebraUtils
+import linear_algebra_utils as linAlg
 from eigenvalue_methods import EigenvalueMethods
 from matrix_generator import MatrixGenerator
 from processing_type import ProcessingType
-from time_measurement import threads_time_accumulator
+from time_measurement import time_measurement, time_accumulator, tests_time
 import time
 import gc
 
 class RichardsonMethod:
+    @time_measurement(time_accumulator)
     def __init__(self, method: ProcessingType, A, b, max_iterations, size: int, x0=None, tol=1e-5):
         self.LinAlg = self.assign_LinAlgType(method)
         self.A = A
@@ -16,7 +15,7 @@ class RichardsonMethod:
         self.x0 = x0 if x0 is not None else [0.0] * len(b)
         self.max_iterations = max_iterations
         self.tol = tol
-        self.I = MatrixGenerator.generate_identity_matrix(size)
+        # self.I = MatrixGenerator.generate_identity_matrix(size)
         self.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.LinAlg, self.A, max_iterations)
         if self.lambda_min < 0:
             raise ValueError("Matrix A is not positive semi-definite.")
@@ -39,20 +38,21 @@ class RichardsonMethod:
     
     @staticmethod
     def assign_LinAlgType(method):
-        metody = {
-            ProcessingType.SEQUENTIAL: SequentialLinearAlgebraUtils,
-            ProcessingType.THREADS: ThreadsLinearAlgebraUtils
+        methods = {
+            ProcessingType.SEQUENTIAL: linAlg.SequentialLinearAlgebraUtils,
+            ProcessingType.THREADS: linAlg.ThreadsLinearAlgebraUtils,
+            ProcessingType.PROCESSES: linAlg.ProcessLinearAlgebraUtils
         }
         
         try:
-            return metody[method]
+            return methods[method]
         except KeyError:
-            raise ValueError("Unknown method, please use 'SEQUENTIAL' or 'THREADS'.")
+            raise ValueError("Unknown method, please use 'SEQUENTIAL', 'THREADS' or 'PROCESSES'.")
 
     def solve(self):
         gc.disable()
-        threads_time_accumulator.total_time = 0
-        start = time.time()
+        time_accumulator.total_time = 0
+        start = time.perf_counter()
         x = self.x0[:]
         #if RichardsonMethod.convergence_norm(self.LinAlg, self.A, self.omega, self.I) >= 1:
         #    return RichardsonMethod.convergence_norm(self.A, self.omega, self.I), "Richardson method for those values will NOT converge", 
@@ -61,14 +61,23 @@ class RichardsonMethod:
             Ax = self.LinAlg.matrix_vector_multiply(self.A, x)
             residual = self.LinAlg.vector_vector_subtraction(self.b, Ax)
             x = self.LinAlg.vector_vector_addition(x, self.LinAlg.scalar_vector_multiply(self.omega, residual))
+            if (linAlg.SequentialLinearAlgebraUtils.vector_norm(residual) < self.tol):
+                break
 
-        end = time.time()
+        end = time.perf_counter()
         total_time = end - start
         gc.enable()
-        if(self.LinAlg == SequentialLinearAlgebraUtils):
-            print(f"Total: {total_time:.3e}s")
-        elif(self.LinAlg == ThreadsLinearAlgebraUtils):
-            sequential_time = total_time - threads_time_accumulator.total_time
-            print(f"Total: {total_time:.3e}s, Seq: {sequential_time:.3e}s, Parallel: {threads_time_accumulator.total_time:.3e}s")
+
+        match self.LinAlg:
+            case linAlg.SequentialLinearAlgebraUtils:
+                print(f"Total: {total_time:.3e}s")
+            case linAlg.ThreadsLinearAlgebraUtils:
+                sequential_time = total_time - time_accumulator.total_time
+                print(f"Total: {total_time:.3e}s, Seq: {sequential_time:.3e}s, Parallel (threads): {time_accumulator.total_time:.3e}s, Tests time: {tests_time.total_time:.3e}s")
+            case linAlg.ProcessLinearAlgebraUtils:
+                sequential_time = total_time - time_accumulator.total_time
+                print(f"Total: {total_time:.3e}s, Seq: {sequential_time:.3e}s, Parallel (processes): {time_accumulator.total_time:.3e}s, Tests time: {tests_time.total_time:.3e}s")
+            case _:
+                print("Unhandled LinAlg type")
             
         return x, 0

code/tests.py

@@ -4,6 +4,7 @@ from scipy.sparse.linalg import cg
 from matrix_generator import MatrixGenerator
 from richardson_method import RichardsonMethod
 from processing_type import ProcessingType
+from time_measurement import time_measurement, tests_time
 
 def calculate_norm_numpy(I, w, A):
     # Calculate the difference between I and w * A
@@ -31,13 +32,27 @@ def calcualte_norm_from_matrix_numpy(A, n, max_iterations):
     return calculate_norm_numpy(I, omega, A)
 
 @pytest.mark.parametrize("n", [2, 3, 4, 5, 10, 20, 50, 100])
-@pytest.mark.parametrize("processing_type", [ProcessingType.SEQUENTIAL, ProcessingType.THREADS])
-def test_richardson_vs_cg(n: int, processing_type: ProcessingType, capsys):
-    print("matrix size: ", n)
+@pytest.mark.parametrize("processing_type", [ProcessingType.SEQUENTIAL, ProcessingType.THREADS, ProcessingType.PROCESSES])
+@pytest.mark.parametrize("matrix_type", ["spd", "nemeth12", "poli3"])
+@time_measurement(tests_time)
+def test_richardson_vs_cg(n: int, processing_type: ProcessingType, matrix_type: str, capsys):
+    print("matrix type: ", matrix_type)
+    print("matrix size: ", n if matrix_type == "spd" else "fixed")
     tolerance = 1e-5
     max_iterations=1000
-    A, b = MatrixGenerator.generate_random_matrix_and_vector(n)
-    richardson_solver = RichardsonMethod(processing_type, A, b, max_iterations, size=n, tol=1e-7)
+    if matrix_type in ["nemeth12", "poli3"] and n != 2:
+        pytest.skip("Fixed matrix size for nemeth12 and poli3, skipping redundant runs.")
+    
+    if matrix_type == "spd":
+        A, b = MatrixGenerator.generate_matrix_and_vector('spd', size=n)
+    elif matrix_type == "poli3":
+        A, b = MatrixGenerator.generate_matrix_and_vector('poli3')
+    elif matrix_type == "nemeth12":
+        A, b = MatrixGenerator.generate_matrix_and_vector('nemeth12')
+    else:
+        raise ValueError("Invalid matrix type specified. Choose 'spd', 'poli3', or 'nemeth12'.")
+    
+    richardson_solver = RichardsonMethod(processing_type, A, b, max_iterations, size=A.shape[0], tol=1e-7)
     # solution_richardson, info_richardson = richardson_solver.solve()
 
     solution_richardson, info_richardson = None, None
@@ -48,7 +63,7 @@ def test_richardson_vs_cg(n: int, processing_type: ProcessingType, capsys):
     captured = capsys.readouterr()
     print("Captured output:", captured.out)
     
-    solution_cg, info = cg(A, b)
+    solution_cg, info = cg(A, b, atol=0.)
     
     if info == 0:  # SciPy CG converged
         assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n)

code/time_measurement.py

@@ -5,15 +5,16 @@ class TimeAccumulator:
     def __init__(self):
         self.total_time = 0
 
-threads_time_accumulator = TimeAccumulator()
+time_accumulator = TimeAccumulator()
+tests_time = TimeAccumulator()
 
 def time_measurement(accumulator):
     def decorator(func):
         @wraps(func)
         def inner(*args, **kwargs):
-            start = time.time()
+            start = time.perf_counter()
             result = func(*args, **kwargs)
-            end = time.time()
+            end = time.perf_counter()
             accumulator.total_time += end - start
             return result
         return inner