diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py index 78b723d1..49f2da8e 100644 --- a/code/eigenvalue_methods.py +++ b/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) diff --git a/code/linear_algebra_utils.py b/code/linear_algebra_utils.py index d16b8436..b54ba9a1 100644 --- a/code/linear_algebra_utils.py +++ b/code/linear_algebra_utils.py @@ -5,20 +5,7 @@ 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 - -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]))] - -def divide_by_scalar(pair): - xi, scalar = pair - return xi / scalar - -def multiply_by_scalar(pair): - element, scalar = pair - return element * scalar +from time_measurement import time_measurement, time_accumulator class LinearAlgebraUtils(ABC): @staticmethod @@ -187,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: @@ -195,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: @@ -204,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) @@ -217,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) @@ -226,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: @@ -234,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: @@ -243,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: @@ -251,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: @@ -260,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) @@ -274,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): @@ -293,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] @@ -327,91 +314,123 @@ class ThreadsLinearAlgebraUtils(ABC): 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 - def scalar_matrix_multiply(omega, vector): + @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)] @@ -421,6 +440,7 @@ class ProcessLinearAlgebraUtils: return result @staticmethod + @time_measurement(time_accumulator) def gaussian_elimination(A, b): try: n = len(A) diff --git a/code/matrix_generator.py b/code/matrix_generator.py index 7451062a..a6107f39 100644 --- a/code/matrix_generator.py +++ b/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 + diff --git a/code/nemeth12.mat b/code/nemeth12.mat new file mode 100644 index 00000000..59b31b2a Binary files /dev/null and b/code/nemeth12.mat differ diff --git a/code/poli3.mat b/code/poli3.mat new file mode 100644 index 00000000..24e4abc2 Binary files /dev/null and b/code/poli3.mat differ diff --git a/code/richardson_method.py b/code/richardson_method.py index 591a80c1..10876990 100644 --- a/code/richardson_method.py +++ b/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.") @@ -52,8 +51,8 @@ class RichardsonMethod: 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", @@ -62,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 diff --git a/code/tests.py b/code/tests.py index f9d861e4..97b5b0d5 100644 --- a/code/tests.py +++ b/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) diff --git a/code/time_measurement.py b/code/time_measurement.py index 899ace07..03ca4c0f 100644 --- a/code/time_measurement.py +++ b/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