diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py deleted file mode 100644 index 40639559..00000000 --- a/code/eigenvalue_methods.py +++ /dev/null @@ -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) diff --git a/code/linear_algebra_utils.py b/code/linear_algebra_utils.py index b54ba9a1..162c3563 100644 --- a/code/linear_algebra_utils.py +++ b/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 class LinearAlgebraUtils(ABC): @staticmethod @@ -58,11 +59,6 @@ class LinearAlgebraUtils(ABC): def matrix_matrix_subtraction(A, B): pass - @staticmethod - @abstractmethod - def gaussian_elimination(A, b): - pass - class SequentialLinearAlgebraUtils(ABC): @staticmethod @@ -83,7 +79,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): @@ -105,39 +101,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) @@ -147,6 +116,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) @@ -160,6 +130,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) @@ -263,8 +234,23 @@ class ThreadsLinearAlgebraUtils(ABC): @staticmethod @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)] + num_rows = len(A) + 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) @@ -279,40 +265,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): @@ -330,6 +282,7 @@ def multiply_by_scalar(pair): element, scalar = pair return element * scalar + class ProcessLinearAlgebraUtils: @staticmethod @time_measurement(time_accumulator) @@ -421,53 +374,27 @@ 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_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 - + def subtract_rows(row_from_A, row_from_B): + return [a - b for a, b in zip(row_from_A, row_from_B)] + @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 + def matrix_matrix_subtraction(A, B): + with Pool() as pool: + result = pool.starmap(ProcessLinearAlgebraUtils.subtract_rows, zip(A, B)) + return result diff --git a/code/matrix_generator.py b/code/matrix_generator.py index 140f1b11..024403a3 100644 --- a/code/matrix_generator.py +++ b/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) - vector = np.random.uniform(-1, 1, size) + try: + 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) - size = matrix.shape[0] - vector = MatrixGenerator.generate_alternating_vector(size) + try: + matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file("nemeth12.npz") + 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) - size = matrix.shape[0] - vector = MatrixGenerator.generate_alternating_vector(size) + try: + matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file("poli3.npz") + 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 \ No newline at end of file + 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 + diff --git a/code/nemeth12_SVD.mat b/code/nemeth12_SVD.mat deleted file mode 100644 index 9a569d11..00000000 Binary files a/code/nemeth12_SVD.mat and /dev/null differ diff --git a/code/poli3_SVD.mat b/code/poli3_SVD.mat deleted file mode 100644 index db51d0e7..00000000 Binary files a/code/poli3_SVD.mat and /dev/null differ diff --git a/code/richardson_method.py b/code/richardson_method.py index 089ac8cd..3e51dda5 100644 --- a/code/richardson_method.py +++ b/code/richardson_method.py @@ -1,5 +1,4 @@ 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 @@ -9,7 +8,7 @@ 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 @@ -17,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) - IMinuswA = LinAlgType.LinAlg.matrix_matrix_subtraction(I, wA) - norm = LinAlgType.LinAlg.matrix_norm(IMinuswA) + wA = LinAlgType.matrix_scalar_multiply(A, omega) + IMinuswA = LinAlgType.matrix_matrix_subtraction(I, wA) + norm = LinAlgType.matrix_norm(IMinuswA) return norm @staticmethod @@ -58,8 +51,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: - # return RichardsonMethod.convergence_norm(self.A, self.omega, self.I), "Richardson method for those values will NOT converge", + # if RichardsonMethod.convergence_norm(self.LinAlg, self.A, self.omega, self.I) >= 1: + # 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) diff --git a/code/tests.py b/code/tests.py index c112f009..45b8f508 100644 --- a/code/tests.py +++ b/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) - # solution_richardson, info_richardson = richardson_solver.solve() + richardson_solver = RichardsonMethod(processing_type, A, b, lambda_min, lambda_max, max_iterations, size=A.shape[0], tol=1e-7) 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_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) + assert_converged(solution_richardson, info_richardson, solution, tolerance, A, n) -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, max_iterations) + numpy_norm = calcualte_norm_from_matrix_numpy(A, n) 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) - print(f"Difference between Richardson and CG solutions: {difference:.8f}") + difference = np.linalg.norm(solution_richardson - solution) + 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()