diff --git a/doc/source/whatsnew/0.2.rst b/doc/source/whatsnew/0.2.rst index e9db472aa9df3ba2572e5009193e91b86935d2d4..81f43b600223d40c8dc6322ff2dd96ea0d6fe40b 100644 --- a/doc/source/whatsnew/0.2.rst +++ b/doc/source/whatsnew/0.2.rst @@ -10,13 +10,15 @@ local density of states. Return value of sparse solver ----------------------------- - `kwant.solvers.sparse.solve` now always returns a single instance of `~kwant.solvers.sparse.BlockResult`. The latter has been generalized to include more information for leads defined as infinite systems. Return value of `~kwant.system.System.hamiltonian_submatrix` ------------------------------------------------------------ - `~kwant.system.System.hamiltonian_submatrix` now always returns a 3-tuple ``hamiltonian_part, to_norb, from_norb``. + +Return value of `~kwant.solvers.sparse.make_linear_sys` has changed +------------------------------------------------------------------- +A namedtuple is used for more clarity. diff --git a/kwant/solvers/sparse.py b/kwant/solvers/sparse.py index cdaa1b31a1ba28a512a0c7d93b7e75ad6f6d31b7..8ce17685ee66743b44e392b8519349b777863f50 100644 --- a/kwant/solvers/sparse.py +++ b/kwant/solvers/sparse.py @@ -1,4 +1,4 @@ -__all__ = ['make_linear_sys', 'solve', 'ldos', 'BlockResult'] +__all__ = ['make_linear_sys', 'LinearSys', 'solve', 'ldos', 'BlockResult'] from functools import reduce from collections import namedtuple @@ -67,6 +67,10 @@ def factorized(A, piv_tol=1.0, sym_piv_tol=1.0): return solve + +LinearSys = namedtuple('LinearSys', ['h_sys', 'rhs', 'keep_vars']) + + def make_linear_sys(sys, out_leads, in_leads, energy=0, force_realspace=False): """ Make a sparse linear system of equations defining a scattering problem. @@ -89,15 +93,17 @@ def make_linear_sys(sys, out_leads, in_leads, energy=0, force_realspace=False): Returns ------- - (h_sys, rhs, keep_vars, lead_info) : tuple of inhomogeneous data + (h_sys, rhs, keep_vars) : LinearSys `h_sys` is a numpy.sparse.csc_matrix, containing the right hand side of the system of equations, `rhs` is the list of matrices with the left hand side, `keep_vars` is a list with numbers of variables in the solution that have to be stored (typically a small part of the - complete solution). Finally, a list `lead_info` is also returned, which - contains mode wave functions in each lead that is defined as a - tight-binding system (and as returned by `kwant.physics.modes`), and - the lead self-energy otherwise. + complete solution). + lead_info : list of objects + Contains one entry for each lead. For a lead defined as a + tight-binding system, this is an instance of `kwant.physics.Modes` (as + returned by `kwant.physics.modes`), otherwise the lead self-energy + matrix. Notes ----- @@ -197,11 +203,10 @@ def make_linear_sys(sys, out_leads, in_leads, energy=0, force_realspace=False): rhs.append(sp.coo_matrix((-np.ones(l), [indices, np.arange(l)]))) - result = (h_sys, rhs, keep_vars, lead_info) - return result + return LinearSys(h_sys, rhs, keep_vars), lead_info -def solve_linsys(a, b, keep_vars=None): +def solve_linear_sys(a, b, keep_vars=None): """ Solve matrix system of equations a x = b with sparse input. @@ -299,10 +304,11 @@ def solve(sys, energy=0, out_leads=None, in_leads=None, force_realspace=False): raise ValueError('Lead lists must be sorted and with unique entries.') if len(in_leads) == 0 or len(out_leads) == 0: raise ValueError('No output is requested.') - linsys = make_linear_sys(sys, out_leads, in_leads, energy, force_realspace) - out_modes = [len(i) for i in linsys[2]] - in_modes = [i.shape[1] for i in linsys[1]] - result = BlockResult(solve_linsys(*linsys[: -1]), linsys[3]) + linsys, lead_info = make_linear_sys(sys, out_leads, in_leads, energy, + force_realspace) + out_modes = [len(i) for i in linsys.keep_vars] + in_modes = [i.shape[1] for i in linsys.rhs] + result = BlockResult(solve_linear_sys(*linsys), lead_info) result.in_leads = in_leads result.out_leads = out_leads return result @@ -395,16 +401,14 @@ def ldos(fsys, e=0): local density of states at each orbital of the system. """ Modes = physics.Modes - linsys = make_linear_sys(fsys, [], [], e) - num_extra_vars = 0 - for i in linsys[-1]: - if isinstance(i, Modes): - num_extra_vars += i.vecs.shape[1] - i.nmodes - a = linsys[0] - slv = factorized(a) - num_orb = a.shape[0] - num_extra_vars - vec = np.zeros(a.shape[0], complex) + linsys, lead_info = make_linear_sys(fsys, [], [], e) + h = linsys.h_sys + num_extra_vars = sum(li.vecs.shape[1] - li.nmodes + for li in lead_info if isinstance(li, Modes)) + num_orb = h.shape[0] - num_extra_vars + vec = np.zeros(h.shape[0], complex) ldos = np.zeros(num_orb, complex) + slv = factorized(h) for i in xrange(num_orb): vec[i] = 1 ldos[i] = slv(vec)[i]