From 219570e4f10813069dc2c83eb39a35aaa4c08b2a Mon Sep 17 00:00:00 2001 From: Niklas Birk Date: Sat, 27 Apr 2024 18:21:33 +0200 Subject: [PATCH] Refactoring and more parallel stuff --- src/cg.py | 73 +++++++++++++------------------------- src/main.py | 21 +++++++---- src/matrix.py | 23 +++++++----- src/matrix_mpi.py | 90 ++++++++++++++++++++++++++++++++--------------- src/vector.py | 13 +++++++ src/vector_mpi.py | 80 +++++++++++++++++------------------------ 6 files changed, 162 insertions(+), 138 deletions(-) diff --git a/src/cg.py b/src/cg.py index e60b501..3302519 100644 --- a/src/cg.py +++ b/src/cg.py @@ -1,58 +1,35 @@ -from mpi4py import MPI - from matrix_mpi import MatrixMPI as Matrix from vector_mpi import VectorMPI as Vector -comm = MPI.COMM_WORLD -size = comm.Get_size() -rank = comm.Get_rank() +# from matrix import Matrix +# from vector import Vector -def cg(n: int, A: Matrix, f: Vector, tol: float): - # Intialisierung des Startvektors x - x = Vector([1] * n) +def cg(A: Matrix, x0: Vector, b: Vector, tolerance: float = 1e-3, max_iterations: int = 1_000): + """ + Solves a system of linear equations of the form Ax = b numerically. + + :param A: The transformation matrix A + :param x0: A vector to start the algorithm with + :param b: The solution vector of the system of linear equations, the right hand side + :param tolerance: The tolerance at which to stop the algorithm, default is 0.001 + :param max_iterations: Maximum number of iterations, default is 1000 + """ + iterations = 0 - # Anzahl der Schritte - count = 0 + x = x0 + r = b - A * x + d = r - # Anfangswerte berechnen - r = f - A * x # Anfangsresiduum - p = r # Anfangsabstiegsrichtung + while r.norm() >= tolerance and iterations < max_iterations: + z = A * d - while r.norm() > tol and count < 1000: - print(f"{count}. Iterationsschritt:\n") - # print("Iterierte:", x) - # print("Residuumsnorm: ", r.norm()) + alpha = (r.T() * d) / (d.T() * z) + x = x + alpha * d + r = r - alpha * z - z = A * p # Matrix-Vektorprodukt berechnen und speichern + beta = -(r.T() * z) / (d.T() * z) + d = r + beta * d - # Minimiere phi in Richung p um neue Iterierte x zu finden - alpha = (r.T() * p) / (p.T() * z) # (np.dot(r , p)) / (np.dot(p , z)) - # print(alpha) - - x = x + alpha * p # neue Itterierte x - r = r - alpha * z # neues Residuum - - # Bestimmung der neuen Suchrichtung - beta = - (r.T() * z) / (p.T() * z) # (np.dot(r , z)) / (np.dot(p , z)) - p = r + beta * p # neue konjugierte Abstiegsrichtung - - count = count + 1 - - print(f"{rank} APFELSTRUDEL") - # if rank == 0: - # # Vergleich mit numpy-interner Lsg - # u = np.linalg.solve(np.array(A.get_data()), np.array(f.get_data())) - # - # print("Lösung mit CG-Verfahren:", x) - # print("Numpy interne Lösung:", u) - # - # if (Vector(u) - x).norm() > eps: - # print("Der CG-Algorithmus hat nicht richtig funktioniert!") - # else: - # print("Der CG-Algorithmus war erfolgreich.") - # - # plt.plot(x.get_data(), linewidth=2) - # plt.plot(u, linewidth=2) - # - # plt.show() + iterations = iterations + 1 + return x diff --git a/src/main.py b/src/main.py index bab7dd2..05314ae 100644 --- a/src/main.py +++ b/src/main.py @@ -1,18 +1,25 @@ -import numpy as np +from mpi4py import MPI import cg from matrix_mpi import MatrixMPI as Matrix from vector_mpi import VectorMPI as Vector +# from matrix import Matrix +# from vector import Vector + +comm = MPI.COMM_WORLD +size = comm.Get_size() +rank = comm.Get_rank() + n = 1_00 h = 1 / (n - 1) -# Initialisierung der Matrix A und des Vektor f für LGS Au = f -A = Matrix(np.diag(-1 * np.ones(n - 1), k=1) + np.diag(2 * np.ones(n), k=0) + np.diag(-1 * np.ones(n - 1), k=-1)) -f = Vector([h ** 2 * 2] * n) +A = Matrix([-1, 2, -1], structure="tridiagonal", n=n) +x0 = Vector([1] * n) +b = Vector([h**2 * 2] * n) -# Toleranz epsilon -tol = 0.001 +x = cg.cg(A, x0, b) -cg.cg(n, A, f, tol) +if rank == 0: + print(x) diff --git a/src/matrix.py b/src/matrix.py index 06f0eb4..a5e776c 100644 --- a/src/matrix.py +++ b/src/matrix.py @@ -22,7 +22,6 @@ class Matrix: - ``Matrix(list, str, int)``: will create a new square matrix of given size and structure of either \"unity\", \"diagonal\" or \"tridiagonal\" - ``Matrix(str, int)``: will create a new square matrix of given size and TODO - :param data: Either a list or an numpy ndarray :param shape: A tuple containing the amount of rows and columns :param structure: Either \"unity\", \"diagonal\" or \"tridiagonal\" @@ -84,6 +83,16 @@ class Matrix: """ return self.__data__ + @staticmethod + def flatten_internal(matrices): + flattened_data = [] + rows = 0 + for matrix in matrices: + flattened_data.extend(matrix.get_data()) + rows += matrix.__shape__[0] + cols = matrices[0].__shape__[1] + return flattened_data, (rows, cols) + @staticmethod def flatten(matrices: list): """ @@ -95,13 +104,8 @@ class Matrix: :return: A ``Matrix`` extended by all matrices in the list. :rtype: ``Matrix`` """ - flattened_data = [] - rows = 0 - for matrix in matrices: - flattened_data.extend(matrix.get_matrix()) - rows += matrix.__shape__[0] - cols = matrices[0].__shape__[1] - return Matrix(flattened_data, (rows, cols)) + flattened_data, shape = Matrix.flatten_internal(matrices) + return Matrix(flattened_data, shape) def shape(self): """ @@ -247,6 +251,9 @@ class Matrix: def __rmul__(self, other): return self * other + def get_abs_sum_of_squares(self): + return self.__abs_sum_of_squares__() + def __abs_sum_of_squares__(self): rows = self.__shape__[0] cols = self.__shape__[1] diff --git a/src/matrix_mpi.py b/src/matrix_mpi.py index 98c436b..c23f53f 100644 --- a/src/matrix_mpi.py +++ b/src/matrix_mpi.py @@ -1,3 +1,5 @@ +import math + import numpy from mpi4py import MPI @@ -9,39 +11,75 @@ class MatrixMPI: __mpi_size__ = __mpi_comm__.Get_size() __mpi_rank__ = __mpi_comm__.Get_rank() - __data__: Matrix = None + __data__ = None + __rank_subdata__ = None __chunk__: list = None def __init__(self, data=None, shape=None, structure=None, model=None, offset=None, n=None): - self.__data__ = Matrix(data=data, shape=shape, structure=structure, model=model, offset=offset, n=n) + """ + Creates a new matrix. + The type of the matrix depends on the signature and arguments. + - ``MatrixMPI(list)``: will create a new matrix with the given data in the list and its shape. + - ``MatrixMPI(numpy.ndarray)``: will create a new matrix with the given data in ndarray and its shape. + - ``MatrixMPI(list, (int,int))``: will create a new nxm matrix with the given rows and columns and data in list. + - ``MatrixMPI(list, str, int, int)``: will create a new square matrix of given size and structure of \"diagonal\" + - ``MatrixMPI(list, str, int)``: will create a new square matrix of given size and structure of either \"unity\", \"diagonal\" or \"tridiagonal\" + - ``MatrixMPI(str, int)``: will create a new square matrix of given size and TODO + + :param data: Either a list or an numpy ndarray + :param shape: A tuple containing the amount of rows and columns + :param structure: Either \"unity\", \"diagonal\" or \"tridiagonal\" + :param model: TODO + :param offset: Offset to diagonal axis + :param n: Amount of rows of a square matrix or offset in case of diagonal structure + + :type data: Matrix | list | numpy.ndarray + :type shape: (int, int) + :type structure: str + :type model: str + :type offset: int + :type n: int + + :rtype: MatrixMPI + """ + if isinstance(data, Matrix): + self.__data__ = data + else: + self.__data__ = Matrix(data=data, shape=shape, structure=structure, model=model, offset=offset, n=n) + + # Calculate how much rows are delegated to the rank total_amount_of_rows = self.__data__.shape()[0] chunks = numpy.array_split(list(range(total_amount_of_rows)), self.__mpi_size__) self.__chunk__ = chunks[self.__mpi_rank__].tolist() + # Store the delegated rows explicitly for calculations + rows = len(self.__chunk__) + cols = self.__data__.shape()[1] + self.__rank_subdata__ = Matrix(self.__data__[self.__chunk__], (rows, cols)) + @staticmethod def of(matrix: Matrix): - return MatrixMPI(matrix.get_data(), matrix.shape()) - - def __str__(self): - return str(self.__data__) + return MatrixMPI(matrix) def shape(self): return self.__data__.shape() - def get_rank_submatrix(self): - rows = len(self.__chunk__) - cols = self.__data__.shape()[1] - return Matrix(self.__data__[self.__chunk__], (rows, cols)) - - def get_matrix(self): + def get_rank_subdata(self): """ - Returns the ``Matrix`` that is used internally + Returns only the delegated rows of the rank as ``Matrix`` + :return: The delegated rows as ``Matrix`` + """ + return self.__rank_subdata__ + + def get_data(self): + """ + Returns the whole ``Matrix`` that is used internally :return: The ``Matrix`` that is used internally """ return self.__data__ - def get_data(self): + def get_internal_data(self): """ Returns the raw data of the internal data structure :return: The raw data of the internal data structure @@ -75,16 +113,12 @@ class MatrixMPI: return self.__data__ == other def __neg__(self): - gathered_data = self.__mpi_comm__.gather(-self.get_rank_submatrix()) - data = self.__mpi_comm__.bcast(gathered_data) - return MatrixMPI.of(Matrix.flatten(data)) + return MatrixMPI.of(Matrix.flatten(self.__mpi_comm__.allgather(-self.__rank_subdata__))) def __add__(self, other): if isinstance(other, MatrixMPI): - other = other.get_rank_submatrix() - gathered_data = self.__mpi_comm__.gather(self.get_rank_submatrix() + other) - data = self.__mpi_comm__.bcast(gathered_data) - return MatrixMPI.of(Matrix.flatten(data)) + other = other.__rank_subdata__ + return MatrixMPI.of(Matrix.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ + other))) def __radd__(self, other): return self + other @@ -96,16 +130,12 @@ class MatrixMPI: return -self + other def __truediv__(self, other): - gathered_data = self.__mpi_comm__.gather(self.get_rank_submatrix() / other) - data = self.__mpi_comm__.bcast(gathered_data) - return MatrixMPI.of(Matrix.flatten(data)) + return MatrixMPI.of(Matrix.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ / other))) def __mul__(self, other): if isinstance(other, MatrixMPI): - other = other.get_matrix() - gathered_data = self.__mpi_comm__.gather(self.get_rank_submatrix() * other) - data = self.__mpi_comm__.bcast(gathered_data) - return MatrixMPI.of(Matrix.flatten(data)) + other = other.get_data() + return MatrixMPI.of(Matrix.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ * other))) def __rmul__(self, other): return self * other @@ -121,6 +151,10 @@ class MatrixMPI: :return: the norm as a number """ + if f == "frobenius": + return math.sqrt(self.__mpi_comm__.allreduce(self.__rank_subdata__.get_abs_sum_of_squares())) + elif f == "row sum": + return max(self.__mpi_comm__.allgather(self.__rank_subdata__.norm(f))) return self.__data__.norm(f) def __getitem__(self, key): diff --git a/src/vector.py b/src/vector.py index 7f0ea93..ad7455a 100644 --- a/src/vector.py +++ b/src/vector.py @@ -24,6 +24,19 @@ class Vector(Matrix): else: raise ValueError("data must be a ``list``, a ``numpy.ndarray`` or an integer for dimension") + @staticmethod + def flatten(vectors: list): + """ + Flattens a list of matrices into one bigger matrix. + The columns must match the first ``Matrix`` in the list and the rows can be arbitrarily. + + :param vectors: A list of vectors. + :type vectors: list + :return: A ``Vector`` extended by all matrices in the list. + """ + flattened_data, shape = Matrix.flatten_internal(vectors) + return Vector(flattened_data, shape) + def __eq__(self, other): """ Return ``self==value`` diff --git a/src/vector_mpi.py b/src/vector_mpi.py index b609288..35341e0 100644 --- a/src/vector_mpi.py +++ b/src/vector_mpi.py @@ -1,45 +1,30 @@ +import math + +import numpy +from mpi4py import MPI + from matrix_mpi import MatrixMPI from vector import Vector class VectorMPI(MatrixMPI): - __data__: Vector = None - def __init__(self, data=None, shape=None): - self.__data__ = Vector(data=data, shape=shape) + if isinstance(data, Vector): + self.__data__ = data + else: + self.__data__ = Vector(data=data, shape=shape) + + # Calculate how much rows are delegated to the rank + total_amount_of_rows = self.__data__.shape()[0] + chunks = numpy.array_split(list(range(total_amount_of_rows)), self.__mpi_size__) + self.__chunk__ = chunks[self.__mpi_rank__].tolist() + + # Store the delegated rows explicitly for calculations + self.__rank_subdata__ = Vector(self.__data__[self.__chunk__]) @staticmethod def of(vector: Vector): - return VectorMPI(vector.get_data(), vector.shape()) - - def get_vector(self): - """ - Returns the ``Vector`` that is used internally - :return: The ``Vector`` that is used internally - """ - return self.__data__ - - def get_data(self): - """ - Returns the raw data of the internal data structure - :return: The raw data of the internal data structure - """ - return self.__data__.get_data() - - def shape(self): - return self.__data__.shape() - - def __eq__(self, other): - """ - Return ``self==value`` - - :param other: The object to compare to; must be either a ``Vector``, a ``list`` or a ``numpy.ndarray`` - :return: True if data in the same-shaped vectors are equal to the given data in other for each component otherwise False - """ - if isinstance(other, VectorMPI): - return self.__data__ == other.__data__ - else: - return self.__data__ == other + return VectorMPI(vector) def transpose(self): """ @@ -50,29 +35,36 @@ class VectorMPI(MatrixMPI): def T(self): return self.transpose() + def __str__(self): + return str(self.__data__) + def __neg__(self): return VectorMPI.of(-self.__data__) def __add__(self, other): if isinstance(other, VectorMPI): - other = other.__data__ - return VectorMPI.of(self.__data__ + other) + other = other.__rank_subdata__ + return VectorMPI.of(Vector.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ + other))) def __mul__(self, other): if isinstance(other, VectorMPI): other = other.__data__ - result = self.__data__ * other + + if isinstance(other, int) or isinstance(other, float): + result = Vector.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ * other)) + else: + result = self.__data__ * other return VectorMPI.of(result) if isinstance(result, Vector) else result def __rmul__(self, other): if isinstance(other, MatrixMPI): - return VectorMPI.of(other.get_matrix() * self.get_vector()) + return VectorMPI.of(Vector.flatten(self.__mpi_comm__.allgather(other.get_rank_subdata() * self.get_data()))) return self * other def __truediv__(self, other): if isinstance(other, VectorMPI): - other = other.__data__ - return VectorMPI.of(self.__data__ / other) + other = other.__rank_subdata__ + return VectorMPI.of(Vector.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ / other))) def norm(self, **kwargs): """ @@ -81,7 +73,7 @@ class VectorMPI(MatrixMPI): :param kwargs: ignored :return: the 2-norm of the vector """ - return self.__data__.norm() + return math.sqrt(self.__mpi_comm__.allreduce(self.__rank_subdata__.get_abs_sum_of_squares())) def normalize(self): """ @@ -90,10 +82,4 @@ class VectorMPI(MatrixMPI): :return: the normalized vector """ - return VectorMPI.of(self.__data__ / self.norm()) - - def __getitem__(self, key): - return self.__data__[key] - - def __setitem__(self, key, value): - self.__data__[key] = value + return VectorMPI.of(Vector.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ / self.norm())))