Refactoring and more parallel stuff

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
+# from matrix import Matrix
-size = comm.Get_size()
+# from vector import Vector
-rank = comm.Get_rank()
 
 
-def cg(n: int, A: Matrix, f: Vector, tol: float):
+def cg(A: Matrix, x0: Vector, b: Vector, tolerance: float = 1e-3, max_iterations: int = 1_000):
-    # Intialisierung des Startvektors x
+    """
-    x = Vector([1] * n)
+    Solves a system of linear equations of the form Ax = b numerically.
     
-    # Anzahl der Schritte
+    :param A: The transformation matrix A
-    count = 0
+    :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
 
-    # Anfangswerte berechnen
+    x = x0
-    r = f - A * x  # Anfangsresiduum
+    r = b - A * x
-    p = r  # Anfangsabstiegsrichtung
+    d = r
 
-    while r.norm() > tol and count < 1000:
+    while r.norm() >= tolerance and iterations < max_iterations:
-        print(f"{count}. Iterationsschritt:\n")
+        z = A * d
-        # print("Iterierte:", x)
-        # print("Residuumsnorm: ", r.norm())
 
-        z = A * p  # Matrix-Vektorprodukt berechnen und speichern
+        alpha = (r.T() * d) / (d.T() * z)
+        x = x + alpha * d
+        r = r - alpha * z
 
-        # Minimiere phi in Richung p um neue Iterierte x zu finden
+        beta = -(r.T() * z) / (d.T() * z)
-        alpha = (r.T() * p) / (p.T() * z)  # (np.dot(r , p)) / (np.dot(p , z))
+        d = r + beta * d
-        # print(alpha)
 
-        x = x + alpha * p  # neue Itterierte x
+        iterations = iterations + 1
-        r = r - alpha * z  # neues Residuum
+    return x
-
-        # 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()

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([-1, 2, -1], structure="tridiagonal", n=n)
-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))
+x0 = Vector([1] * n)
-f = Vector([h ** 2 * 2] * n)
+b = Vector([h**2 * 2] * n)
 
-# Toleranz epsilon
+x = cg.cg(A, x0, b)
-tol = 0.001
 
-cg.cg(n, A, f, tol)
+if rank == 0:
+    print(x)

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 = []
+        flattened_data, shape = Matrix.flatten_internal(matrices)
-        rows = 0
+        return Matrix(flattened_data, shape)
-        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))
 
     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]

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):
+        """
+        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())
+        return MatrixMPI(matrix)
-
-    def __str__(self):
-        return str(self.__data__)
 
     def shape(self):
         return self.__data__.shape()
 
-    def get_rank_submatrix(self):
+    def get_rank_subdata(self):
-        rows = len(self.__chunk__)
-        cols = self.__data__.shape()[1]
-        return Matrix(self.__data__[self.__chunk__], (rows, cols))
-
-    def get_matrix(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())
+        return MatrixMPI.of(Matrix.flatten(self.__mpi_comm__.allgather(-self.__rank_subdata__)))
-        data = self.__mpi_comm__.bcast(gathered_data)
-        return MatrixMPI.of(Matrix.flatten(data))
 
     def __add__(self, other):
         if isinstance(other, MatrixMPI):
-            other = other.get_rank_submatrix()
+            other = other.__rank_subdata__
-        gathered_data = self.__mpi_comm__.gather(self.get_rank_submatrix() + other)
+        return MatrixMPI.of(Matrix.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ + other)))
-        data = self.__mpi_comm__.bcast(gathered_data)
-        return MatrixMPI.of(Matrix.flatten(data))
 
     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)
+        return MatrixMPI.of(Matrix.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ / other)))
-        data = self.__mpi_comm__.bcast(gathered_data)
-        return MatrixMPI.of(Matrix.flatten(data))
 
     def __mul__(self, other):
         if isinstance(other, MatrixMPI):
-            other = other.get_matrix()
+            other = other.get_data()
-        gathered_data = self.__mpi_comm__.gather(self.get_rank_submatrix() * other)
+        return MatrixMPI.of(Matrix.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ * other)))
-        data = self.__mpi_comm__.bcast(gathered_data)
-        return MatrixMPI.of(Matrix.flatten(data))
 
     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):

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``

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):
+        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())
+        return VectorMPI(vector)
-
-    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
 
     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__
+            other = other.__rank_subdata__
-        return VectorMPI.of(self.__data__ + other)
+        return VectorMPI.of(Vector.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ + other)))
 
     def __mul__(self, other):
         if isinstance(other, VectorMPI):
             other = other.__data__
+
+        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__
+            other = other.__rank_subdata__
-        return VectorMPI.of(self.__data__ / other)
+        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())
+        return VectorMPI.of(Vector.flatten(self.__mpi_comm__.allgather(self.__rank_subdata__ / self.norm())))
-
-    def __getitem__(self, key):
-        return self.__data__[key]
-
-    def __setitem__(self, key, value):
-        self.__data__[key] = value