123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136 |
- from numpy.testing import assert_array_almost_equal, assert_allclose, assert_
- from numpy import (array, eye, zeros, empty_like, empty, tril_indices_from,
- tril, triu_indices_from, spacing, float32, float64,
- complex64, complex128)
- from numpy.random import rand, randint, seed
- from scipy.linalg import ldl
- import pytest
- from pytest import raises as assert_raises, warns
- from numpy import ComplexWarning
- def test_args():
- A = eye(3)
- # Nonsquare array
- assert_raises(ValueError, ldl, A[:, :2])
- # Complex matrix with imaginary diagonal entries with "hermitian=True"
- with warns(ComplexWarning):
- ldl(A*1j)
- def test_empty_array():
- a = empty((0, 0), dtype=complex)
- l, d, p = ldl(empty((0, 0)))
- assert_array_almost_equal(l, empty_like(a))
- assert_array_almost_equal(d, empty_like(a))
- assert_array_almost_equal(p, array([], dtype=int))
- def test_simple():
- a = array([[-0.39-0.71j, 5.14-0.64j, -7.86-2.96j, 3.80+0.92j],
- [5.14-0.64j, 8.86+1.81j, -3.52+0.58j, 5.32-1.59j],
- [-7.86-2.96j, -3.52+0.58j, -2.83-0.03j, -1.54-2.86j],
- [3.80+0.92j, 5.32-1.59j, -1.54-2.86j, -0.56+0.12j]])
- b = array([[5., 10, 1, 18],
- [10., 2, 11, 1],
- [1., 11, 19, 9],
- [18., 1, 9, 0]])
- c = array([[52., 97, 112, 107, 50],
- [97., 114, 89, 98, 13],
- [112., 89, 64, 33, 6],
- [107., 98, 33, 60, 73],
- [50., 13, 6, 73, 77]])
- d = array([[2., 2, -4, 0, 4],
- [2., -2, -2, 10, -8],
- [-4., -2, 6, -8, -4],
- [0., 10, -8, 6, -6],
- [4., -8, -4, -6, 10]])
- e = array([[-1.36+0.00j, 0+0j, 0+0j, 0+0j],
- [1.58-0.90j, -8.87+0j, 0+0j, 0+0j],
- [2.21+0.21j, -1.84+0.03j, -4.63+0j, 0+0j],
- [3.91-1.50j, -1.78-1.18j, 0.11-0.11j, -1.84+0.00j]])
- for x in (b, c, d):
- l, d, p = ldl(x)
- assert_allclose(l.dot(d).dot(l.T), x, atol=spacing(1000.), rtol=0)
- u, d, p = ldl(x, lower=False)
- assert_allclose(u.dot(d).dot(u.T), x, atol=spacing(1000.), rtol=0)
- l, d, p = ldl(a, hermitian=False)
- assert_allclose(l.dot(d).dot(l.T), a, atol=spacing(1000.), rtol=0)
- u, d, p = ldl(a, lower=False, hermitian=False)
- assert_allclose(u.dot(d).dot(u.T), a, atol=spacing(1000.), rtol=0)
- # Use upper part for the computation and use the lower part for comparison
- l, d, p = ldl(e.conj().T, lower=0)
- assert_allclose(tril(l.dot(d).dot(l.conj().T)-e), zeros((4, 4)),
- atol=spacing(1000.), rtol=0)
- def test_permutations():
- seed(1234)
- for _ in range(10):
- n = randint(1, 100)
- # Random real/complex array
- x = rand(n, n) if randint(2) else rand(n, n) + rand(n, n)*1j
- x = x + x.conj().T
- x += eye(n)*randint(5, 1e6)
- l_ind = tril_indices_from(x, k=-1)
- u_ind = triu_indices_from(x, k=1)
- # Test whether permutations lead to a triangular array
- u, d, p = ldl(x, lower=0)
- # lower part should be zero
- assert_(not any(u[p, :][l_ind]), 'Spin {} failed'.format(_))
- l, d, p = ldl(x, lower=1)
- # upper part should be zero
- assert_(not any(l[p, :][u_ind]), 'Spin {} failed'.format(_))
- @pytest.mark.parametrize("dtype", [float32, float64])
- @pytest.mark.parametrize("n", [30, 150])
- def test_ldl_type_size_combinations_real(n, dtype):
- seed(1234)
- msg = ("Failed for size: {}, dtype: {}".format(n, dtype))
- x = rand(n, n).astype(dtype)
- x = x + x.T
- x += eye(n, dtype=dtype)*dtype(randint(5, 1e6))
- l, d1, p = ldl(x)
- u, d2, p = ldl(x, lower=0)
- rtol = 1e-4 if dtype is float32 else 1e-10
- assert_allclose(l.dot(d1).dot(l.T), x, rtol=rtol, err_msg=msg)
- assert_allclose(u.dot(d2).dot(u.T), x, rtol=rtol, err_msg=msg)
- @pytest.mark.parametrize("dtype", [complex64, complex128])
- @pytest.mark.parametrize("n", [30, 150])
- def test_ldl_type_size_combinations_complex(n, dtype):
- seed(1234)
- msg1 = ("Her failed for size: {}, dtype: {}".format(n, dtype))
- msg2 = ("Sym failed for size: {}, dtype: {}".format(n, dtype))
- # Complex hermitian upper/lower
- x = (rand(n, n)+1j*rand(n, n)).astype(dtype)
- x = x+x.conj().T
- x += eye(n, dtype=dtype)*dtype(randint(5, 1e6))
- l, d1, p = ldl(x)
- u, d2, p = ldl(x, lower=0)
- rtol = 1e-4 if dtype is complex64 else 1e-10
- assert_allclose(l.dot(d1).dot(l.conj().T), x, rtol=rtol, err_msg=msg1)
- assert_allclose(u.dot(d2).dot(u.conj().T), x, rtol=rtol, err_msg=msg1)
- # Complex symmetric upper/lower
- x = (rand(n, n)+1j*rand(n, n)).astype(dtype)
- x = x+x.T
- x += eye(n, dtype=dtype)*dtype(randint(5, 1e6))
- l, d1, p = ldl(x, hermitian=0)
- u, d2, p = ldl(x, lower=0, hermitian=0)
- assert_allclose(l.dot(d1).dot(l.T), x, rtol=rtol, err_msg=msg2)
- assert_allclose(u.dot(d2).dot(u.T), x, rtol=rtol, err_msg=msg2)
|