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


[docs]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}