From b02117c11b22a09992b4ed53a819259459f1a44d Mon Sep 17 00:00:00 2001 From: Gromiusz <104679774+Gromiusz@users.noreply.github.com> Date: Fri, 25 Oct 2024 17:01:19 +0200 Subject: [PATCH] Add Interface for Extending Parallel Processing Methods (#4) * Add choice between sequential or parallel(threads) processing The LinearAlgebraUtils class has been changed to abstract class. * remove needless imports and move the functionality from constructor to the function * add ProcessingType class for a more transparent selection of methods * the refactor of the assign_method * rename the assign_method method --- code/eigenvalue_methods.py | 16 +++--- code/linear_algebra_utils.py | 108 ++++++++++++++++++++++++++++++++++- code/processing_type.py | 5 ++ code/richardson_method.py | 42 +++++++++----- code/tests.py | 3 +- 5 files changed, 149 insertions(+), 25 deletions(-) create mode 100644 code/processing_type.py diff --git a/code/eigenvalue_methods.py b/code/eigenvalue_methods.py index 0d97defa..78b723d1 100644 --- a/code/eigenvalue_methods.py +++ b/code/eigenvalue_methods.py @@ -1,16 +1,14 @@ -from linear_algebra_utils import LinearAlgebraUtils - class EigenvalueMethods: @staticmethod - def power_method(A, max_iter, tol=1e-6): + def power_method(LinAlgType, A, max_iter, tol=1e-6): n = len(A) x = [1] * n lambda_old = 0 for _ in range(max_iter): - x = LinearAlgebraUtils.matrix_vector_multiply(A, x) - lambda_new = LinearAlgebraUtils.vector_norm(x) - x = LinearAlgebraUtils.vector_scalar_divide(x, lambda_new) + 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 @@ -18,10 +16,10 @@ class EigenvalueMethods: return lambda_new @staticmethod - def inverse_power_method(A, max_iter, tol=1e-6): + def inverse_power_method(LinAlgType, A, max_iter, tol=1e-6): n = len(A) I = [[1 if i == j else 0 for j in range(n)] for i in range(n)] - A_inv = [LinearAlgebraUtils.gaussian_elimination(A.tolist(), I_col) for I_col in I] + 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(A_inv, max_iter, tol) + 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 5ff15f83..a159f42f 100644 --- a/code/linear_algebra_utils.py +++ b/code/linear_algebra_utils.py @@ -1,13 +1,71 @@ import math +from abc import ABC, abstractmethod -class LinearAlgebraUtils: +class LinearAlgebraUtils(ABC): + @staticmethod + @abstractmethod + def dot_product(v1, v2): + pass + + @staticmethod + @abstractmethod + def matrix_vector_multiply(A, x): + pass + + @staticmethod + @abstractmethod + def vector_norm(v): + pass + + @staticmethod + @abstractmethod + def vector_scalar_divide(x, scalar): + pass + + @staticmethod + @abstractmethod + def matrix_scalar_multiply(A, w): + pass + + @staticmethod + @abstractmethod + def vector_vector_subtraction(v1, v2): + pass + + @staticmethod + @abstractmethod + def vector_vector_addition(v1, v2): + pass + + @staticmethod + @abstractmethod + def scalar_matrix_multiply(omega, vector): + pass + + @staticmethod + @abstractmethod + def matrix_norm(A): + pass + + @staticmethod + @abstractmethod + def matrix_matrix_subtraction(A, B): + pass + + @staticmethod + @abstractmethod + def gaussian_elimination(A, b): + pass + + +class SequentialLinearAlgebraUtils(ABC): @staticmethod def dot_product(v1, v2): return sum(x*y for x, y in zip(v1, v2)) @staticmethod def matrix_vector_multiply(A, x): - return [LinearAlgebraUtils.dot_product(row, x) for row in A] + return [SequentialLinearAlgebraUtils.dot_product(row, x) for row in A] @staticmethod def vector_norm(v): @@ -69,3 +127,49 @@ class LinearAlgebraUtils: M[k][-1] -= M[k][i] * x[i] return x + + +class ThreadsLinearAlgebraUtils(ABC): + @staticmethod + def dot_product(v1, v2): + pass + + @staticmethod + def matrix_vector_multiply(A, x): + pass + + @staticmethod + def vector_norm(v): + pass + + @staticmethod + def vector_scalar_divide(x, scalar): + pass + + @staticmethod + def matrix_scalar_multiply(A, w): + pass + + @staticmethod + def vector_vector_subtraction(v1, v2): + pass + + @staticmethod + def vector_vector_addition(v1, v2): + pass + + @staticmethod + def scalar_matrix_multiply(omega, vector): + pass + + @staticmethod + def matrix_norm(A): + pass + + @staticmethod + def matrix_matrix_subtraction(A, B): + pass + + @staticmethod + def gaussian_elimination(A, b): + pass diff --git a/code/processing_type.py b/code/processing_type.py new file mode 100644 index 00000000..daad3e54 --- /dev/null +++ b/code/processing_type.py @@ -0,0 +1,5 @@ +from enum import Enum, auto + +class ProcessingType(Enum): + SEQUENTIAL = auto() + THREADS = auto() \ No newline at end of file diff --git a/code/richardson_method.py b/code/richardson_method.py index 47b08394..8fb02275 100644 --- a/code/richardson_method.py +++ b/code/richardson_method.py @@ -1,43 +1,59 @@ from linear_algebra_utils import LinearAlgebraUtils +from linear_algebra_utils import SequentialLinearAlgebraUtils +from linear_algebra_utils import ThreadsLinearAlgebraUtils from eigenvalue_methods import EigenvalueMethods from matrix_generator import MatrixGenerator +from processing_type import ProcessingType class RichardsonMethod: - def __init__(self, A, b, max_iterations, size: int, x0=None, tol=1e-5): + 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 self.b = b 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.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.A, max_iterations) + 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.") self.omega = RichardsonMethod.calculate_omega(self.lambda_min, self.lambda_max) - + @staticmethod - def calculate_eigenvalues(A, max_iterations): - return EigenvalueMethods.inverse_power_method(A, max_iterations), EigenvalueMethods.power_method(A, max_iterations) + def calculate_eigenvalues(LinAlgType, A, max_iterations): + return EigenvalueMethods.inverse_power_method(LinAlgType, A, max_iterations), EigenvalueMethods.power_method(LinAlgType, A, max_iterations) @staticmethod def calculate_omega(lambda_min, lambda_max): return 2 / (lambda_min + lambda_max) @staticmethod - def convergence_norm(A, omega, I) -> bool: - wA = LinearAlgebraUtils.matrix_scalar_multiply(A, omega) - IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(I, wA) - norm = LinearAlgebraUtils.matrix_norm(IMinuswA) + 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) return norm + + @staticmethod + def assign_LinAlgType(method): + metody = { + ProcessingType.SEQUENTIAL: SequentialLinearAlgebraUtils, + ProcessingType.THREADS: ThreadsLinearAlgebraUtils + } + + try: + return metody[method] + except KeyError: + raise ValueError("Unknown method, please use 'SEQUENTIAL' or 'THREADS'.") def solve(self): x = self.x0[:] - #if RichardsonMethod.convergence_norm(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", for iteration in range(self.max_iterations): - Ax = LinearAlgebraUtils.matrix_vector_multiply(self.A, x) - residual = LinearAlgebraUtils.vector_vector_subtraction(self.b, Ax) - x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(self.omega, residual)) + 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_matrix_multiply(self.omega, residual)) return x, 0 diff --git a/code/tests.py b/code/tests.py index 1bdf7961..675916c0 100644 --- a/code/tests.py +++ b/code/tests.py @@ -3,6 +3,7 @@ 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 def calculate_norm_numpy(I, w, A): # Calculate the difference between I and w * A @@ -37,7 +38,7 @@ def test_richardson_vs_cg(n: int): tolerance = 1e-5 max_iterations=1000 A, b = MatrixGenerator.generate_random_matrix_and_vector(n) - richardson_solver = RichardsonMethod(A, b, max_iterations, size=n, tol=1e-7) + richardson_solver = RichardsonMethod(ProcessingType.SEQUENTIAL , A, b, max_iterations, size=n, tol=1e-7) solution_richardson, info_richardson = richardson_solver.solve() solution_cg, info = cg(A, b)