Reapply "Merge branch 'main' of https://github.com/kuhyx/PORR"

This reverts commit 87d6bb0d97.

code/eigenvalue_methods.py

@@ -1,43 +0,0 @@
-import numpy as np
-import scipy
-class EigenvalueMethods:
-    @staticmethod
-    def get_sing_vals(file_path):
-        mat_contents = scipy.io.loadmat(file_path)
-        A = mat_contents['S'][0][0]  # Pobranie pierwszego elementu z pola 'S'
-        singular_values = A['s'].flatten()
-        return singular_values
-     
-    @staticmethod
-    def power_method(LinAlgType, A, type, max_iter=100, tol=1e-6):
-        if (type == 'nemeth12'):
-            singular_vals = EigenvalueMethods.get_sing_vals("nemeth12_SVD.mat")
-            return np.max(singular_vals)
-        
-        n = len(A)
-        x = [1] * n
-        lambda_old = 0
-        
-        for _ in range(max_iter):
-            x = LinAlgType.matrix_vector_multiply(A, x)
-            lambda_new = LinAlgType.vector_norm(x)
-            x = LinAlgType.vector_scalar_divide(x, lambda_new)
-            if abs(lambda_new - lambda_old) < tol:
-                break
-            lambda_old = lambda_new
-            
-        return lambda_new
-
-    @staticmethod
-    def inverse_power_method(LinAlgType, A, type, max_iter=100, tol=1e-6):
-
-        if (type == 'nemeth12'):
-            singular_vals = EigenvalueMethods.get_sing_vals("nemeth12_SVD.mat")
-            return np.min(singular_vals)
-        
-        n = len(A)
-        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, type, max_iter, tol)

code/linear_algebra_utils.py

@@ -6,6 +6,7 @@ from abc import ABC, abstractmethod
 from concurrent.futures import ThreadPoolExecutor
 from functools import partial
 from time_measurement import time_measurement, time_accumulator
+import numpy as np
 import dask.array as da
 
 class LinearAlgebraUtils(ABC):
@@ -59,11 +60,6 @@ class LinearAlgebraUtils(ABC):
     def matrix_matrix_subtraction(A, B):
         pass
 
-    @staticmethod
-    @abstractmethod
-    def gaussian_elimination(A, b):
-        pass
-
 
 class SequentialLinearAlgebraUtils(ABC):
     @staticmethod
@@ -84,7 +80,7 @@ class SequentialLinearAlgebraUtils(ABC):
 
     @staticmethod
     def matrix_scalar_multiply(A, w):
-        return [[w * A[i][j] for j in range(len(A[0]))] for i in range(len(A))]
+        return A * w
 
     @staticmethod
     def vector_vector_subtraction(v1, v2):
@@ -106,39 +102,12 @@ class SequentialLinearAlgebraUtils(ABC):
     def matrix_matrix_subtraction(A, B):
         return [[A[i][j] - B[i][j] for j in range(len(A[0]))] for i in range(len(A))]
 
-    @staticmethod
-    def gaussian_elimination(A, b):
-        n = len(A)
-        M = [row[:] for row in A]
-        
-        for i in range(n):
-            M[i].append(b[i])
-        
-        for k in range(n):
-            if M[k][k] == 0:
-                for i in range(k + 1, n):
-                    if M[i][k] != 0:
-                        M[k], M[i] = M[i], M[k]
-                        break
-            
-            for i in range(k + 1, n):
-                factor = M[i][k] / M[k][k]
-                for j in range(k, n + 1):
-                    M[i][j] -= factor * M[k][j]
-        
-        x = [0] * n
-        for i in range(n - 1, -1, -1):
-            x[i] = M[i][-1] / M[i][i]
-            for k in range(i - 1, -1, -1):
-                M[k][-1] -= M[k][i] * x[i]
-        
-        return x
-    
 
 class ThreadsLinearAlgebraUtils(ABC):
     NUM_THREADS = 4
 
     @staticmethod
+    @time_measurement(time_accumulator)
     def get_chunk_size(data):
         num_elements = len(data)
         num_threads = min(ThreadsLinearAlgebraUtils.NUM_THREADS, num_elements)
@@ -148,6 +117,7 @@ class ThreadsLinearAlgebraUtils(ABC):
 
 
     @staticmethod
+    @time_measurement(time_accumulator)
     def divide_vectors_to_chunks(v1, v2):
         chunk_size, num_threads, remainder = ThreadsLinearAlgebraUtils.get_chunk_size(v1)
 
@@ -161,6 +131,7 @@ class ThreadsLinearAlgebraUtils(ABC):
         return chunks
     
     @staticmethod
+    @time_measurement(time_accumulator)
     def divide_vector_or_matrix_to_chunks(v):
         chunk_size, num_threads, remainder = ThreadsLinearAlgebraUtils.get_chunk_size(v)
 
@@ -264,8 +235,23 @@ class ThreadsLinearAlgebraUtils(ABC):
     @staticmethod
     @time_measurement(time_accumulator)
     def divide_matrixes_to_chunks(A, B):
-        chunk_size = len(A) // ThreadsLinearAlgebraUtils.NUM_THREADS
+        num_rows = len(A)
-        return [(A[i:i + chunk_size], B[i:i + chunk_size]) for i in range(0, len(A), chunk_size)]
+        num_threads = ThreadsLinearAlgebraUtils.NUM_THREADS
+        if num_threads > num_rows:
+            num_threads = num_rows
+        if num_rows == 0:
+            return []
+        chunk_size = num_rows // num_threads
+        remainder = num_rows % num_threads
+        chunks = []
+        start = 0
+        for _ in range(num_threads):
+            end = start + chunk_size + (1 if remainder > 0 else 0)
+            chunks.append((A[start:end], B[start:end]))
+            start = end
+            if remainder > 0:
+                remainder -= 1
+        return chunks
 
     @staticmethod
     @time_measurement(time_accumulator)
@@ -280,40 +266,6 @@ class ThreadsLinearAlgebraUtils(ABC):
             results = executor.map(subtract_chunk, chunks)
         return [row for chunk in results for row in chunk]
 
-    @staticmethod
-    @time_measurement(time_accumulator)
-    def gaussian_elimination(A, b):
-        n = len(A)
-        M = [row[:] for row in A]
-
-        for i in range(n):
-            M[i].append(b[i])
-
-        for k in range(n):
-            # Pivoting
-            if M[k][k] == 0:
-                for i in range(k + 1, n):
-                    if M[i][k] != 0:
-                        M[k], M[i] = M[i], M[k]
-                        break
-
-            # Threads
-            def eliminate_row(i):
-                factor = M[i][k] / M[k][k]
-                for j in range(k, n + 1):
-                    M[i][j] -= factor * M[k][j]
-
-            with ThreadPoolExecutor(max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS) as executor:
-                rows_to_eliminate = range(k + 1, n)
-                executor.map(eliminate_row, rows_to_eliminate)
-
-        x = [0] * n
-        for i in range(n - 1, -1, -1):
-            x[i] = M[i][-1] / M[i][i]
-            for k in range(i - 1, -1, -1):
-                M[k][-1] -= M[k][i] * x[i]
-
-        return x
 
 @time_measurement(time_accumulator)
 def process_row(params):
@@ -331,6 +283,7 @@ def multiply_by_scalar(pair):
     element, scalar = pair
     return element * scalar
 
+
 class ProcessLinearAlgebraUtils:
     @staticmethod
     @time_measurement(time_accumulator)
@@ -422,56 +375,30 @@ class ProcessLinearAlgebraUtils:
         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 sum_of_squares(row):
+        return sum(x ** 2 for x in row)
 
     @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)
+            row_sums = pool.map(ProcessLinearAlgebraUtils.sum_of_squares, A)
         return math.sqrt(sum(row_sums))
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_matrix_subtraction(A, B):
+    def subtract_rows(row_from_A, row_from_B):
-        def subtract_rows(row_pair):
+        return [a - b for a, b in zip(row_from_A, row_from_B)]
-            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):
+    def matrix_matrix_subtraction(A, B):
-        try:
+        with Pool() as pool:
-            n = len(A)
+            result = pool.starmap(ProcessLinearAlgebraUtils.subtract_rows, zip(A, B))
-            A = [list(row) + [b_i] for row, b_i in zip(A, b)] 
+        return result
-
-            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
 
 class DistributedArraysLinearAlgebraUtils(ABC):
     @staticmethod

code/matrix_generator.py

@@ -1,20 +1,12 @@
 import numpy as np
 import scipy.io
+import os
 
 class MatrixGenerator:
     @staticmethod
     def generate_spd_matrix(n: int) -> np.ndarray:
-        """
-        Generates a random symmetric positive definite matrix of size n x n.
-
-        Parameters:
-        n (int): The size of the matrix to generate.
-
-        Returns:
-        np.ndarray: A symmetric positive definite matrix of size n x n.
-        """
         A = np.random.rand(n, n)
-        spd_matrix = np.dot(A, A.T) + n * np.eye(n)  # Adding n*I ensures positive definiteness
+        spd_matrix = np.dot(A, A.T) + n * MatrixGenerator.generate_identity_matrix(n)  # Adding n*I ensures positive definiteness
         return spd_matrix
 
     @staticmethod
@@ -41,17 +33,54 @@ class MatrixGenerator:
         if type == 'spd':
             if size is None:
                 raise ValueError("Size must be provided for SPD matrix generation.")
-            matrix = MatrixGenerator.generate_spd_matrix(size)
+            try:
-            vector = np.random.uniform(-1, 1, size)
+                matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file("spd_"+str(size)+".npz")
+            except FileNotFoundError as e:
+                matrix = MatrixGenerator.generate_spd_matrix(size)
+                vector = np.random.uniform(-1, 1, size)
+                lambda_min, lambda_max = MatrixGenerator.calculate_eigenvalues(matrix)
+                MatrixGenerator.save_to_file(matrix, vector, lambda_min, lambda_max, "spd_"+str(size)+".npz")
         elif type == 'nemeth12':
-            matrix = -1 * MatrixGenerator.get_matrix_from_file("nemeth12.mat", 1)
+            try:
-            size = matrix.shape[0]
+                matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file("nemeth12.npz")
-            vector = MatrixGenerator.generate_alternating_vector(size)
+            except FileNotFoundError as e:
+                matrix = -1 * MatrixGenerator.get_matrix_from_file("nemeth12.mat", 1)
+                size = matrix.shape[0]
+                vector = MatrixGenerator.generate_alternating_vector(size)
+                lambda_min, lambda_max = MatrixGenerator.calculate_eigenvalues(matrix)
+                MatrixGenerator.save_to_file(matrix, vector, lambda_min, lambda_max, "nemeth12.npz")         
         elif type == 'poli3':
-            matrix = MatrixGenerator.get_matrix_from_file("poli3.mat", 2)
+            try:
-            size = matrix.shape[0]
+                matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file("poli3.npz")
-            vector = MatrixGenerator.generate_alternating_vector(size)
+            except FileNotFoundError as e:
+                matrix = MatrixGenerator.get_matrix_from_file("poli3.mat", 2)
+                size = matrix.shape[0]
+                vector = MatrixGenerator.generate_alternating_vector(size)
+                lambda_min, lambda_max = MatrixGenerator.calculate_eigenvalues(matrix)
+                MatrixGenerator.save_to_file(matrix, vector, lambda_min, lambda_max, "poli3.npz")
         else:
             raise ValueError("Invalid type specified. Choose 'spd', 'nemeth12', or 'poli3'.")
         
-        return matrix, vector
+        return matrix, vector, lambda_min, lambda_max
+    
+    @staticmethod
+    def calculate_eigenvalues(A):
+        eigenvalues = np.linalg.eigvals(A)
+        lambda_min = np.min(eigenvalues)
+        lambda_max = np.max(eigenvalues)
+        return lambda_min, lambda_max
+    
+    @staticmethod
+    def save_to_file(matrix, vector, lambda_min, lambda_max, file_path):
+        np.savez(file_path, matrix=matrix, vector=vector, lambda_min=lambda_min, lambda_max=lambda_max)
+
+    def load_from_file(file_path):
+        if not os.path.exists(file_path):
+            raise FileNotFoundError(f"The file {file_path} does not exist.")
+        data = np.load(file_path)
+        matrix = data['matrix']
+        vector = data['vector']
+        lambda_min = data['lambda_min']
+        lambda_max = data['lambda_max']
+        return matrix, vector, lambda_min, lambda_max
+

code/richardson_method.py

@@ -1,14 +1,14 @@
 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 time_measurement, time_accumulator, tests_time
 import time
 import gc
+import numpy as np
 
 class RichardsonMethod:
     @time_measurement(time_accumulator)
-    def __init__(self, method: ProcessingType, A, type, b, max_iterations, size: int, x0=None, tol=1e-5):
+    def __init__(self, method: ProcessingType, A, b, lambda_min, lambda_max, max_iterations, size: int, x0=None, tol=1e-5):
         self.LinAlg = self.assign_LinAlgType(method)
         self.A = A
         self.b = b
@@ -16,27 +16,21 @@ class RichardsonMethod:
         self.max_iterations = max_iterations
         self.tol = tol
         # self.I = MatrixGenerator.generate_identity_matrix(size)
-        self.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.LinAlg, self.A, type)
+        self.lambda_min = lambda_min
+        self.lambda_max = lambda_max
         if self.lambda_min < 0:
             raise ValueError("Matrix A is not positive semi-definite.")
         self.omega = RichardsonMethod.calculate_omega(self.lambda_min, self.lambda_max)
         
-    @staticmethod
-    def calculate_eigenvalues(LinAlgType, A, type):
-        eigenvalues = np.linalg.eigvals(A)
-        lambda_min = np.min(eigenvalues)
-        lambda_max = np.max(eigenvalues)
-        return lambda_min, lambda_max
-
     @staticmethod
     def calculate_omega(lambda_min, lambda_max):
         return 2 / (lambda_min + lambda_max)
     
     @staticmethod
     def convergence_norm(LinAlgType, A, omega, I) -> bool:
-        wA = LinAlgType.LinAlg.matrix_scalar_multiply(A, omega)
+        wA = LinAlgType.matrix_scalar_multiply(A, omega)
-        IMinuswA = LinAlgType.LinAlg.matrix_matrix_subtraction(I, wA)
+        IMinuswA = LinAlgType.matrix_matrix_subtraction(I, wA)
-        norm = LinAlgType.LinAlg.matrix_norm(IMinuswA)
+        norm = LinAlgType.matrix_norm(IMinuswA)
         return norm
     
     @staticmethod
@@ -58,8 +52,8 @@ class RichardsonMethod:
         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:
+        # 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", 
+        #    return RichardsonMethod.convergence_norm(self.LinAlg, self.A, self.omega, self.I), "Richardson method for those values will NOT converge", 
 
         for iteration in range(self.max_iterations):
             Ax = self.LinAlg.matrix_vector_multiply(self.A, x)

code/tests.py

@@ -1,59 +1,52 @@
 import pytest
 import numpy as np
-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
     difference = I - w * A
-    
-    # Calculate the Euclidean norm of the difference
     norm = np.linalg.norm(difference)
-    
     return norm
 
 def calculate_eigenvalues(A):
-    # Calculate the eigenvalues of matrix A
     eigenvalues = np.linalg.eigvals(A)
-    
-    # Find the minimum and maximum eigenvalues
     lambda_min = np.min(eigenvalues)
     lambda_max = np.max(eigenvalues)
-    
     return lambda_min, lambda_max
 
-def calcualte_norm_from_matrix_numpy(A, n, max_iterations):
+def calcualte_norm_from_matrix_numpy(A, n):
     lambda_min, lambda_max = calculate_eigenvalues(A)
     omega = 2 / (lambda_min + lambda_max)
     I = np.eye(n)
     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, ProcessingType.PROCESSES])
-@pytest.mark.parametrize("matrix_type", ["spd", "nemeth12"])#, "poli3"])
 @time_measurement(tests_time)
+def solution_lib(A, b):
+    return np.linalg.solve(A, b)
+
+@pytest.mark.parametrize("n", [2, 5, 10, 50, 100, 300, 500, 750, 1000, 5000, 10000])
+@pytest.mark.parametrize("processing_type", [ProcessingType.SEQUENTIAL, ProcessingType.THREADS, ProcessingType.PROCESSES])
+@pytest.mark.parametrize("matrix_type", ["spd", "nemeth12", "poli3"])
 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 = 7e-3
+    tolerance = 8e-3
     max_iterations=100
     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)
+        A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('spd', size=n)
     elif matrix_type == "poli3":
-        A, b = MatrixGenerator.generate_matrix_and_vector('poli3')
+        A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('poli3')
     elif matrix_type == "nemeth12":
-        A, b = MatrixGenerator.generate_matrix_and_vector('nemeth12')
+        A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('nemeth12')
     else:
         raise ValueError("Invalid matrix type specified. Choose 'spd', 'poli3', or 'nemeth12'.")
     
-    richardson_solver = RichardsonMethod(processing_type, A, matrix_type, b, max_iterations, size=A.shape[0], tol=1e-7)
+    richardson_solver = RichardsonMethod(processing_type, A, b, lambda_min, lambda_max, max_iterations, size=A.shape[0], tol=1e-7)
-    # solution_richardson, info_richardson = richardson_solver.solve()
 
     solution_richardson, info_richardson = None, None
     with capsys.disabled():
@@ -63,37 +56,24 @@ def test_richardson_vs_cg(n: int, processing_type: ProcessingType, matrix_type:
     captured = capsys.readouterr()
     print("Captured output:", captured.out)
     
-    solution_cg, info = cg(A, b, atol=0.)
+    solution = solution_lib(A,b)
     
-    if info == 0:  # SciPy CG converged
+    assert_converged(solution_richardson, info_richardson, solution, tolerance, A, n)
-        assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n)
-    else:  # SciPy CG did not converge
-        assert_scipy_not_converged(solution_richardson, info_richardson, A, b)
 
-def assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n):
+def assert_converged(solution_richardson, info_richardson, solution, tolerance, A, n):
     if info_richardson == "Richardson method for those values will NOT converge":
-        print("Richardson did not converge, while SciPy did")
+        numpy_norm = calcualte_norm_from_matrix_numpy(A, n)
-        numpy_norm = calcualte_norm_from_matrix_numpy(A, n, max_iterations)
         print("Numpy norm: ", numpy_norm, " Richardson norm: ", solution_richardson)
-        assert False, "Richardson did not converge, while SciPy did"
+        assert False, "Richardson did not converge"
     else:
-        difference = np.linalg.norm(solution_richardson - solution_cg)
+        difference = np.linalg.norm(solution_richardson - solution)
-        print(f"Difference between Richardson and CG solutions: {difference:.8f}")
+        print(f"Difference between Richardson and numpy solutions: {difference:.8f}")
         if difference < tolerance:
-            print("Both Richardson and CG converged and calculated correct values.")
+            print("Both Richardson and numpy method converged and calculated correct values.")
         else:
-            print("Solution CG:\n", solution_cg)
+            print("Solution numpy:\n", solution)
             print("Solution Richardson:\n", solution_richardson)
         assert difference < tolerance, f"The solutions are different! Difference: {difference:.8f}"
 
-def assert_scipy_not_converged(solution_richardson, info_richardson, A, b):
-    if info_richardson == "Richardson method for those values will NOT converge":
-        print("Richardson and SciPy did not converge")
-    else:
-        print("Richardson converged while SciPy did not:", solution_richardson)
-        print("Matrix A:\n", A)
-        print("Vector b:\n", b)
-        assert False, "Richardson converged while SciPy did not"
-        
 if __name__ == "__main__":
     pytest.main()