Source code for tryalgo.gauss_jordan

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""\
Linear equation system Ax=b by Gauss-Jordan

jill-jenn vie et christoph durr - 2014-2018
"""

__all__ = ["gauss_jordan", "GJ_ZERO_SOLUTIONS", "GJ_SINGLE_SOLUTION",
           "GJ_SEVERAL_SOLUTIONS"]


# snip{
# pylint: disable=chained-comparison
def is_zero(x):                    # tolerance
    """error tolerant zero test
    """
    return -1e-6 < x and x < 1e-6
    # replace with x == 0 si we are handling Fraction elements


GJ_ZERO_SOLUTIONS = 0
GJ_SINGLE_SOLUTION = 1
GJ_SEVERAL_SOLUTIONS = 2




def gauss_jordan(A, x, b):
    """Linear equation system Ax=b by Gauss-Jordan

    :param A: m by n matrix
    :param x: table of size n
    :param b: table of size m
    :modifies: x will contain solution if any
    :returns int:
          0 if no solution,
          1 if solution unique,
          2 otherwise
    :complexity: :math:`O(n^2m)`
    """
    n = len(x)
    m = len(b)
    assert len(A) == m and len(A[0]) == n
    S = []                        # put linear system in a single matrix S
    for i in range(m):
        S.append(A[i][:] + [b[i]])
    S.append(list(range(n)))      # indices in x
    k = diagonalize(S, n, m)
    if k < m:
        for i in range(k, m):
            if not is_zero(S[i][n]):
                return GJ_ZERO_SOLUTIONS
    for j in range(k):
        x[S[m][j]] = S[j][n]
    if k < n:
        for j in range(k, n):
            x[S[m][j]] = 0
        return GJ_SEVERAL_SOLUTIONS
    return GJ_SINGLE_SOLUTION



def diagonalize(S, n, m):
    """diagonalize """
    for k in range(min(n, m)):
        val, i, j = max((abs(S[i][j]), i, j)
                        for i in range(k, m) for j in range(k, n))
        if is_zero(val):
            return k
        S[i], S[k] = S[k], S[i]    # swap lines k, i
        for r in range(m + 1):     # swap columns k, j
            S[r][j], S[r][k] = S[r][k], S[r][j]
        pivot = float(S[k][k])     # without float if Fraction elements
        for j in range(k, n + 1):
            S[k][j] /= pivot       # divide line k by pivot
        for i in range(m):         # remove line k scaled by line i
            if i != k:
                fact = S[i][k]
                for j in range(k, n + 1):
                    S[i][j] -= fact * S[k][j]
    return min(n, m)
# snip}