Coverage for gpaw/test/matrix/test_eighg.py: 100%

1import numpy as np 

2import pytest 

3from gpaw.core.matrix import Matrix 

4from gpaw.mpi import broadcast_exception, world 

5import scipy.linalg as linalg 

6 

7 

8@pytest.mark.parametrize('dtype', [float, complex]) 

9def test_eighg(dtype): 

10 n = 5 

11 S0 = Matrix(n, n, dist=(world, 1, 1), dtype=dtype) 

12 S = S0.new(dist=(world, world.size, 1)) 

13 H0 = S0.new() 

14 H = S.new() 

15 if world.rank == 0: 

16 S0.data[:] = np.eye(n) 

17 H0.data[:] = 0.0 

18 H0.data.ravel()[::n + 1] = np.arange(n) + 1 

19 if dtype == float: 

20 S0.data[-1, 0] = 0.001 

21 H0.data[-1, 0] = 0.001 

22 else: 

23 S0.data[-1, 0] = 0.001j 

24 H0.data[-1, 0] = 0.001j 

25 H0.tril2full() 

26 S0.tril2full() 

27 H00 = H0.copy() 

28 S0.redist(S) 

29 H0.redist(H) 

30 L = S.copy() 

31 L.invcholesky() 

32 if world.rank == 0: 

33 eigs0, C0 = linalg.eigh(H0.data, S0.data) 

34 print(eigs0) 

35 print(C0) 

36 error = H0.data @ C0 - S0.data @ C0 @ np.diag(eigs0) 

37 print(error) 

38 

39 L0 = S0.new() 

40 L.redist(L0) 

41 eigs = H.eighg(L0) 

42 H.redist(H0) 

43 print(eigs) 

44 print(H0.data) 

45 with broadcast_exception(world): 

46 if world.rank == 0: 

47 assert abs(eigs - eigs0).max() < 1e-14 

48 error = H00.data @ H0.data - S0.data @ H0.data @ np.diag(eigs) 

49 assert abs(error).max() < 1e-14