Coverage for gpaw/test/parallel/test_mkl_scalapack_eig.py: 32%

1from gpaw.utilities.scalapack import mkl_scalapack_diagonalize_non_symmetric 

2from gpaw.utilities.scalapack import have_mkl 

3from gpaw.blacs import BlacsGrid, Redistributor 

4from gpaw.mpi import world 

5from gpaw.matrix import suggest_blocking 

6 

7import numpy as np 

8 

9import pytest 

10 

11 

12def test_mkl_eig(): 

13 if not have_mkl(): 

14 pytest.skip('MKL Scalapack functions only testable with ' 

15 'the intelMKL toolchain active and ' 

16 'intelmkl = True in siteconfig.py') 

17 

18 gsize = 200 

19 

20 scgrid = BlacsGrid(world, 1, 1) 

21 scdesc = scgrid.new_descriptor(gsize, gsize, gsize, gsize) 

22 asc = scdesc.empty(dtype=complex) 

23 

24 if world.rank == 0: 

25 np.random.seed(1337) 

26 asc[:] = np.random.rand(gsize, gsize) 

27 

28 ga, gb, bsize = suggest_blocking(gsize, world.size) 

29 grid = BlacsGrid(world, ga, gb) 

30 desc = grid.new_descriptor(gsize, gsize, bsize, bsize) 

31 redist_sc_to_full = Redistributor(world, scdesc, desc) 

32 a = redist_sc_to_full.redistribute(asc) 

33 z = desc.empty(dtype=complex) 

34 

35 eps = np.empty(gsize, dtype=complex) 

36 mkl_scalapack_diagonalize_non_symmetric(desc, a, z, eps) 

37 

38 redist_full_to_sc = Redistributor(world, desc, scdesc) 

39 zsc = redist_full_to_sc.redistribute(z) 

40 

41 if world.rank == 0: 

42 assert np.linalg.norm( 

43 asc - zsc @ np.diag(eps) @ np.linalg.inv(zsc) 

44 ) == pytest.approx(0, abs=1e-8) 

45 

46 # Test against lapack 

47 [eps2, z2] = np.linalg.eig(asc) 

48 

49 assert np.sum(np.abs(np.subtract.outer(eps, eps2)) < 1e-10) == gsize