diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py new file mode 100644 index 00000000..40639559 --- /dev/null +++ b/code/eigenvalue_methods.py @@ -0,0 +1,43 @@ +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 162c3563..b54ba9a1 100644 --- a/code/linear_algebra_utils.py +++ b/code/linear_algebra_utils.py @@ -6,7 +6,6 @@ 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 @@ -59,6 +58,11 @@ class LinearAlgebraUtils(ABC): def matrix_matrix_subtraction(A, B): pass + @staticmethod + @abstractmethod + def gaussian_elimination(A, b): + pass + class SequentialLinearAlgebraUtils(ABC): @staticmethod @@ -79,7 +83,7 @@ class SequentialLinearAlgebraUtils(ABC): @staticmethod def matrix_scalar_multiply(A, w): - return A * w + return [[w * A[i][j] for j in range(len(A[0]))] for i in range(len(A))] @staticmethod def vector_vector_subtraction(v1, v2): @@ -101,12 +105,39 @@ 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) @@ -116,7 +147,6 @@ 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) @@ -130,7 +160,6 @@ 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) @@ -234,23 +263,8 @@ class ThreadsLinearAlgebraUtils(ABC): @staticmethod @time_measurement(time_accumulator) def divide_matrixes_to_chunks(A, B): - 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 + 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(time_accumulator) @@ -265,6 +279,40 @@ 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): @@ -282,7 +330,6 @@ def multiply_by_scalar(pair): element, scalar = pair return element * scalar - class ProcessLinearAlgebraUtils: @staticmethod @time_measurement(time_accumulator) @@ -374,27 +421,53 @@ 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(ProcessLinearAlgebraUtils.sum_of_squares, A) + 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 subtract_rows(row_from_A, row_from_B): - return [a - b for a, b in zip(row_from_A, row_from_B)] - + 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 matrix_matrix_subtraction(A, B): - with Pool() as pool: - result = pool.starmap(ProcessLinearAlgebraUtils.subtract_rows, zip(A, B)) - return result + 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 diff --git a/code/matrix_generator.py b/code/matrix_generator.py index 024403a3..140f1b11 100644 --- a/code/matrix_generator.py +++ b/code/matrix_generator.py @@ -1,12 +1,20 @@ 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 * MatrixGenerator.generate_identity_matrix(n) # Adding n*I ensures positive definiteness + spd_matrix = np.dot(A, A.T) + n * np.eye(n) # Adding n*I ensures positive definiteness return spd_matrix @staticmethod @@ -33,54 +41,17 @@ class MatrixGenerator: if type == 'spd': if size is None: raise ValueError("Size must be provided for SPD matrix generation.") - 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") + matrix = MatrixGenerator.generate_spd_matrix(size) + vector = np.random.uniform(-1, 1, size) elif type == 'nemeth12': - 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") + matrix = -1 * MatrixGenerator.get_matrix_from_file("nemeth12.mat", 1) + size = matrix.shape[0] + vector = MatrixGenerator.generate_alternating_vector(size) elif type == 'poli3': - 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") + 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, 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 - + return matrix, vector \ No newline at end of file diff --git a/code/nemeth12_SVD.mat b/code/nemeth12_SVD.mat new file mode 100644 index 00000000..9a569d11 Binary files /dev/null and b/code/nemeth12_SVD.mat differ diff --git a/code/poli3_SVD.mat b/code/poli3_SVD.mat new file mode 100644 index 00000000..db51d0e7 Binary files /dev/null and b/code/poli3_SVD.mat differ diff --git a/code/richardson_method.py b/code/richardson_method.py index 3e51dda5..34e836ea 100644 --- a/code/richardson_method.py +++ b/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, b, lambda_min, lambda_max, max_iterations, size: int, x0=None, tol=1e-5): + def __init__(self, method: ProcessingType, A, type, b, max_iterations, size: int, x0=None, tol=1e-5): self.LinAlg = self.assign_LinAlgType(method) self.A = A self.b = b @@ -16,21 +16,24 @@ class RichardsonMethod: self.max_iterations = max_iterations self.tol = tol # self.I = MatrixGenerator.generate_identity_matrix(size) - self.lambda_min = lambda_min - self.lambda_max = lambda_max + self.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.LinAlg, self.A, type) 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): + return EigenvalueMethods.inverse_power_method(LinAlgType, A, type), EigenvalueMethods.power_method(LinAlgType, A, type) + @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.matrix_scalar_multiply(A, omega) - IMinuswA = LinAlgType.matrix_matrix_subtraction(I, wA) - norm = LinAlgType.matrix_norm(IMinuswA) + wA = LinAlgType.LinAlg.matrix_scalar_multiply(A, omega) + IMinuswA = LinAlgType.LinAlg.matrix_matrix_subtraction(I, wA) + norm = LinAlgType.LinAlg.matrix_norm(IMinuswA) return norm @staticmethod @@ -51,8 +54,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.LinAlg, 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.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 45b8f508..c112f009 100644 --- a/code/tests.py +++ b/code/tests.py @@ -1,52 +1,59 @@ 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): +def calcualte_norm_from_matrix_numpy(A, n, max_iterations): lambda_min, lambda_max = calculate_eigenvalues(A) omega = 2 / (lambda_min + lambda_max) I = np.eye(n) return calculate_norm_numpy(I, omega, A) -@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("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"]) +@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 = 8e-3 + tolerance = 7e-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, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('spd', size=n) + A, b = MatrixGenerator.generate_matrix_and_vector('spd', size=n) elif matrix_type == "poli3": - A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('poli3') + A, b = MatrixGenerator.generate_matrix_and_vector('poli3') elif matrix_type == "nemeth12": - A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('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, lambda_min, lambda_max, max_iterations, size=A.shape[0], tol=1e-7) + 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() solution_richardson, info_richardson = None, None with capsys.disabled(): @@ -56,24 +63,37 @@ def test_richardson_vs_cg(n: int, processing_type: ProcessingType, matrix_type: captured = capsys.readouterr() print("Captured output:", captured.out) - solution = solution_lib(A,b) + solution_cg, info = cg(A, b, atol=0.) - assert_converged(solution_richardson, info_richardson, solution, tolerance, A, n) + 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) -def 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): if info_richardson == "Richardson method for those values will NOT converge": - numpy_norm = calcualte_norm_from_matrix_numpy(A, n) + print("Richardson did not converge, while SciPy did") + 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" + assert False, "Richardson did not converge, while SciPy did" else: - difference = np.linalg.norm(solution_richardson - solution) - print(f"Difference between Richardson and numpy solutions: {difference:.8f}") + difference = np.linalg.norm(solution_richardson - solution_cg) + print(f"Difference between Richardson and CG solutions: {difference:.8f}") if difference < tolerance: - print("Both Richardson and numpy method converged and calculated correct values.") + print("Both Richardson and CG converged and calculated correct values.") else: - print("Solution numpy:\n", solution) + print("Solution CG:\n", solution_cg) 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()