## "Fossies" - the Fresh Open Source Software Archive

### Source code changes of the file "mpmath/tests/test_linalg.py" betweenmpmath-1.0.0.tar.gz and mpmath-1.1.0.tar.gz

About: mpmath is a Python library for arbitrary-precision floating-point arithmetic.

test_linalg.py  (mpmath-1.0.0):test_linalg.py  (mpmath-1.1.0)
# TODO: don't use round # TODO: don't use round
from __future__ import division from __future__ import division
import pytest
from mpmath import * from mpmath import *
xrange = libmp.backend.xrange xrange = libmp.backend.xrange
# XXX: these shouldn't be visible(?) # XXX: these shouldn't be visible(?)
LU_decomp = mp.LU_decomp LU_decomp = mp.LU_decomp
L_solve = mp.L_solve L_solve = mp.L_solve
U_solve = mp.U_solve U_solve = mp.U_solve
householder = mp.householder householder = mp.householder
improve_solution = mp.improve_solution improve_solution = mp.improve_solution
skipping to change at line 135 skipping to change at line 136
H, p, x, r = householder(extend(A, y)) H, p, x, r = householder(extend(A, y))
x = matrix(x) x = matrix(x)
y = matrix(y) y = matrix(y)
residuals.append(norm(r, 2)) residuals.append(norm(r, 2))
refres.append(norm(residual(A, x, y), 2)) refres.append(norm(residual(A, x, y), 2))
assert [round(res, 10) for res in residuals] == [15.1733888877, assert [round(res, 10) for res in residuals] == [15.1733888877,
0.82378073210000002, 0.302645887, 0.0260109244, 0.82378073210000002, 0.302645887, 0.0260109244,
0.00058653999999999998] 0.00058653999999999998]
assert norm(matrix(residuals) - matrix(refres), inf) < 1.e-13 assert norm(matrix(residuals) - matrix(refres), inf) < 1.e-13
def hilbert_cmplx(n):
# Complexified Hilbert matrix
A = hilbert(2*n,n)
v = randmatrix(2*n, 2, min=-1, max=1)
v = v.apply(lambda x: exp(1J*pi()*x))
A = diag(v[:,0])*A*diag(v[:n,1])
return A
residuals_cmplx = []
refres_cmplx = []
for n in range(2, 10):
A = hilbert_cmplx(n)
H, p, x, r = householder(A.copy())
residuals_cmplx.append(norm(r, 2))
refres_cmplx.append(norm(residual(A[:,:n-1], x, A[:,n-1]), 2))
assert norm(matrix(residuals_cmplx) - matrix(refres_cmplx), inf) < 1.e-13
def test_factorization(): def test_factorization():
A = randmatrix(5) A = randmatrix(5)
P, L, U = lu(A) P, L, U = lu(A)
assert mnorm(P*A - L*U, 1) < 1.e-15 assert mnorm(P*A - L*U, 1) < 1.e-15
def test_solve(): def test_solve():
assert norm(residual(A6, lu_solve(A6, b6), b6), inf) < 1.e-10 assert norm(residual(A6, lu_solve(A6, b6), b6), inf) < 1.e-10
assert norm(residual(A7, lu_solve(A7, b7), b7), inf) < 1.5 assert norm(residual(A7, lu_solve(A7, b7), b7), inf) < 1.5
assert norm(residual(A8, lu_solve(A8, b8), b8), inf) <= 3 + 1.e-10 assert norm(residual(A8, lu_solve(A8, b8), b8), inf) <= 3 + 1.e-10
assert norm(residual(A6, qr_solve(A6, b6)[0], b6), inf) < 1.e-10 assert norm(residual(A6, qr_solve(A6, b6)[0], b6), inf) < 1.e-10
skipping to change at line 160 skipping to change at line 178
def test_solve_overdet_complex(): def test_solve_overdet_complex():
A = matrix([[1, 2j], [3, 4j], [5, 6]]) A = matrix([[1, 2j], [3, 4j], [5, 6]])
b = matrix([1 + j, 2, -j]) b = matrix([1 + j, 2, -j])
assert norm(residual(A, lu_solve(A, b), b)) < 1.0208 assert norm(residual(A, lu_solve(A, b), b)) < 1.0208
def test_singular(): def test_singular():
mp.dps = 15 mp.dps = 15
A = [[5.6, 1.2], [7./15, .1]] A = [[5.6, 1.2], [7./15, .1]]
B = repr(zeros(2)) B = repr(zeros(2))
b = [1, 2] b = [1, 2]
for i in ['lu_solve(%s, %s)' % (A, b), 'lu_solve(%s, %s)' % (B, b), for i in ['lu_solve(%s, %s)' % (A, b), 'lu_solve(%s, %s)' % (B, b),
'qr_solve(%s, %s)' % (A, b), 'qr_solve(%s, %s)' % (B, b)]: 'qr_solve(%s, %s)' % (A, b), 'qr_solve(%s, %s)' % (B, b)]:
) pytest.raises((ZeroDivisionError, ValueError), lambda: eval(i))
def test_cholesky(): def test_cholesky():
assert fp.cholesky(fp.matrix(A9)) == fp.matrix([[2, 0, 0], [1, 2, 0], [-1, - 3/2, 3/2]]) assert fp.cholesky(fp.matrix(A9)) == fp.matrix([[2, 0, 0], [1, 2, 0], [-1, - 3/2, 3/2]])
x = fp.cholesky_solve(A9, b9) x = fp.cholesky_solve(A9, b9)
assert fp.norm(fp.residual(A9, x, b9), fp.inf) == 0 assert fp.norm(fp.residual(A9, x, b9), fp.inf) == 0
def test_det(): def test_det():
assert det(A1) == 1 assert det(A1) == 1
assert round(det(A2), 14) == 8 assert round(det(A2), 14) == 8
assert round(det(A3)) == 1834 assert round(det(A3)) == 1834
skipping to change at line 301 skipping to change at line 313
Q, R = qr(A, mode, edps = exdps) Q, R = qr(A, mode, edps = exdps)
#print('\n\n A = \n', nstr(A, 4)) #print('\n\n A = \n', nstr(A, 4))
#print('\n Q = \n', nstr(Q, 4)) #print('\n Q = \n', nstr(Q, 4))
#print('\n R = \n', nstr(R, 4)) #print('\n R = \n', nstr(R, 4))
#print('\n Q*R = \n', nstr(Q*R, 4)) #print('\n Q*R = \n', nstr(Q*R, 4))
maxnorm = mpf('1.0E-11') maxnorm = mpf('1.0E-11')
n1 = norm(A - Q * R) n1 = norm(A - Q * R)
#print '\n Norm of A - Q * R = ', n1 #print '\n Norm of A - Q * R = ', n1
n1 assert n1 <= maxnorm
if dtype == 'real': if dtype == 'real':
n1 = norm(eye(m) - Q.T * Q) n1 = norm(eye(m) - Q.T * Q)
#print ' Norm of I - Q.T * Q = ', n1 #print ' Norm of I - Q.T * Q = ', n1
n1 assert n1 <= maxnorm
n1 = norm(eye(m) - Q * Q.T) n1 = norm(eye(m) - Q * Q.T)
#print ' Norm of I - Q * Q.T = ', n1 #print ' Norm of I - Q * Q.T = ', n1
n1 assert n1 <= maxnorm
if dtype == 'complex': if dtype == 'complex':
n1 = norm(eye(m) - Q.T * Q.conjugate()) n1 = norm(eye(m) - Q.T * Q.conjugate())
#print ' Norm of I - Q.T * Q.conjugate() = ', n1 #print ' Norm of I - Q.T * Q.conjugate() = ', n1
n1 assert n1 <= maxnorm
n1 = norm(eye(m) - Q.conjugate() * Q.T) n1 = norm(eye(m) - Q.conjugate() * Q.T)
#print ' Norm of I - Q.conjugate() * Q.T = ', n1 #print ' Norm of I - Q.conjugate() * Q.T = ', n1
n1 assert n1 <= maxnorm
End of changes. 9 change blocks.
15 lines changed or deleted 23 lines changed or added