实现高斯约旦消元法

from .Matrix import Matrix
from .Vector import Vector


class LinearSystem:

    def __init__(self, A, b):

        assert A.row_num() == len(b), "row number of A must be equal to the length of b"
        self._m = A.row_num()
        self._n = A.col_num()
        assert self._m == self._n  # TODO: no this restriction

        self.Ab = [Vector(A.row_vector(i).underlying_list() + [b[i]])
                   for i in range(self._m)]

    def _max_row(self, index_i, index_j, n):

        best, ret = abs(self.Ab[index_i][index_j]), index_i
        for i in range(index_i + 1, n):
            if abs(self.Ab[i][index_j]) > best:
                best, ret = abs(self.Ab[i][index_j]), i
        return ret

    def _forward(self):

        n = self._m
        for i in range(n):
            # Ab[i][i]为主元
            max_row = self._max_row(i, n)
            self.Ab[i], self.Ab[max_row] = self.Ab[max_row], self.Ab[i]

            # 将主元归为一
            self.Ab[i] = self.Ab[i] / self.Ab[i][i]  # TODO: self.Ab[i][i] == 0?
            for j in range(i + 1, n):
                self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i]

    def _backward(self):

        n = self._m
        for i in range(n - 1, -1, -1):
            # Ab[i][i]为主元
            for j in range(i - 1, -1, -1):
                self.Ab[j] = self.Ab[j] - self.Ab[j][i] * self.Ab[i]

    def gauss_jordan_elimination(self):

        self._forward()
        self._backward()

    def fancy_print(self):

        for i in range(self._m):
            print(" ".join(str(self.Ab[i][j]) for j in range(self._n)), end=" ")
            print("|", self.Ab[i][-1])

测试

from playLA.Matrix import Matrix
from playLA.Vector import Vector
from playLA.LinearSystem import LinearSystem


if __name__ == "__main__":

    A = Matrix([[1, 2, 4], [3, 7, 2], [2, 3, 3]])
    b = Vector([7, -11, 1])
    ls = LinearSystem(A, b)
    ls.gauss_jordan_elimination()
    ls.fancy_print()
    print()
    #
    # A2 = Matrix([[1, -3, 5], [2, -1, -3], [3, 1, 4]])
    # b2 = Vector([-9, 19, -13])
    # ls2 = LinearSystem(A2, b2)
    # ls2.gauss_jordan_elimination()
    # ls2.fancy_print()
    # print()
    #
    # A3 = Matrix([[1, 2, -2], [2, -3, 1], [3, -1, 3]])
    # b3 = Vector([6, -10, -16])
    # ls3 = LinearSystem(A3, b3)
    # ls3.gauss_jordan_elimination()
    # ls3.fancy_print()
    # print()
    #
    # A4 = Matrix([[3, 1, -2], [5, -3, 10], [7, 4, 16]])
    # b4 = Vector([4, 32, 13])
    # ls4 = LinearSystem(A4, b4)
    # ls4.gauss_jordan_elimination()
    # ls4.fancy_print()
    # print()
    #
    # A5 = Matrix([[6, -3, 2], [5, 1, 12], [8, 5, 1]])
    # b5 = Vector([31, 36, 11])
    # ls5 = LinearSystem(A5, b5)
    # ls5.gauss_jordan_elimination()
    # ls5.fancy_print()
    # print()
    #
    # A6 = Matrix([[1, 1, 1], [1, -1, -1], [2, 1, 5]])
    # b6 = Vector([3, -1, 8])
    # ls6 = LinearSystem(A6, b6)
    # ls6.gauss_jordan_elimination()
    # ls6.fancy_print()
    # print()