diff --git a/mpmath/matrices/linalg.py b/mpmath/matrices/linalg.py index 2e9503cc..a6a3571e 100644 --- a/mpmath/matrices/linalg.py +++ b/mpmath/matrices/linalg.py @@ -134,6 +134,9 @@ class LinearAlgebraMethods: if current > biggest: # TODO: what if equal? biggest = current p[j] = k + # without pivot LU fails + if p[j] is None: + raise ZeroDivisionError('matrix is numerically singular') # swap rows according to p ctx.swap_row(A, j, p[j]) if ctx.absmin(A[j,j]) <= tol: @@ -539,8 +542,6 @@ class LinearAlgebraMethods: try: # do not overwrite A A = ctx.matrix(A).copy() - if any(all(ctx.almosteq(0, x) for x in A[:, c]) for c in range(A.cols)): - return 0 * A[0, 0] # use LU factorization to calculate determinant try: R, p = ctx.LU_decomp(A) diff --git a/mpmath/tests/test_linalg.py b/mpmath/tests/test_linalg.py index 86e92c3d..b9dc72a9 100644 --- a/mpmath/tests/test_linalg.py +++ b/mpmath/tests/test_linalg.py @@ -69,10 +69,15 @@ A10 = matrix([[1.0 + 1.0j, 2.0, 2.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]) b10 = [1.0, 1.0 + 1.0j, 1.0] + A11 = matrix([[4, 0, -2], [2, 0, -4], [2, 0, 5.5]]) +A12 = matrix([[1,0,0], + [0,1,0], + [0,0,1.0j]]) + def test_LU_decomp(): A = A3.copy() @@ -96,6 +101,12 @@ def test_LU_decomp(): LU_decomp(A, overwrite=1) assert A != bak + try: + LU_decomp(A11) + assert False + except ZeroDivisionError: + assert True + def test_inverse(): for A in [A1, A2, A5]: inv = inverse(A) @@ -197,7 +208,7 @@ def test_det(): assert round(det(A6)) == 78356463 assert det(zeros(3)) == 0 assert det(A11) == 0 - + assert absmin(det(A12*1e-30) - 1e-30) < eps def test_cond(): A = matrix([[1.2969, 0.8648], [0.2161, 0.1441]])