From a5c24b91577c3c8f3e711aafe288078363487000 Mon Sep 17 00:00:00 2001
From: Michael Wimmer <wimmer@lorentz.leidenuniv.nl>
Date: Tue, 15 May 2012 17:28:18 +0200
Subject: [PATCH] add new solver based on MUMPS; reorganization of sparse
solver code in general
---
TODO.txt | 4 -
doc/source/images/tutorial1a.py | 2 +-
doc/source/reference/kwant.solvers.common.rst | 10 +
doc/source/reference/kwant.solvers.mumps.rst | 9 +
doc/source/reference/kwant.solvers.rst | 77 ++-
doc/source/reference/kwant.solvers.sparse.rst | 7 +-
kwant/__init__.py | 9 +-
kwant/physics/noise.py | 26 +
kwant/physics/tests/test_noise.py | 60 ++
kwant/solvers/common.py | 472 +++++++++++++++
kwant/solvers/mumps.py | 163 +++++
kwant/solvers/sparse.py | 562 +++++-------------
kwant/solvers/tests/_test_sparse.py | 364 ++++++++++++
kwant/solvers/tests/test_mumps.py | 92 +++
kwant/solvers/tests/test_sparse.py | 292 +--------
15 files changed, 1447 insertions(+), 702 deletions(-)
create mode 100644 doc/source/reference/kwant.solvers.common.rst
create mode 100644 doc/source/reference/kwant.solvers.mumps.rst
create mode 100644 kwant/physics/noise.py
create mode 100644 kwant/physics/tests/test_noise.py
create mode 100644 kwant/solvers/common.py
create mode 100644 kwant/solvers/mumps.py
create mode 100644 kwant/solvers/tests/_test_sparse.py
create mode 100644 kwant/solvers/tests/test_mumps.py
diff --git a/TODO.txt b/TODO.txt
index 3a9a154..42d1050 100644
--- a/TODO.txt
+++ b/TODO.txt
@@ -17,10 +17,6 @@ Roughly in order of importance. -*-org-*-
* Provide nice support for graphene double layers
This could be done by allowing lattices to be shifted, or in some other way.
-* Wrap mumps, umfpack, or some other sparse linear algebra library with Cython.
- Use it directly in sparse solver. This will allow to fine-tune the solution
- of sparse systems.
-
* Benchmark mumps and check whether nested dissection would be useful.
If yes, implement it.
diff --git a/doc/source/images/tutorial1a.py b/doc/source/images/tutorial1a.py
index 0e14b6e..41c5fef 100644
--- a/doc/source/images/tutorial1a.py
+++ b/doc/source/images/tutorial1a.py
@@ -93,7 +93,7 @@ for ie in xrange(100):
energy = ie * 0.01
# compute the scattering matrix at energy energy
- smatrix = kwant.solvers.sparse.solve(fsys, energy)
+ smatrix = kwant.solve(fsys, energy)
# compute the transmission probability from lead 0 to
# lead 1
diff --git a/doc/source/reference/kwant.solvers.common.rst b/doc/source/reference/kwant.solvers.common.rst
new file mode 100644
index 0000000..282b8dd
--- /dev/null
+++ b/doc/source/reference/kwant.solvers.common.rst
@@ -0,0 +1,10 @@
+:mod:`kwant.solvers.common` -- common functionality used by solvers
+===================================================================
+
+.. automodule:: kwant.solvers.common
+
+.. autosummary::
+ :toctree: generated/
+
+ BlockResult
+ SparseSolver
diff --git a/doc/source/reference/kwant.solvers.mumps.rst b/doc/source/reference/kwant.solvers.mumps.rst
new file mode 100644
index 0000000..418e54a
--- /dev/null
+++ b/doc/source/reference/kwant.solvers.mumps.rst
@@ -0,0 +1,9 @@
+:mod:`kwant.solvers.mumps` -- High performance sparse solver based on MUMPS
+===========================================================================
+
+.. automodule:: kwant.solvers.mumps
+
+.. autosummary::
+ :toctree: generated/
+
+ Solver
diff --git a/doc/source/reference/kwant.solvers.rst b/doc/source/reference/kwant.solvers.rst
index bde4ca8..2282ec8 100644
--- a/doc/source/reference/kwant.solvers.rst
+++ b/doc/source/reference/kwant.solvers.rst
@@ -1,8 +1,79 @@
:mod:`kwant.solvers` -- Library of solvers
==========================================
-The following solvers are available. Note that the solvers (with the exception
-of `kwant.solvers.sparse`) have to be imported explicitly.
+Overview of the solver infrastructure
+-------------------------------------
+
+kwant offers a variety of solvers for computing the solution
+to a quantum transport problem. These different solvers may either
+depend on different external libraries, or use very different internal
+algorithms for doing the actual computation.
+
+Fortunately, all these different solvers still use the same interface
+conventions. Thus, all solvers can be used without explicit knowledge
+of the internals. On the other hand, many solvers allow for some
+controls that can affect performance.
+
+All of the solver modules implement the following functions:
+
+- `~kwant.solvers.common.SparseSolver.solve`: Compute the Scattering matrix
+ or Green's function between different leads. Can be used to compute
+ conductances
+- `~kwant.solvers.common.SparseSolver.ldos`: Compute the local density of
+ states of the system.
+
+For details see the respective function.
+
+Due to the common interface you can thus write code that allows you to
+switch solvers by changing one `import`-line::
+
+ from kwant.solvers.some_solver import solve
+
+ ...
+
+ smatrix = solve(fsys)
+
+In addition to the common interface, some solvers allow for additional
+control parameters for performance-tuning. In this case they
+have a function `options` on the module level. Code could then
+look like this::
+
+ import kwant.solvers.some_solver as solver
+
+ solver.options(...)
+
+ ...
+
+ smatrix - solver.solve(fsys)
+
+Right now, the following solvers are implemented in kwant:
+
+- `~kwant.solvers.sparse`: A solver based on solving a sparse linear
+ system, using the direct sparse solvers provided by scipy.
+- `~kwant.solvers.mumps`: A solver based on solving a sparse linear
+ system using the direct sparse solver MUMPS. To use it, the MUMPS
+ library must be installed on the system. This solver typically
+ gives the best performance of all sparse solvers for a single core.
+ it allows for various solve options using `~kwant.solvers.mumps.options`.
+
+The default solver
+------------------
+
+Since computing conductance is the most basic task of kwant, the
+`solve`-function of one of the solvers is provided via `kwant.solve`.
+kwant chooses the solver which it considers best amongst the
+available solvers. You can see by calling
+
+>>> help(kwant.solve)
+
+from whoch module it has been imported.
+
+List of solver modules
+----------------------
+
+The modules of the solver package are listed in detail below. Note
+that the solvers (with the exception of the module providing
+`kwant.solve`) have to be imported explicitly.
.. module:: kwant.solvers
@@ -10,3 +81,5 @@ of `kwant.solvers.sparse`) have to be imported explicitly.
:maxdepth: 1
kwant.solvers.sparse
+ kwant.solvers.mumps
+ kwant.solvers.common
diff --git a/doc/source/reference/kwant.solvers.sparse.rst b/doc/source/reference/kwant.solvers.sparse.rst
index f0ea19a..ff10c25 100644
--- a/doc/source/reference/kwant.solvers.sparse.rst
+++ b/doc/source/reference/kwant.solvers.sparse.rst
@@ -1,12 +1,9 @@
:mod:`kwant.solvers.sparse` -- Basic sparse matrix solver
=========================================================
-.. module:: kwant.solvers.sparse
+.. automodule:: kwant.solvers.sparse
.. autosummary::
:toctree: generated/
- solve
- ldos
- make_linear_sys
- BlockResult
+ Solver
diff --git a/kwant/__init__.py b/kwant/__init__.py
index 63c4933..3b34cce 100644
--- a/kwant/__init__.py
+++ b/kwant/__init__.py
@@ -20,5 +20,10 @@ except:
else:
__all__.extend(['plotter', 'plot'])
-from .solvers.sparse import solve
-__all__.append('solve')
+# Mumps usally works best, use scipy_sparse as fallback
+try:
+ from .solvers.mumps import solve
+ __all__.append('solve')
+except ImportError:
+ from .solvers.sparse import solve
+ __all__.append('solve')
diff --git a/kwant/physics/noise.py b/kwant/physics/noise.py
new file mode 100644
index 0000000..9809b4e
--- /dev/null
+++ b/kwant/physics/noise.py
@@ -0,0 +1,26 @@
+import numpy as np
+
+def two_terminal_shotnoise(smatrix):
+ """Compute the shot-noise in a two-terminal setup.
+ [Given by tr((1 - t*t^\dagger) * t*t^\dagger)]."""
+
+ if len(smatrix.lead_info) != 2:
+ raise ValueError("Only works for two-terminal systems!")
+
+ if 1 in smatrix.out_leads and 0 in smatrix.in_leads:
+ ttdag = smatrix._a_ttdagger_a_inv(1, 0)
+ if 0 in smatrix.out_leads and 1 in smatrix.in_leads:
+ ttdag = smatrix._a_ttdagger_a_inv(0, 1)
+ else:
+ raise ValueError("Need S-matrix block for transmission!")
+
+ ttdag -= np.dot(ttdag, ttdag)
+ return np.trace(ttdag).real
+
+
+# A general multi-terminal routine for noise would need to also have the
+# voltages at various leads as input. (See
+# http://arxiv.org/abs/cond-mat/9910158) It could still be based on
+# smatrix._a_ttdagger_a_inv, i.e. be also valid also for self-energy leads,
+# provided that only true transmission blocks are used As long as nobody needs
+# it though, it does make little sense to make such a routine.
diff --git a/kwant/physics/tests/test_noise.py b/kwant/physics/tests/test_noise.py
new file mode 100644
index 0000000..11ccea1
--- /dev/null
+++ b/kwant/physics/tests/test_noise.py
@@ -0,0 +1,60 @@
+import numpy as np
+from nose.tools import assert_raises
+from numpy.testing import assert_almost_equal
+import kwant
+from kwant.physics.noise import two_terminal_shotnoise
+
+n = 5
+chain = kwant.lattice.Chain()
+
+def _twoterminal_system():
+ np.random.seed(11)
+ system = kwant.Builder()
+ lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
+ h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ h += h.conjugate().transpose()
+ h *= 0.8
+ t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+ lead[chain(0)] = h
+ system[chain(0)] = h * 1.2
+ lead[chain(0), chain(1)] = t
+ system.attach_lead(lead)
+ system.attach_lead(lead.reversed())
+ return system.finalized()
+
+
+def test_twoterminal_input():
+ """Input checks for two_terminal_shotnoise"""
+
+ fsys = _twoterminal_system()
+ sol = kwant.solve(fsys, out_leads=[0], in_leads=[0])
+ assert_raises(ValueError, two_terminal_shotnoise, sol)
+
+
+def test_twoterminal():
+ """Shot noise in a two-terminal conductor"""
+
+ fsys = _twoterminal_system()
+
+ sol = kwant.solve(fsys)
+ t = sol.submatrix(1, 0)
+ Tn = np.linalg.eigvalsh(np.dot(t, t.conj().T))
+ noise_should_be = np.sum(Tn * (1 - Tn))
+
+ assert_almost_equal(noise_should_be, two_terminal_shotnoise(sol))
+
+ # replace leads successively with self-energy
+ class LeadWithOnlySelfEnergy(object):
+ def __init__(self, lead):
+ self.lead = lead
+
+ def self_energy(self, energy):
+ return self.lead.self_energy(energy)
+
+ fsys.leads[0] = LeadWithOnlySelfEnergy(fsys.leads[0])
+ sol = kwant.solve(fsys)
+ assert_almost_equal(noise_should_be, two_terminal_shotnoise(sol))
+
+ fsys.leads[1] = LeadWithOnlySelfEnergy(fsys.leads[1])
+ sol = kwant.solve(fsys)
+ assert_almost_equal(noise_should_be, two_terminal_shotnoise(sol))
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
new file mode 100644
index 0000000..f1d8395
--- /dev/null
+++ b/kwant/solvers/common.py
@@ -0,0 +1,472 @@
+"""Collection of things commonly used by the different solvers.
+
+Typically needs not be called by the user, but is rather used by the
+individual solver modules
+"""
+
+__all__ = ['SparseSolver', 'BlockResult']
+
+from collections import namedtuple
+import abc
+import numpy as np
+import scipy.sparse as sp
+from kwant import physics, system
+
+# Currently, scipy.sparse does not support matrices with one dimension being
+# zero: http://projects.scipy.org/scipy/ticket/1602 We use numpy dense matrices
+# as a replacement.
+# TODO: Once this issue is fixed, code for the special cases can be removed
+# from make_linear_sys, solve_linear_sys and possibly other places marked by
+# the line "See comment about zero-shaped sparse matrices at the top of
+# _sparse.py".
+
+
+LinearSys = namedtuple('LinearSys', ['lhs', 'rhs', 'kept_vars'])
+
+class SparseSolver(object):
+ """Solver class for computing physical quantities based on solving
+ a liner system of equations.
+
+ `SparseSolver` provides the following functionality:
+
+ - `solve`: Compute a scattering matrix or Green's function
+ - `ldos`: Compute the local density of states
+
+ `SparseSolver` is an abstract base class. In particular, a derived class
+ must implement a function that does the actual linear solve step. To be
+ more precise, a derived class must implement the method
+
+ - `solve_linear_sys`, that solves a linear system of equations
+
+ and the properties
+
+ - `lhsformat`, `rhsformat` that give the sparse matrix format of
+ left and right hand sides of the linear system, repsectively;
+ - `nrhs` which determines how many right hand side vectors should be
+ solved at once for efficiency.
+
+ For details see the documentation of the respective abstract methods and
+ properties.
+ """
+ __metaclass__ = abc.ABCMeta
+
+ @abc.abstractmethod
+ def solve_linear_sys(self, a, b, kept_vars, factored=None):
+ """
+ Solve the linar system `a x = b`, returning the part of
+ the result indicated in kept_vars.
+
+ A particular implementation of a sparse solver must
+ provide an implementation of this abstract method.
+ """
+ pass
+
+ @abc.abstractproperty
+ def lhsformat(self):
+ """Sparse matrix format of the left hand side in the sparse linear
+ system. Must be one of {'coo', 'csc', 'csr'}.
+ """
+ pass
+
+ @abc.abstractproperty
+ def rhsformat(self):
+ """Sparse matrix format of the right hand side in the sparse linear
+ system. Must be one of {'coo', 'csc', 'csr'}.
+ """
+ pass
+
+ @abc.abstractproperty
+ def nrhs(self):
+ """Number of rhs vectors that should be solved at once in a call
+ to `solve_linear_sys`, when the full solution is computed (i.e.
+ kept_vars covers all entries in the solution). This should be not too
+ big too avoid excessive memory usage, but for some solvers not too
+ small for performance reasons. Only used for `ldos`.
+ """
+ pass
+
+
+ def make_linear_sys(self, sys, out_leads, in_leads, energy=0,
+ force_realspace=False, check_hermiticity=True):
+ """
+ Make a sparse linear system of equations defining a scattering
+ problem.
+
+ Parameters
+ ----------
+ sys : kwant.system.FiniteSystem
+ low level system, containing the leads and the Hamiltonian of a
+ scattering region.
+ out_leads : list of integers
+ numbers of leads where current or wave function is extracted
+ in_leads : list of integers
+ numbers of leads in which current or wave function is injected.
+ energy : number
+ excitation energy at which to solve the scattering problem.
+ force_realspace : bool
+ calculate Green's function between the outermost lead
+ sites, instead of lead modes. This is almost always
+ more computationally expensive and less stable.
+ check_hermiticity : bool
+ check if Hamiltonian matrices are in fact Hermitian.
+
+ Returns
+ -------
+ (lhs, rhs, kept_vars) : LinearSys
+ `lhs` is a numpy.sparse.csc_matrix, containing the left hand
+ side of the system of equations. `rhs` is a list of matrices
+ with the right hand side, with each matrix corresponding to one
+ lead mentioned in `in_leads`. `kept_vars` is a list of numbers
+ of variables in the solution that have to be stored (typically
+ a small part of the 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
+ -----
+ Both incoming and outgoing leads can be defined via either self-energy,
+ or a low-level translationally invariant system.
+ The system of equations that is created is described in
+ kwant/doc/other/linear_system.pdf
+ """
+
+ splhsmat = getattr(sp, self.lhsformat + '_matrix')
+ sprhsmat = getattr(sp, self.rhsformat + '_matrix')
+
+ if not sys.lead_neighbor_seqs:
+ raise ValueError('System contains no leads.')
+ lhs, norb = sys.hamiltonian_submatrix(sparse=True)[:2]
+ lhs = getattr(lhs, 'to' + self.lhsformat)()
+ lhs = lhs - energy * sp.identity(lhs.shape[0], format=self.lhsformat)
+
+ if check_hermiticity:
+ if np.any(abs((lhs - lhs.T.conj()).data) > 1e-13):
+ raise ValueError('System Hamiltonian is not Hermitian.')
+
+ offsets = np.zeros(norb.shape[0] + 1, int)
+ offsets[1 :] = np.cumsum(norb)
+
+ # Process the leads, generate the eigenvector matrices and lambda
+ # vectors. Then create blocks of the linear system and add them
+ # step by step.
+ kept_vars = []
+ rhs = []
+ lead_info = []
+ for leadnum, lead_neighbors in enumerate(sys.lead_neighbor_seqs):
+ lead = sys.leads[leadnum]
+ if isinstance(lead, system.InfiniteSystem) and not force_realspace:
+ h = lead.slice_hamiltonian()
+
+ if check_hermiticity:
+ if not np.allclose(h, h.T.conj(), rtol=1e-13):
+ msg = "Lead number {0} has a non-Hermitian " \
+ "slice Hamiltonian."
+ raise ValueError(msg.format(leadnum))
+
+ h -= energy * np.identity(h.shape[0])
+ v = lead.inter_slice_hopping()
+ modes = physics.modes(h, v)
+ lead_info.append(modes)
+ if not np.any(v):
+ # See comment about zero-shaped sparse matrices at the top.
+ rhs.append(np.zeros((lhs.shape[1], 0)))
+ continue
+ u, ulinv, nprop, svd = modes
+
+ if leadnum in out_leads:
+ kept_vars.extend(range(lhs.shape[0], lhs.shape[0] + nprop))
+
+ u_out, ulinv_out = u[:, nprop:], ulinv[:, nprop:]
+ u_in, ulinv_in = u[:, : nprop], ulinv[:, : nprop]
+
+ # Construct a matrix of 1's that translates the
+ # inter-slice hopping to a proper hopping
+ # from the system to the lead.
+ neighbors = np.r_[tuple(np.arange(offsets[i], offsets[i + 1])
+ for i in lead_neighbors)]
+ coords = np.r_[[np.arange(neighbors.size)], [neighbors]]
+ tmp = sp.csc_matrix((np.ones(neighbors.size), coords),
+ shape=(neighbors.size, lhs.shape[0]))
+
+ if svd is not None:
+ v_sp = sp.csc_matrix(svd[2].T.conj()) * tmp
+ vdaguout_sp = tmp.T * sp.csc_matrix(np.dot(svd[2] * svd[1],
+ u_out))
+ lead_mat = - ulinv_out
+ else:
+ v_sp = tmp
+ vdaguout_sp = tmp.T * sp.csc_matrix(np.dot(v.T.conj(),
+ u_out))
+ lead_mat = - ulinv_out
+
+ lhs = sp.bmat([[lhs, vdaguout_sp], [v_sp, lead_mat]],
+ format=self.lhsformat)
+
+ if leadnum in in_leads and nprop > 0:
+ if svd:
+ vdaguin_sp = tmp.T * sp.csc_matrix(
+ -np.dot(svd[2] * svd[1], u_in))
+ else:
+ vdaguin_sp = tmp.T * sp.csc_matrix(
+ -np.dot(v.T.conj(), u_in))
+
+ # defer formation of the real matrix until the proper
+ # system size is known
+ rhs.append((vdaguin_sp, ulinv_in))
+ else:
+ # See comment about zero-shaped sparse matrices at the top.
+ rhs.append(np.zeros((lhs.shape[1], 0)))
+ else:
+ sigma = lead.self_energy(energy)
+ lead_info.append(sigma)
+ indices = np.r_[tuple(range(offsets[i], offsets[i + 1]) for i
+ in lead_neighbors)]
+ assert sigma.shape == 2 * indices.shape
+ y, x = np.meshgrid(indices, indices)
+ sig_sparse = splhsmat((sigma.flat, [x.flat, y.flat]),
+ lhs.shape)
+ lhs = lhs + sig_sparse # __iadd__ is not implemented in v0.7
+ if leadnum in out_leads:
+ kept_vars.extend(list(indices))
+ if leadnum in in_leads:
+ # defer formation of true rhs until the proper system
+ # size is known
+ rhs.append((indices, ))
+
+
+ # Resize the right-hand sides to be compatible with the full lhs
+ for i, mats in enumerate(rhs):
+ if isinstance(mats, tuple):
+ if len(mats) == 1:
+ # self-energy lead
+ l = mats[0].shape[0]
+ rhs[i] = sprhsmat((-np.ones(l), [mats[0], np.arange(l)]),
+ shape=(lhs.shape[0], l))
+ else:
+ # lead with modes
+ zero_rows = (lhs.shape[0] - mats[0].shape[0] -
+ mats[1].shape[0])
+
+ if zero_rows:
+ zero_mat = sprhsmat((zero_rows, mats[0].shape[1]))
+ bmat = [[mats[0]], [mats[1]], [zero_mat]]
+ else:
+ bmat = [[mats[0]], [mats[1]]]
+
+ rhs[i] = sp.bmat(bmat, format=self.rhsformat)
+
+ return LinearSys(lhs, rhs, kept_vars), lead_info
+
+
+ def solve(self, sys, energy=0, out_leads=None, in_leads=None,
+ force_realspace=False, check_hermiticity=True):
+ """
+ Calculate a Green's function of a system.
+
+ Parameters
+ ----------
+ sys : `kwant.system.FiniteSystem`
+ low level system, containing the leads and the Hamiltonian of a
+ scattering region.
+ energy : number
+ excitation energy at which to solve the scattering problem.
+ out_leads : list of integers or None
+ numbers of leads where current or wave function is extracted.
+ None is interpreted as all leads. Default is None.
+ in_leads : list of integers or None
+ numbers of leads in which current or wave function is injected.
+ None is interpreted as all leads. Default is None.
+ force_realspace : bool
+ calculate Green's function between the outermost lead
+ sites, instead of lead modes. This is almost always
+ more computationally expensive and less stable.
+ check_hermiticity : bool
+ check if the Hamiltonian matrices are Hermitian.
+
+ Returns
+ -------
+ output : `BlockResult`
+ see notes below and `BlockResult` docstring for more
+ information about the output format.
+
+ Notes
+ -----
+ Both in_leads and out_leads must be sorted and must only contain
+ unique entries.
+
+ Returns the Green's function elements between in_leads and
+ out_leads. If the leads are defined as a self-energy, the result is
+ just the real space retarded Green's function between from in_leads
+ to out_leads. If the leads are defined as tight-binding systems,
+ then Green's function from incoming to outgoing modes is
+ returned. Also returned is a list containing the output of
+ `kwant.physics.modes` for the leads which are defined as builders,
+ and self-energies for leads defined via self-energy. This list
+ allows to split the Green's function into blocks corresponding to
+ different leads. The Green's function elements between incoming and
+ outgoing modes form the scattering matrix of the system. If some
+ leads are defined via self-energy, and some as tight-binding
+ systems, result has Green's function's elements between modes and
+ sites.
+
+ Alternatively, if force_realspace=True is used, G^R is returned
+ always in real space, however this option is more computationally
+ expensive and can be less stable.
+ """
+
+ n = len(sys.lead_neighbor_seqs)
+ if in_leads is None:
+ in_leads = range(n)
+ if out_leads is None:
+ out_leads = range(n)
+ if sorted(in_leads) != in_leads or sorted(out_leads) != out_leads or \
+ len(set(in_leads)) != len(in_leads) or \
+ len(set(out_leads)) != len(out_leads):
+ 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, lead_info = self.make_linear_sys(sys, out_leads, in_leads,
+ energy, force_realspace,
+ check_hermiticity)
+
+ result = BlockResult(self.solve_linear_sys(*linsys)[0], lead_info)
+
+ result.in_leads = in_leads
+ result.out_leads = out_leads
+
+ return result
+
+
+ def ldos(self, fsys, energy=0):
+ """
+ Calculate the local density of states of a system at a given energy.
+
+ Parameters
+ ----------
+ sys : `kwant.system.FiniteSystem`
+ low level system, containing the leads and the Hamiltonian of the
+ scattering region.
+ energy : number
+ excitation energy at which to solve the scattering problem.
+
+ Returns
+ -------
+ ldos : a numpy array
+ local density of states at each orbital of the system.
+ """
+ for lead in fsys.leads:
+ if not isinstance(lead, system.InfiniteSystem):
+ # TODO: fix this
+ raise ValueError("ldos only works when all leads are "
+ "tight binding systems.")
+
+ (h, rhs, kept_vars), lead_info = \
+ self.make_linear_sys(fsys, [], xrange(len(fsys.leads)), energy)
+
+ Modes = physics.Modes
+ 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
+
+ ldos = np.zeros(num_orb, float)
+ factored = None
+
+ for mat in rhs:
+ if mat.shape[1] == 0:
+ continue
+
+ for j in xrange(0, mat.shape[1], self.nrhs):
+ jend = min(j + self.nrhs, mat.shape[1])
+
+ psi, factored = self.solve_linear_sys(h, [mat[:, j:jend]],
+ slice(num_orb),
+ factored)
+
+ ldos += np.sum(np.square(abs(psi)), axis=1)
+
+ return ldos * (0.5 / np.pi)
+
+
+class BlockResult(namedtuple('BlockResultTuple', ['data', 'lead_info'])):
+ """
+ Solution of a transport problem, subblock of retarded Green's function.
+
+ This class is derived from ``namedtuple('BlockResultTuple', ['data',
+ 'lead_info'])``. In addition to direct access to `data` and `lead_info`,
+ this class also supports a higher level interface via its methods.
+
+ Instance Variables
+ ------------------
+ data : numpy matrix
+ a matrix containing all the requested matrix elements of Green's
+ function.
+ lead_info : list of data
+ a list with output of `kwant.physics.modes` for each lead defined as a
+ builder, and self-energy for each lead defined as self-energy term.
+ """
+ def block_coords(self, lead_out, lead_in):
+ """
+ Return slices corresponding to the block from lead_in to lead_out.
+ """
+ lead_out = self.out_leads.index(lead_out)
+ lead_in = self.in_leads.index(lead_in)
+ if not hasattr(self, '_sizes'):
+ sizes = []
+ for i in self.lead_info:
+ if isinstance(i, tuple):
+ sizes.append(i[2])
+ else:
+ sizes.append(i.shape[0])
+ self._sizes = np.array(sizes)
+ self._in_offsets = np.zeros(len(self.in_leads) + 1, int)
+ self._in_offsets[1 :] = np.cumsum(self._sizes[self.in_leads])
+ self._out_offsets = np.zeros(len(self.out_leads) + 1, int)
+ self._out_offsets[1 :] = np.cumsum(self._sizes[self.out_leads])
+ return slice(self._out_offsets[lead_out],
+ self._out_offsets[lead_out + 1]), \
+ slice(self._in_offsets[lead_in], self._in_offsets[lead_in + 1])
+
+ def submatrix(self, lead_out, lead_in):
+ """Return the matrix elements from lead_in to lead_out."""
+ return self.data[self.block_coords(lead_out, lead_in)]
+
+ def _a_ttdagger_a_inv(self, lead_out, lead_in):
+ gf = self.submatrix(lead_out, lead_in)
+ factors = []
+ for lead, gf2 in ((lead_out, gf), (lead_in, gf.conj().T)):
+ possible_se = self.lead_info[lead]
+ if not isinstance(possible_se, tuple):
+ # Lead is a "self energy lead": multiply gf2 with a gamma
+ # matrix.
+ factors.append(1j * (possible_se - possible_se.conj().T))
+ factors.append(gf2)
+ return reduce(np.dot, factors)
+
+ def transmission(self, lead_out, lead_in):
+ """Return transmission from lead_in to lead_out."""
+ if isinstance(self.lead_info[lead_out], tuple) and \
+ isinstance(self.lead_info[lead_in], tuple):
+ return np.linalg.norm(self.submatrix(lead_out, lead_in))**2
+ else:
+ result = np.trace(self._a_ttdagger_a_inv(lead_out, lead_in)).real
+ if lead_out == lead_in:
+ # For reflection we have to be more careful
+ gamma = 1j * (self.lead_info[lead_in] -
+ self.lead_info[lead_in].conj().T)
+ gf = self.submatrix(lead_out, lead_in)
+
+ # The number of channels is given by the number of
+ # nonzero eigenvalues of Gamma
+ # rationale behind the threshold from
+ # Golub; van Loan, chapter 5.5.8
+ eps = np.finfo(gamma.dtype).eps * 1000
+ N = np.sum(np.linalg.eigvalsh(gamma) >
+ eps * np.linalg.norm(gamma, np.inf))
+
+ result += 2 * np.trace(np.dot(gamma, gf)).imag + N
+
+ return result
diff --git a/kwant/solvers/mumps.py b/kwant/solvers/mumps.py
new file mode 100644
index 0000000..2524771
--- /dev/null
+++ b/kwant/solvers/mumps.py
@@ -0,0 +1,163 @@
+"""Implementation of the sparse solver framework using the direct sparse solver
+MUMPS (http://graal.ens-lyon.fr/MUMPS/, only the sequential, single core
+version is used).
+
+MUMPS is a very efficient direct sparse solver that can take advantage of
+memory beyond 3GB for the solution of large problems. Furthermore, it offers a
+choice of several orderings of the input matrix from which can speed up a
+calculation significantly.
+
+Compared to the generic sparse solver framework, `mumps` adds the following
+control options that may affect performance:
+
+- `ordering`: a fill-in reducing ordering of the matrix
+- `nrhs`: number of right hand sides that should be solved simultaneously
+- `sparse_rhs`: whether to use dense or sparse right hand sides
+
+For more details see `~Solver.options`.
+"""
+
+__all__ = ['solve', 'ldos', 'options', 'Solver']
+
+import numpy as np
+import scipy.sparse as sp
+import kwant.solvers.common
+import kwant.linalg.mumps as mumps
+
+class Solver(kwant.solvers.common.SparseSolver):
+ """Sparse Solver class based on the sparse direct solver MUMPS.
+ """
+
+ lhsformat = 'coo'
+ rhsformat = 'csc'
+ nrhs = 1 # later governed by options(), reset_options()
+
+ def __init__(self):
+ self.nrhs = self.ordering = self.sparse_rhs = None
+ self.reset_options()
+
+ def reset_options(self):
+ """Set the options to default values."""
+ self.options(nrhs=6, ordering='kwant_decides', sparse_rhs=False)
+
+ def options(self, nrhs=None, ordering=None, sparse_rhs=None):
+ """
+ Parameters
+ ----------
+ nrhs : number
+ number of right hand sides that should be solved simultaneously. A
+ value around 5-10 gives optimal performance on many machines. If
+ memory is an issue, it can be set to 1, to minimize memory usage
+ (at the cost of slower performance). Default value is 6.
+ ordering : string
+ one of the ordering methods supported by the MUMPS solver (see
+ `~kwant.linalg.mumps`. The availability of certain orderings
+ depends on the MUMPS installation.), or 'kwant_decides'. If
+ ``ordering=='kwant_decides'``, the ordering that typically gives
+ the best performance is chosen from the available ones. One can
+ also defer the choice of ordering to MUMPS by specifying 'auto', in
+ some cases MUMPS however chooses poorly.
+
+ The choice of ordering can significantly influence the performance
+ and memory impact of the solve phase. Typically the nested
+ dissection orderings 'metis' and 'scotch' are most suited for
+ physical systems. Default is 'kwant_decides'
+ sparse_rhs : True or False
+ whether to use a sparse right hand side in the solve phase of
+ MUMPS. Preliminary tests have not shown a significant performance
+ increase when this feature is used, but this needs more looking
+ into. Default value is False.
+
+ Returns
+ -------
+ old_options: dict
+ dictionary containing the previous options.
+ """
+
+ old_opts = {'nrhs': self.nrhs,
+ 'ordering': self.ordering,
+ 'sparse_rhs': self.sparse_rhs}
+
+ if nrhs is not None:
+ if nrhs < 1 and int(nrhs) != nrhs:
+ raise ValueError("nrhs must be an integer bigger than zero")
+ nrhs = int(nrhs)
+ self.nrhs = nrhs
+
+ if ordering is not None:
+ if ordering not in mumps.orderings.keys() + ['kwant_decides']:
+ raise ValueError("Invalid ordering: " + ordering)
+ if ordering == 'kwant_decides':
+ # Choose what is considered to be the best ordering.
+ sorted_orderings = [ordering
+ for ordering in ['metis', 'scotch', 'auto']
+ if ordering in mumps.possible_orderings()]
+ ordering = sorted_orderings[0]
+ self.ordering = ordering
+
+ if sparse_rhs is not None:
+ self.sparse_rhs = bool(sparse_rhs)
+
+ return old_opts
+
+ def solve_linear_sys(self, a, b, kept_vars=None, factored=None):
+
+ """
+ Solve matrix system of equations a x = b with sparse input,
+ using MUMPS.
+
+ Parameters
+ ----------
+ a : a scipy.sparse.coo_matrix sparse matrix.
+ b : a list of scipy.sparse.csc_matrices.
+ kept_vars : list of integers
+ a list of numbers of variables to keep in the solution.
+ factored : a factorized lhs as returned by a previous call
+ to solve_linear_sys with the same lhs.
+
+ Returns
+ -------
+ output : a numpy matrix
+ solution to the system of equations.
+ factored : factorized lhs. Can be reused in later solves with
+ the same lhs, but different rhs.
+ """
+
+ if kept_vars == None:
+ kept_vars = [range(a.shape[1])]
+
+ sols = []
+
+ if not factored:
+ inst = mumps.MUMPSContext()
+ inst.factor(a, ordering=self.ordering)
+ else:
+ inst = factored
+
+ for mat in b:
+ if mat.shape[1] != 0:
+ # See comment about zero-shaped sparse matrices at the top
+ # of sparse.
+ mat = sp.csr_matrix(mat)
+
+ for j in xrange(0, mat.shape[1], self.nrhs):
+ jend = min(j + self.nrhs, mat.shape[1])
+
+ if self.sparse_rhs:
+ sols.append(inst.solve(mat[:, j:jend])[kept_vars, :])
+ else:
+ sols.append(inst.solve(mat[:, j:jend].todense())
+ [kept_vars, :])
+
+ if len(sols):
+ return np.concatenate(sols, axis=1), inst
+ else:
+ return np.zeros(shape=(len(kept_vars), 0)), inst
+
+
+default_solver = Solver()
+
+solve = default_solver.solve
+ldos = default_solver.ldos
+options = default_solver.options
+reset_options = default_solver.reset_options
diff --git a/kwant/solvers/sparse.py b/kwant/solvers/sparse.py
index e0f6304..566b203 100644
--- a/kwant/solvers/sparse.py
+++ b/kwant/solvers/sparse.py
@@ -1,442 +1,164 @@
-__all__ = ['make_linear_sys', 'LinearSys', 'solve', 'ldos', 'BlockResult']
+"""Implementation of the sparse solver framework using the direct sparse solver
+provided by `scipy.sparse.linalg
+<http://docs.scipy.org/doc/scipy/reference/sparse.linalg.html>`.
+
+Scipy currently uses internally either the direct sparse solver Umfpack or if
+that is not installed, SuperLU. Often, scipy's SuperLU will give quite poor
+performance and you will be warned if only SuperLU is found. The module
+variable `uses_umfpack` can be checked to determine if Umfpack is being used.
+
+`sparse` does not introduce any additional options as compared to the generic
+sparse solver framework.
+"""
+
+__all__ = ['solve', 'ldos', 'Solver']
from functools import reduce
-from collections import namedtuple
+import warnings
+import kwant.solvers.common
import numpy as np
-import scipy.linalg as la
import scipy.sparse as sp
-import scipy.sparse.linalg.dsolve.umfpack as umfpack
-from kwant import physics, system
-
-
-# Currently, scipy.sparse does not support matrices with one dimension being
-# zero: http://projects.scipy.org/scipy/ticket/1602 We use numpy dense matrices
-# as a replacement.
-# TODO: Once this issue is fixed, code for the special cases can be removed
-# from make_linear_sys, solve_linear_sys and possibly other places marked by
-# the line "See comment about zero-shaped sparse matrices at the top."
-
-
-# This patches a memory leak in scipy:
-# http://projects.scipy.org/scipy/ticket/1597
-#
-# TODO: Remove this code once it is likely that the official bug fix has
-# reached all of our users.
-def del_for_umfpackcontext(self):
- self.free()
-if not hasattr(umfpack.UmfpackContext, '__del__'):
- umfpack.UmfpackContext.__del__ = del_for_umfpackcontext
-del del_for_umfpackcontext
-
-def factorized(A, piv_tol=1.0, sym_piv_tol=1.0):
- """
- Return a fuction for solving a sparse linear system, with A pre-factorized.
-
- Example:
- solve = factorized( A ) # Makes LU decomposition.
- x1 = solve( rhs1 ) # Uses the LU factors.
- x2 = solve( rhs2 ) # Uses again the LU factors.
-
- Parameters
- ----------
- A : csc_matrix
- matrix to be factorized
- piv_tol : float, 0 <= piv_tol <= 1.0
- sym_piv_tol : float, 0 <= piv_tol <= 1.0
- thresholds used by umfpack for pivoting. 0 means no pivoting,
- 1.0 means full pivoting as in dense matrices (guaranteeing
- stability, but reducing possibly sparsity). Defaults of umfpack
- are 0.1 and 0.001 respectively. Whether piv_tol or sym_piv_tol
- are used is decided internally by umfpack, depending on whether
- the matrix is "symmetric" enough.
- """
- if not sp.isspmatrix_csc(A):
- A = sp.csc_matrix(A)
+# Note: previous code would have failed if umfpack was provided by scikit
+import scipy.sparse.linalg.dsolve.linsolve as linsolve
+umfpack = linsolve.umfpack
+uses_umfpack = linsolve.isUmfpack
+
+import kwant.system as system
+import kwant.physics as physics
+
+# check if we are actually using umfpack or rather SuperLU
+
+if uses_umfpack:
+ # This patches a memory leak in scipy:
+ # http://projects.scipy.org/scipy/ticket/1597
+ #
+ # TODO: Remove this code once it is likely that the official bug fix has
+ # reached all of our users.
+ def del_for_umfpackcontext(self):
+ self.free()
+ if not hasattr(umfpack.UmfpackContext, '__del__'):
+ umfpack.UmfpackContext.__del__ = del_for_umfpackcontext
+ del del_for_umfpackcontext
+
+ def factorized(A, piv_tol=1.0, sym_piv_tol=1.0):
+ """
+ Return a fuction for solving a sparse linear system, with A
+ pre-factorized.
+
+ Example:
+ solve = factorized( A ) # Makes LU decomposition.
+ x1 = solve( rhs1 ) # Uses the LU factors.
+ x2 = solve( rhs2 ) # Uses again the LU factors.
+
+ Parameters
+ ----------
+ A : csc_matrix
+ matrix to be factorized
+ piv_tol : float, 0 <= piv_tol <= 1.0
+ sym_piv_tol : float, 0 <= piv_tol <= 1.0
+ thresholds used by umfpack for pivoting. 0 means no pivoting, 1.0
+ means full pivoting as in dense matrices (guaranteeing stability,
+ but reducing possibly sparsity). Defaults of umfpack are 0.1 and
+ 0.001 respectively. Whether piv_tol or sym_piv_tol are used is
+ decided internally by umfpack, depending on whether the matrix is
+ "symmetric" enough.
+ """
- A.sort_indices()
- A = A.asfptype() #upcast to a floating point format
+ if not sp.isspmatrix_csc(A):
+ A = sp.csc_matrix(A)
- if A.dtype.char not in 'dD':
- raise ValueError("convert matrix data to double, please, using"
- " .astype()")
+ A.sort_indices()
+ A = A.asfptype() #upcast to a floating point format
- family = {'d' : 'di', 'D' : 'zi'}
- umf = umfpack.UmfpackContext( family[A.dtype.char] )
+ if A.dtype.char not in 'dD':
+ raise ValueError("convert matrix data to double, please, using"
+ " .astype()")
- # adjust pivot thresholds
- umf.control[umfpack.UMFPACK_PIVOT_TOLERANCE] = piv_tol
- umf.control[umfpack.UMFPACK_SYM_PIVOT_TOLERANCE] = sym_piv_tol
+ family = {'d' : 'di', 'D' : 'zi'}
+ umf = umfpack.UmfpackContext( family[A.dtype.char] )
- # Make LU decomposition.
- umf.numeric( A )
+ # adjust pivot thresholds
+ umf.control[umfpack.UMFPACK_PIVOT_TOLERANCE] = piv_tol
+ umf.control[umfpack.UMFPACK_SYM_PIVOT_TOLERANCE] = sym_piv_tol
- def solve( b ):
- return umf.solve( umfpack.UMFPACK_A, A, b, autoTranspose = True )
+ # Make LU decomposition.
+ umf.numeric( A )
- return solve
+ def solve( b ):
+ return umf.solve( umfpack.UMFPACK_A, A, b, autoTranspose = True )
+ return solve
+else:
+ # no Umfpack found. SuperLU is being used, but usually abysmally slow
+ # (SuperLu is not bad per se, somehow the scipy version isn't good)
+ warnings.warn("The installed scipy does not use Umfpack. Instead, "
+ "scipy will use the version of SuperLu it is shipped with. "
+ "Performance can be very poor in this case.", RuntimeWarning)
-LinearSys = namedtuple('LinearSys', ['lhs', 'rhs', 'kept_vars'])
+ factorized = linsolve.factorized
-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.
-
- Parameters
- ----------
- sys : kwant.system.FiniteSystem
- low level system, containing the leads and the Hamiltonian of a
- scattering region.
- out_leads : list of integers
- numbers of leads where current or wave function is extracted
- in_leads : list of integers
- numbers of leads in which current or wave function is injected.
- energy : number
- excitation energy at which to solve the scattering problem.
- force_realspace : bool
- calculate Green's function between the outermost lead
- sites, instead of lead modes. This is almost always
- more computationally expensive and less stable.
-
- Returns
- -------
- (lhs, rhs, kept_vars) : LinearSys
- `lhs` is a numpy.sparse.csc_matrix, containing the left hand side of
- the system of equations. `rhs` is a list of matrices with the right
- hand side, with each matrix corresponding to one lead mentioned in
- `in_leads`. `kept_vars` is a list of numbers of variables in the
- solution that have to be stored (typically a small part of the 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
- -----
- Both incoming and outgoing leads can be defined via either self-energy,
- or a low-level translationally invariant system.
- The system of equations that is created is described in
- kwant/doc/other/linear_system.pdf
- """
- if not sys.lead_neighbor_seqs:
- raise ValueError('System contains no leads.')
- lhs, norb = sys.hamiltonian_submatrix(sparse=True)[:2]
- lhs = lhs.tocsc()
- lhs = lhs - energy * sp.identity(lhs.shape[0], format='csc')
-
- # Hermiticity check.
- if np.any(np.abs((lhs - lhs.T.conj()).data) > 1e-13):
- raise ValueError('System Hamiltonian is not Hermitian.')
-
- offsets = np.zeros(norb.shape[0] + 1, int)
- offsets[1 :] = np.cumsum(norb)
-
- # Process the leads, generate the eigenvector matrices and lambda vectors.
- # Then create blocks of the linear system and add them step by step.
- kept_vars = []
- rhs = []
- lead_info = []
- for leadnum, lead_neighbors in enumerate(sys.lead_neighbor_seqs):
- lead = sys.leads[leadnum]
- if isinstance(lead, system.InfiniteSystem) and not force_realspace:
- h = lead.slice_hamiltonian()
-
- # Hermiticity check.
- if not np.allclose(h, h.T.conj(), rtol=1e-13):
- msg = 'Lead number {0} has a non-Hermitian slice Hamiltonian.'
- raise ValueError(msg.format(leadnum))
-
- h -= energy * np.identity(h.shape[0])
- v = lead.inter_slice_hopping()
- modes = physics.modes(h, v)
- lead_info.append(modes)
- if not np.any(v):
- # See comment about zero-shaped sparse matrices at the top.
- rhs.append(np.zeros((lhs.shape[1], 0)))
- continue
- u, ulinv, nprop, svd = modes
-
- if leadnum in out_leads:
- kept_vars.append(range(lhs.shape[0], lhs.shape[0] + nprop))
-
- u_out, ulinv_out = u[:, nprop:], ulinv[:, nprop:]
- u_in, ulinv_in = u[:, : nprop], ulinv[:, : nprop]
-
- # Construct a matrix of 1's that translates the
- # inter-slice hopping to a proper hopping
- # from the system to the lead.
- neighbors = np.r_[tuple(np.arange(offsets[i], offsets[i + 1])
- for i in lead_neighbors)]
- coords = np.r_[[np.arange(neighbors.size)], [neighbors]]
- tmp = sp.csc_matrix((np.ones(neighbors.size), coords),
- (neighbors.size, lhs.shape[0]))
-
- if svd is not None:
- v_sp = sp.csc_matrix(svd[2].T.conj()) * tmp
- vdaguout_sp = tmp.T * sp.csr_matrix(np.dot(svd[2] * svd[1],
- u_out))
- lead_mat = - ulinv_out
- else:
- v_sp = tmp
- vdaguout_sp = tmp.T * sp.csr_matrix(np.dot(v.T.conj(), u_out))
- lead_mat = - ulinv_out
-
- lhs = sp.bmat([[lhs, vdaguout_sp], [v_sp, lead_mat]])
-
- if leadnum not in in_leads:
- continue
- if nprop > 0:
- if svd:
- vdaguin_sp = tmp.T * sp.csr_matrix(
- -np.dot(svd[2] * svd[1], u_in))
- else:
- vdaguin_sp = tmp.T * sp.csr_matrix(
- -np.dot(v.T.conj(), u_in))
- rhs.append(sp.bmat([[vdaguin_sp], [ulinv_in]]))
- else:
- # See comment about zero-shaped sparse matrices at the top.
- rhs.append(np.zeros((lhs.shape[1], 0)))
- else:
- sigma = lead.self_energy(energy)
- lead_info.append(sigma)
- indices = np.r_[tuple(range(offsets[i], offsets[i + 1]) for i in
- lead_neighbors)]
- assert sigma.shape == 2 * indices.shape
- y, x = np.meshgrid(indices, indices)
- sig_sparse = sp.coo_matrix((sigma.flat, [x.flat, y.flat]),
- lhs.shape)
- lhs = lhs + sig_sparse # __iadd__ is not implemented in v0.7
- if leadnum in out_leads:
- kept_vars.append(list(indices))
- if leadnum in in_leads:
- l = indices.shape[0]
- rhs.append(sp.coo_matrix((-np.ones(l), [indices,
- np.arange(l)])))
-
- return LinearSys(lhs, rhs, kept_vars), lead_info
-
-
-def solve_linear_sys(a, b, kept_vars=None):
- """
- Solve matrix system of equations a x = b with sparse input.
-
- Parameters
- ----------
- a : a scipy.sparse.csc_matrix sparse matrix
- b : a list of matrices.
- Sizes of these matrices may be smaller than needed, the missing
- entries at the end are padded with zeros.
- kept_vars : list of lists of integers
- a list of numbers of variables to keep in the solution
-
- Returns
- -------
- output : a numpy matrix
- solution to the system of equations.
-
- Notes
- -----
- This function is largely a wrapper to `factorized`.
- """
- a = sp.csc_matrix(a)
-
- if kept_vars == None:
- kept_vars = [range(a.shape[1])]
- slv = factorized(a)
- keeptot = sum(kept_vars, [])
- sols = []
- vec = np.empty(a.shape[0], complex)
- for mat in b:
- if mat.shape[1] != 0:
- # See comment about zero-shaped sparse matrices at the top.
- mat = sp.csr_matrix(mat)
- for j in xrange(mat.shape[1]):
- vec[: mat.shape[0]] = mat[:, j].todense().flatten()
- vec[mat.shape[0] :] = 0
- sols.append(slv(vec)[keeptot])
- return np.mat(sols).T
-
-
-def solve(sys, energy=0, out_leads=None, in_leads=None, force_realspace=False):
+class Solver(kwant.solvers.common.SparseSolver):
+ """Sparse Solver class based on the sparse direct solvers provided
+ by scipy.
"""
- Calculate a Green's function of a system.
-
- Parameters
- ----------
- sys : `kwant.system.FiniteSystem`
- low level system, containing the leads and the Hamiltonian of a
- scattering region.
- energy : number
- excitation energy at which to solve the scattering problem.
- out_leads : list of integers
- numbers of leads where current or wave function is extracted
- in_leads : list of integers
- numbers of leads in which current or wave function is injected.
- force_realspace : bool
- calculate Green's function between the outermost lead
- sites, instead of lead modes. This is almost always
- more computationally expensive and less stable.
-
- Returns
- -------
- output : `BlockResult`
- see notes below and `BlockResult` docstring for more information about
- the output format.
-
- Notes
- -----
- Both in_leads and out_leads must be sorted and must only contain
- unique entries.
-
- Returns the Green's function elements between in_leads and out_leads. If
- the leads are defined as a self-energy, the result is just the real
- space retarded Green's function between from in_leads to out_leads. If the
- leads are defined as tight-binding systems, then Green's function from
- incoming to outgoing modes is returned. Also returned is a list containing
- the output of `kwant.physics.modes` for the leads which are defined as
- builders, and self-energies for leads defined via self-energy. This list
- allows to split the Green's function into blocks corresponding to different
- leads. The Green's function elements between incoming and outgoing modes
- form the scattering matrix of the system. If some leads are defined via
- self-energy, and some as tight-binding systems, result has Green's
- function's elements between modes and sites.
-
- Alternatively, if force_realspace=True is used, G^R is returned
- always in real space, however this option is more computationally
- expensive and can be less stable.
- """
- n = len(sys.lead_neighbor_seqs)
- if in_leads is None:
- in_leads = range(n)
- if out_leads is None:
- out_leads = range(n)
- if sorted(in_leads) != in_leads or sorted(out_leads) != out_leads or \
- len(set(in_leads)) != len(in_leads) or \
- len(set(out_leads)) != len(out_leads):
- 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, lead_info = make_linear_sys(sys, out_leads, in_leads, energy,
- force_realspace)
- out_modes = [len(i) for i in linsys.kept_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
-
-
-class BlockResult(namedtuple('BlockResultTuple', ['data', 'lead_info'])):
- """
- Solution of a transport problem, subblock of retarded Green's function.
-
- This class is derived from ``namedtuple('BlockResultTuple', ['data',
- 'lead_info'])``. In addition to direct access to `data` and `lead_info`,
- this class also supports a higher level interface via its methods.
-
- Instance Variables
- ------------------
- data : numpy matrix
- a matrix containing all the requested matrix elements of Green's
- function.
- lead_info : list of data
- a list with output of `kwant.physics.modes` for each lead defined as a
- builder, and self-energy for each lead defined as self-energy term.
- """
- def block_coords(self, lead_out, lead_in):
+
+ lhsformat = 'csc'
+ rhsformat = 'csc'
+ nrhs = 1
+
+ def solve_linear_sys(self, a, b, kept_vars=None, factored=None):
"""
- Return slices corresponding to the block from lead_in to lead_out.
+ Solve matrix system of equations a x = b with sparse input,
+ using UMFPACK.
+
+ Parameters
+ ----------
+ a : a scipy.sparse.csc_matrix sparse matrix
+ b : a list of matrices.
+ Sizes of these matrices may be smaller than needed, the missing
+ entries at the end are padded with zeros.
+ kept_vars : slice object or list of integers
+ a list of numbers of variables to keep in the solution
+
+ Returns
+ -------
+ output : a numpy matrix
+ solution to the system of equations.
+
+ Notes
+ -----
+ This function is largely a wrapper to `factorized`.
"""
- lead_out = self.out_leads.index(lead_out)
- lead_in = self.in_leads.index(lead_in)
- if not hasattr(self, '_sizes'):
- sizes = []
- for i in self.lead_info:
- if isinstance(i, tuple):
- sizes.append(i[2])
- else:
- sizes.append(i.shape[0])
- self._sizes = np.array(sizes)
- self._in_offsets = np.zeros(len(self.in_leads) + 1, int)
- self._in_offsets[1 :] = np.cumsum(self._sizes[self.in_leads])
- self._out_offsets = np.zeros(len(self.out_leads) + 1, int)
- self._out_offsets[1 :] = np.cumsum(self._sizes[self.out_leads])
- return slice(self._out_offsets[lead_out],
- self._out_offsets[lead_out + 1]), \
- slice(self._in_offsets[lead_in], self._in_offsets[lead_in + 1])
-
- def submatrix(self, lead_out, lead_in):
- """Return the matrix elements from lead_in to lead_out."""
- return self.data[self.block_coords(lead_out, lead_in)]
-
- def _a_ttdagger_a_inv(self, lead_out, lead_in):
- gf = self.submatrix(lead_out, lead_in)
- factors = []
- for lead, gf2 in ((lead_out, gf), (lead_in, gf.conj().T)):
- possible_se = self.lead_info[lead]
- if not isinstance(possible_se, tuple):
- # Lead is a "self energy lead": multiply gf2 with a gamma
- # matrix.
- factors.append(1j * (possible_se - possible_se.conj().T))
- factors.append(gf2)
- return reduce(np.dot, factors)
-
- def transmission(self, lead_out, lead_in):
- """Return transmission from lead_in to lead_out."""
- if isinstance(self.lead_info[lead_out], tuple) and \
- isinstance(self.lead_info[lead_in], tuple):
- return la.norm(self.submatrix(lead_out, lead_in))**2
+ a = sp.csc_matrix(a)
+
+ if kept_vars == None:
+ kept_vars = [range(a.shape[1])]
+
+ if not factored:
+ slv = factorized(a)
+ else:
+ slv = factored
+
+ sols = []
+ vec = np.empty(a.shape[0], complex)
+ for mat in b:
+ if mat.shape[1] != 0:
+ # See comment about zero-shaped sparse matrices at the top of
+ # _sparse.py.
+ mat = sp.csr_matrix(mat)
+ for j in xrange(mat.shape[1]):
+ vec[:] = mat[:, j].todense().flatten()
+ sols.append(slv(vec)[kept_vars])
+
+ if len(sols):
+ return np.asarray(sols).transpose(), factored
else:
- return np.trace(self._a_ttdagger_a_inv(lead_out, lead_in)).real
+ return np.asarray(np.zeros(shape=(len(kept_vars), 0))), factored
- def noise(self, lead_out, lead_in):
- """Return shot noise from lead_in to lead_out."""
- ttdag = self._a_ttdagger_a_inv(lead_out, lead_in)
- ttdag -= np.dot(ttdag, ttdag)
- return np.trace(ttdag).real
-def ldos(fsys, energy=0):
- """
- Calculate the local density of states of a system at a given energy.
-
- Parameters
- ----------
- sys : `kwant.system.FiniteSystem`
- low level system, containing the leads and the Hamiltonian of the
- scattering region.
- energy : number
- excitation energy at which to solve the scattering problem.
-
- Returns
- -------
- ldos : a numpy array
- local density of states at each orbital of the system.
- """
- for lead in fsys.leads:
- if not isinstance(lead, system.InfiniteSystem):
- # TODO: fix this
- msg = 'ldos only works when all leads are tight binding systems.'
- raise ValueError(msg)
- (h, rhs, kept_vars), lead_info = \
- make_linear_sys(fsys, [], xrange(len(fsys.leads)), energy)
- Modes = physics.Modes
- 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
- ldos = np.zeros(num_orb, float)
- slv = factorized(h)
- vec = np.empty(h.shape[0], complex)
- for mat in rhs:
- if mat.shape[1] == 0:
- continue
- mat = sp.csr_matrix(mat)
- for j in xrange(mat.shape[1]):
- vec[: mat.shape[0]] = mat[:, j].todense().flatten()
- vec[mat.shape[0] :] = 0
- psi = slv(vec)[:num_orb]
- ldos += abs(psi)**2
- return ldos * (0.5 / np.pi)
+default_solver = Solver()
+
+solve = default_solver.solve
+ldos = default_solver.ldos
diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py
new file mode 100644
index 0000000..ed3aeb2
--- /dev/null
+++ b/kwant/solvers/tests/_test_sparse.py
@@ -0,0 +1,364 @@
+from __future__ import division
+import numpy as np
+from nose.tools import assert_raises
+from numpy.testing import assert_equal, assert_almost_equal
+import kwant
+
+n = 5
+chain = kwant.lattice.Chain()
+square = kwant.lattice.Square()
+
+# Test output sanity: that an error is raised if no output is requested,
+# and that solving for a subblock of a scattering matrix is the same as taking
+# a subblock of the full scattering matrix.
+def test_output(solve):
+ np.random.seed(3)
+ system = kwant.Builder()
+ left_lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
+ right_lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
+ for b, site in [(system, chain(0)), (system, chain(1)),
+ (left_lead, chain(0)), (right_lead, chain(0))]:
+ h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ h += h.conjugate().transpose()
+ b[site] = h
+ for b, hopp in [(system, (chain(0), chain(1))),
+ (left_lead, (chain(0), chain(1))),
+ (right_lead, (chain(0), chain(1)))]:
+ b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ system.attach_lead(left_lead)
+ system.attach_lead(right_lead)
+ fsys = system.finalized()
+
+ result1 = solve(fsys)
+ s, modes1 = result1
+ assert s.shape == 2 * (sum(i[2] for i in modes1),)
+ s1 = result1.submatrix(1, 0)
+ s2, modes2 = solve(fsys, 0, [1], [0])
+ assert s2.shape == (modes2[1][2], modes2[0][2])
+ assert_almost_equal(s1, s2)
+ assert_almost_equal(np.dot(s.conjugate().transpose(), s),
+ np.identity(s.shape[0]))
+ assert_raises(ValueError, solve, fsys, 0, [])
+ modes = solve(fsys)[1]
+ h = fsys.leads[0].slice_hamiltonian()
+ t = fsys.leads[0].inter_slice_hopping()
+ modes1 = kwant.physics.modes(h, t)
+ h = fsys.leads[1].slice_hamiltonian()
+ t = fsys.leads[1].inter_slice_hopping()
+ modes2 = kwant.physics.modes(h, t)
+ assert_almost_equal(modes1[0], modes[0][0])
+ assert_almost_equal(modes2[1], modes[1][1])
+
+
+# Test that a system with one lead has unitary scattering matrix.
+def test_one_lead(solve):
+ np.random.seed(3)
+ system = kwant.Builder()
+ lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
+ for b, site in [(system, chain(0)), (system, chain(1)),
+ (system, chain(2)), (lead, chain(0))]:
+ h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ h += h.conjugate().transpose()
+ b[site] = h
+ for b, hopp in [(system, (chain(0), chain(1))),
+ (system, (chain(1), chain(2))),
+ (lead, (chain(0), chain(1)))]:
+ b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ system.attach_lead(lead)
+ fsys = system.finalized()
+
+ s = solve(fsys)[0]
+ assert_almost_equal(np.dot(s.conjugate().transpose(), s),
+ np.identity(s.shape[0]))
+
+
+# Test that a system with one lead with no propagating modes has a
+# 0x0 S-matrix.
+def test_smatrix_shape(solve):
+ system = kwant.Builder()
+ lead0 = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
+ lead1 = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
+ for b, site in [(system, chain(0)), (system, chain(1)),
+ (system, chain(2))]:
+ b[site] = 2
+ lead0[chain(0)] = lambda site: lead0_val
+ lead1[chain(0)] = lambda site: lead1_val
+
+ for b, hopp in [(system, (chain(0), chain(1))),
+ (system, (chain(1), chain(2))),
+ (lead0, (chain(0), chain(1))),
+ (lead1, (chain(0), chain(1)))]:
+ b[hopp] = -1
+ system.attach_lead(lead0)
+ system.attach_lead(lead1)
+ fsys = system.finalized()
+
+ lead0_val = 4
+ lead1_val = 4
+ s = solve(fsys, energy=1.0, out_leads=[1], in_leads=[0])[0]
+ assert s.shape == (0, 0)
+
+ lead0_val = 2
+ lead1_val = 2
+ s = solve(fsys, energy=1.0, out_leads=[1], in_leads=[0])[0]
+ assert s.shape == (1, 1)
+
+ lead0_val = 4
+ lead1_val = 2
+ s = solve(fsys, energy=1.0, out_leads=[1], in_leads=[0])[0]
+ assert s.shape == (1, 0)
+
+ lead0_val = 2
+ lead1_val = 4
+ s = solve(fsys, energy=1.0, out_leads=[1], in_leads=[0])[0]
+ assert s.shape == (0, 1)
+
+# Test that a translationally invariant system with two leads has only
+# transmission and that transmission does not mix modes.
+def test_two_equal_leads(solve):
+ def check_fsys():
+ sol = solve(fsys)
+ s, leads = sol[:2]
+ assert_almost_equal(np.dot(s.conjugate().transpose(), s),
+ np.identity(s.shape[0]))
+ n_modes = leads[0][2]
+ assert leads[1][2] == n_modes
+ assert_almost_equal(s[: n_modes, : n_modes], 0)
+ t_elements = np.sort(abs(np.asarray(s[n_modes :, : n_modes])),
+ axis=None)
+ t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
+ assert_almost_equal(t_elements, t_el_should_be)
+ assert_almost_equal(sol.transmission(1,0), n_modes)
+ np.random.seed(11)
+ system = kwant.Builder()
+ lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
+ h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ h += h.conjugate().transpose()
+ h *= 0.8
+ t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+ lead[chain(0)] = system[chain(0)] = h
+ lead[chain(0), chain(1)] = t
+ system.attach_lead(lead)
+ system.attach_lead(lead.reversed())
+ fsys = system.finalized()
+ check_fsys()
+
+ # Test the same, but with a larger scattering region.
+ system = kwant.Builder()
+ system[[chain(0), chain(1)]] = h
+ system[chain(0), chain(1)] = t
+ system.attach_lead(lead)
+ system.attach_lead(lead.reversed())
+ fsys = system.finalized()
+ check_fsys()
+
+
+# Test a more complicated graph with non-singular hopping.
+def test_graph_system(solve):
+ np.random.seed(11)
+ system = kwant.Builder()
+ lead = kwant.Builder(kwant.TranslationalSymmetry([(-1, 0)]))
+ lead.default_site_group = system.default_site_group = square
+
+ h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ h += h.conjugate().transpose()
+ h *= 0.8
+ t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+ t1 = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+ lead[0, 0] = system[[(0, 0), (1, 0)]] = h
+ lead[0, 1] = system[[(0, 1), (1, 1)]] = 4 * h
+ for builder in [system, lead]:
+ builder[(0, 0), (1, 0)] = t
+ builder[(0, 1), (1, 0)] = t1
+ builder[(0, 1), (1, 1)] = 1.1j * t1
+ system.attach_lead(lead)
+ system.attach_lead(lead.reversed())
+ fsys = system.finalized()
+
+ s, leads = solve(fsys)[: 2]
+ assert_almost_equal(np.dot(s.conjugate().transpose(), s),
+ np.identity(s.shape[0]))
+ n_modes = leads[0][2]
+ assert_equal(leads[1][2], n_modes)
+ assert_almost_equal(s[: n_modes, : n_modes], 0)
+ t_elements = np.sort(abs(np.asarray(s[n_modes :, : n_modes])),
+ axis=None)
+ t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
+ assert_almost_equal(t_elements, t_el_should_be)
+
+
+# Test a system with singular hopping.
+def test_singular_graph_system(solve):
+ np.random.seed(11)
+
+ system = kwant.Builder()
+ lead = kwant.Builder(kwant.TranslationalSymmetry([(-1, 0)]))
+ lead.default_site_group = system.default_site_group = square
+ h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ h += h.conjugate().transpose()
+ h *= 0.8
+ t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+ t1 = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
+ lead[0, 0] = system[[(0, 0), (1, 0)]] = h
+ lead[0, 1] = system[[(0, 1), (1, 1)]] = 4 * h
+ for builder in [system, lead]:
+ builder[(0, 0), (1, 0)] = t
+ builder[(0, 1), (1, 0)] = t1
+ system.attach_lead(lead)
+ system.attach_lead(lead.reversed())
+ fsys = system.finalized()
+
+ s, leads = solve(fsys)[: 2]
+ assert_almost_equal(np.dot(s.conjugate().transpose(), s),
+ np.identity(s.shape[0]))
+ n_modes = leads[0][2]
+ assert leads[1][2] == n_modes
+ assert_almost_equal(s[: n_modes, : n_modes], 0)
+ t_elements = np.sort(abs(np.asarray(s[n_modes :, : n_modes])),
+ axis=None)
+ t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
+ assert_almost_equal(t_elements, t_el_should_be)
+
+
+# This test features inside the onslice Hamiltonian a hopping matrix with more
+# zero eigenvalues than the lead hopping matrix. The previous version of the
+# sparse solver failed here.
+
+def test_tricky_singular_hopping(solve):
+ system = kwant.Builder()
+ lead = kwant.Builder(kwant.TranslationalSymmetry([(4, 0)]))
+ lead.default_site_group = system.default_site_group = square
+
+ neighbors = []
+ for i in xrange(n):
+ site = square(-1, i)
+ neighbors.append(site)
+ system[site] = 0
+ for j in xrange(4):
+ lead[j, i] = 0
+ for i in xrange(n-1):
+ system[(-1, i), (-1, i+1)] = -1
+ for j in xrange(4):
+ lead[(j, i), (j, i+1)] = -1
+ for i in xrange(n):
+ for j in xrange(4):
+ lead[(j, i), (j+1, i)] = -1
+ del lead[(1, 0), (2, 0)]
+
+ system.leads.append(kwant.builder.BuilderLead(lead, neighbors))
+ fsys = system.finalized()
+
+ s = solve(fsys, -1.3)[0]
+ assert_almost_equal(np.dot(s.conjugate().transpose(), s),
+ np.identity(s.shape[0]))
+
+
+# Test equivalence between self-energy and scattering matrix representations.
+# Also check that transmission works.
+
+def test_self_energy(solve):
+ class LeadWithOnlySelfEnergy(object):
+ def __init__(self, lead):
+ self.lead = lead
+
+ def self_energy(self, energy):
+ return self.lead.self_energy(energy)
+
+ np.random.seed(4)
+ system = kwant.Builder()
+ left_lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
+ right_lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
+ for b, site in [(system, chain(0)), (system, chain(1)),
+ (left_lead, chain(0)), (right_lead, chain(0))]:
+ h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ h += h.conjugate().transpose()
+ b[site] = h
+ for b, hopp in [(system, (chain(0), chain(1))),
+ (left_lead, (chain(0), chain(1))),
+ (right_lead, (chain(0), chain(1)))]:
+ b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ system.attach_lead(left_lead)
+ system.attach_lead(right_lead)
+ fsys = system.finalized()
+
+ t = solve(fsys, 0, [1], [0]).data
+ eig_should_be = np.linalg.eigvals(t * t.conjugate().transpose())
+ n_eig = len(eig_should_be)
+
+ fsys.leads[1] = LeadWithOnlySelfEnergy(fsys.leads[1])
+ sol = solve(fsys, 0, [1], [0])
+ ttdagnew = sol._a_ttdagger_a_inv(1, 0)
+ eig_are = np.linalg.eigvals(ttdagnew)
+ t_should_be = np.sum(eig_are)
+ assert_almost_equal(eig_are.imag, 0)
+ assert_almost_equal(np.sort(eig_are.real)[-n_eig :],
+ np.sort(eig_should_be.real))
+ assert_almost_equal(t_should_be, sol.transmission(1, 0))
+
+ fsys.leads[0] = LeadWithOnlySelfEnergy(fsys.leads[0])
+ sol = solve(fsys, 0, [1], [0])
+ ttdagnew = sol._a_ttdagger_a_inv(1, 0)
+ eig_are = np.linalg.eigvals(ttdagnew)
+ t_should_be = np.sum(eig_are)
+ assert_almost_equal(eig_are.imag, 0)
+ assert_almost_equal(np.sort(eig_are.real)[-n_eig :],
+ np.sort(eig_should_be.real))
+ assert_almost_equal(t_should_be, sol.transmission(1, 0))
+
+def test_self_energy_reflection(solve):
+ class LeadWithOnlySelfEnergy(object):
+ def __init__(self, lead):
+ self.lead = lead
+
+ def self_energy(self, energy):
+ return self.lead.self_energy(energy)
+
+ np.random.seed(4)
+ system = kwant.Builder()
+ left_lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
+ for b, site in [(system, chain(0)), (system, chain(1)),
+ (left_lead, chain(0))]:
+ h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ h += h.conjugate().transpose()
+ b[site] = h
+ for b, hopp in [(system, (chain(0), chain(1))),
+ (left_lead, (chain(0), chain(1)))]:
+ b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
+ system.attach_lead(left_lead)
+ fsys = system.finalized()
+
+ t = solve(fsys, 0, [0], [0])
+
+ fsys.leads[0] = LeadWithOnlySelfEnergy(fsys.leads[0])
+ sol = solve(fsys, 0, [0], [0])
+
+ assert_almost_equal(sol.transmission(0,0), t.transmission(0,0))
+
+def test_very_singular_leads(solve):
+ sys = kwant.Builder()
+ gr = kwant.lattice.Chain()
+ left_lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
+ right_lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
+ sys.default_site_group = gr
+ left_lead.default_site_group = right_lead.default_site_group = gr
+ sys[(0,)] = left_lead[(0,)] = right_lead[(0,)] = np.identity(2)
+ left_lead[(0,), (1,)] = np.zeros((2, 2))
+ right_lead[(0,), (1,)] = np.identity(2)
+ sys.attach_lead(left_lead)
+ sys.attach_lead(right_lead)
+ fsys = sys.finalized()
+ result = solve(fsys)
+ assert [i[2] for i in result[1]] == [0, 2]
+
+def test_ldos(ldos):
+ sys = kwant.Builder()
+ gr = kwant.lattice.Chain()
+ lead = kwant.Builder(kwant.TranslationalSymmetry((gr.vec((1,)),)))
+ sys.default_site_group = lead.default_site_group = gr
+ sys[(0,)] = sys[(1,)] = lead[(0,)] = 0
+ sys[(0,), (1,)] = lead[(0,), (1,)] = 1
+ sys.attach_lead(lead)
+ sys.attach_lead(lead.reversed())
+ fsys = sys.finalized()
+ assert_almost_equal(ldos(fsys, 0),
+ np.array([1, 1]) / (2 * np.pi))
diff --git a/kwant/solvers/tests/test_mumps.py b/kwant/solvers/tests/test_mumps.py
new file mode 100644
index 0000000..3b7712a
--- /dev/null
+++ b/kwant/solvers/tests/test_mumps.py
@@ -0,0 +1,92 @@
+from nose.plugins.skip import Skip, SkipTest
+from numpy.testing.decorators import skipif
+try:
+ from kwant.solvers.mumps import solve, ldos, options
+ import _test_sparse
+ _no_mumps = False
+except ImportError:
+ _no_mumps = True
+
+
+opt_list=[{},
+ {'nrhs' : 1},
+ {'nrhs' : 10},
+ {'nrhs' : 1, 'ordering' : 'amd'},
+ {'nrhs' : 10, 'sparse_rhs' : True},
+ {'nrhs' : 2, 'ordering' : 'amd', 'sparse_rhs' : True}]
+
+
+@skipif(_no_mumps)
+def test_output():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_output(solve)
+
+
+@skipif(_no_mumps)
+def test_one_lead():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_one_lead(solve)
+
+
+@skipif(_no_mumps)
+def test_smatrix_shape():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_smatrix_shape(solve)
+
+
+@skipif(_no_mumps)
+def test_two_equal_leads():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_two_equal_leads(solve)
+
+@skipif(_no_mumps)
+def test_graph_system():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_graph_system(solve)
+
+
+@skipif(_no_mumps)
+def test_singular_graph_system():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_singular_graph_system(solve)
+
+
+@skipif(_no_mumps)
+def test_tricky_singular_hopping():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_tricky_singular_hopping(solve)
+
+
+@skipif(_no_mumps)
+def test_self_energy():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_self_energy(solve)
+
+
+@skipif(_no_mumps)
+def test_self_energy_reflection():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_self_energy_reflection(solve)
+
+
+@skipif(_no_mumps)
+def test_very_singular_leads():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_very_singular_leads(solve)
+
+
+@skipif(_no_mumps)
+def test_ldos():
+ for opts in opt_list:
+ options(**opts)
+ _test_sparse.test_ldos(ldos)
diff --git a/kwant/solvers/tests/test_sparse.py b/kwant/solvers/tests/test_sparse.py
index 586fafa..9e8080d 100644
--- a/kwant/solvers/tests/test_sparse.py
+++ b/kwant/solvers/tests/test_sparse.py
@@ -1,298 +1,54 @@
-from __future__ import division
-import numpy as np
-from nose.tools import assert_raises
-from numpy.testing import assert_equal, assert_almost_equal
-import kwant
+from nose.plugins.skip import Skip, SkipTest
+from kwant.solvers.sparse import solve, ldos
+import kwant.solvers.sparse
+import _test_sparse
-# The solver has to provide full scattering matrix and labels with lead numbers
-# of each mode of the output.
-from kwant.solvers.sparse import solve
-
-n = 5
-chain = kwant.lattice.Chain()
-square = kwant.lattice.Square()
-
-# Test output sanity: that an error is raised if no output is requested,
-# and that solving for a subblock of a scattering matrix is the same as taking
-# a subblock of the full scattering matrix.
def test_output():
- np.random.seed(3)
- system = kwant.Builder()
- left_lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
- right_lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
- for b, site in [(system, chain(0)), (system, chain(1)),
- (left_lead, chain(0)), (right_lead, chain(0))]:
- h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
- h += h.conjugate().transpose()
- b[site] = h
- for b, hopp in [(system, (chain(0), chain(1))),
- (left_lead, (chain(0), chain(1))),
- (right_lead, (chain(0), chain(1)))]:
- b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
- system.attach_lead(left_lead)
- system.attach_lead(right_lead)
- fsys = system.finalized()
-
- result1 = solve(fsys)
- s, modes1 = result1
- assert s.shape == 2 * (sum(i[2] for i in modes1),)
- s1 = result1.submatrix(1, 0)
- s2, modes2 = solve(fsys, 0, [1], [0])
- assert s2.shape == (modes2[1][2], modes2[0][2])
- assert_almost_equal(s1, s2)
- assert_almost_equal(s.conjugate().transpose() * s, np.identity(s.shape[0]))
- assert_raises(ValueError, solve, fsys, 0, [])
- modes = solve(fsys)[1]
- h = fsys.leads[0].slice_hamiltonian()
- t = fsys.leads[0].inter_slice_hopping()
- modes1 = kwant.physics.modes(h, t)
- h = fsys.leads[1].slice_hamiltonian()
- t = fsys.leads[1].inter_slice_hopping()
- modes2 = kwant.physics.modes(h, t)
- assert_almost_equal(modes1[0], modes[0][0])
- assert_almost_equal(modes2[1], modes[1][1])
+ _test_sparse.test_output(solve)
-# Test that a system with one lead has unitary scattering matrix.
def test_one_lead():
- np.random.seed(3)
- system = kwant.Builder()
- lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
- for b, site in [(system, chain(0)), (system, chain(1)), (system, chain(2)),
- (lead, chain(0))]:
- h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
- h += h.conjugate().transpose()
- b[site] = h
- for b, hopp in [(system, (chain(0), chain(1))),
- (system, (chain(1), chain(2))),
- (lead, (chain(0), chain(1)))]:
- b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
- system.attach_lead(lead)
- fsys = system.finalized()
+ _test_sparse.test_one_lead(solve)
- s = solve(fsys)[0]
- assert_almost_equal(s.conjugate().transpose() * s, np.identity(s.shape[0]))
+def test_smatrix_shape():
+ _test_sparse.test_smatrix_shape(solve)
-# Test that a translationally invariant system with two leads has only
-# transmission and that transmission does not mix modes.
-def test_two_equal_leads():
- def check_fsys():
- s, leads = solve(fsys)[: 2]
- assert_almost_equal(s.conjugate().transpose() * s, np.identity(s.shape[0]))
- n_modes = leads[0][2]
- assert leads[1][2] == n_modes
- assert_almost_equal(s[: n_modes, : n_modes], 0)
- t_elements = np.sort(np.abs(np.asarray(s[n_modes :, : n_modes])),
- axis=None)
- t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
- assert_almost_equal(t_elements, t_el_should_be)
-
- np.random.seed(11)
- system = kwant.Builder()
- lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
- h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
- h += h.conjugate().transpose()
- h *= 0.8
- t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
- lead[chain(0)] = system[chain(0)] = h
- lead[chain(0), chain(1)] = t
- system.attach_lead(lead)
- system.attach_lead(lead.reversed())
- fsys = system.finalized()
- check_fsys()
- # Test the same, but with a larger scattering region.
- system = kwant.Builder()
- system[[chain(0), chain(1)]] = h
- system[chain(0), chain(1)] = t
- system.attach_lead(lead)
- system.attach_lead(lead.reversed())
- fsys = system.finalized()
- check_fsys()
+def test_two_equal_leads():
+ _test_sparse.test_two_equal_leads(solve)
-# Test a more complicated graph with non-singular hopping.
def test_graph_system():
- np.random.seed(11)
- system = kwant.Builder()
- lead = kwant.Builder(kwant.TranslationalSymmetry([(-1, 0)]))
- lead.default_site_group = system.default_site_group = square
-
- h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
- h += h.conjugate().transpose()
- h *= 0.8
- t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
- t1 = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
- lead[0, 0] = system[[(0, 0), (1, 0)]] = h
- lead[0, 1] = system[[(0, 1), (1, 1)]] = 4 * h
- for builder in [system, lead]:
- builder[(0, 0), (1, 0)] = t
- builder[(0, 1), (1, 0)] = t1
- builder[(0, 1), (1, 1)] = 1.1j * t1
- system.attach_lead(lead)
- system.attach_lead(lead.reversed())
- fsys = system.finalized()
+ _test_sparse.test_graph_system(solve)
- s, leads = solve(fsys)[: 2]
- assert_almost_equal(s.conjugate().transpose() * s, np.identity(s.shape[0]))
- n_modes = leads[0][2]
- assert_equal(leads[1][2], n_modes)
- assert_almost_equal(s[: n_modes, : n_modes], 0)
- t_elements = np.sort(np.abs(np.asarray(s[n_modes :, : n_modes])),
- axis=None)
- t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
- assert_almost_equal(t_elements, t_el_should_be)
-
-# Test a system with singular hopping.
def test_singular_graph_system():
- np.random.seed(11)
-
- system = kwant.Builder()
- lead = kwant.Builder(kwant.TranslationalSymmetry([(-1, 0)]))
- lead.default_site_group = system.default_site_group = square
- h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
- h += h.conjugate().transpose()
- h *= 0.8
- t = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
- t1 = 4 * np.random.rand(n, n) + 4j * np.random.rand(n, n)
- lead[0, 0] = system[[(0, 0), (1, 0)]] = h
- lead[0, 1] = system[[(0, 1), (1, 1)]] = 4 * h
- for builder in [system, lead]:
- builder[(0, 0), (1, 0)] = t
- builder[(0, 1), (1, 0)] = t1
- system.attach_lead(lead)
- system.attach_lead(lead.reversed())
- fsys = system.finalized()
+ _test_sparse.test_singular_graph_system(solve)
- s, leads = solve(fsys)[: 2]
- assert_almost_equal(s.conjugate().transpose() * s, np.identity(s.shape[0]))
- n_modes = leads[0][2]
- assert leads[1][2] == n_modes
- assert_almost_equal(s[: n_modes, : n_modes], 0)
- t_elements = np.sort(np.abs(np.asarray(s[n_modes :, : n_modes])),
- axis=None)
- t_el_should_be = n_modes * (n_modes - 1) * [0] + n_modes * [1]
- assert_almost_equal(t_elements, t_el_should_be)
-
-
-# This test features inside the onslice Hamiltonian a hopping matrix with more
-# zero eigenvalues than the lead hopping matrix. The previous version of the
-# sparse solver failed here.
def test_tricky_singular_hopping():
- system = kwant.Builder()
- lead = kwant.Builder(kwant.TranslationalSymmetry([(4, 0)]))
- lead.default_site_group = system.default_site_group = square
-
- neighbors = []
- for i in xrange(n):
- site = square(-1, i)
- neighbors.append(site)
- system[site] = 0
- for j in xrange(4):
- lead[j, i] = 0
- for i in xrange(n-1):
- system[(-1, i), (-1, i+1)] = -1
- for j in xrange(4):
- lead[(j, i), (j, i+1)] = -1
- for i in xrange(n):
- for j in xrange(4):
- lead[(j, i), (j+1, i)] = -1
- del lead[(1, 0), (2, 0)]
-
- system.leads.append(kwant.builder.BuilderLead(lead, neighbors))
- fsys = system.finalized()
-
- s = solve(fsys, -1.3)[0]
- assert_almost_equal(s.conjugate().transpose() * s, np.identity(s.shape[0]))
+ _test_sparse.test_tricky_singular_hopping(solve)
-# Test equivalence between self-energy and scattering matrix representations.
-# Also check that transmission and noise work.
-
def test_self_energy():
- class LeadWithOnlySelfEnergy(object):
- def __init__(self, lead):
- self.lead = lead
-
- def self_energy(self, energy):
- return self.lead.self_energy(energy)
+ _test_sparse.test_self_energy(solve)
- np.random.seed(4)
- system = kwant.Builder()
- left_lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
- right_lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
- for b, site in [(system, chain(0)), (system, chain(1)),
- (left_lead, chain(0)), (right_lead, chain(0))]:
- h = np.random.rand(n, n) + 1j * np.random.rand(n, n)
- h += h.conjugate().transpose()
- b[site] = h
- for b, hopp in [(system, (chain(0), chain(1))),
- (left_lead, (chain(0), chain(1))),
- (right_lead, (chain(0), chain(1)))]:
- b[hopp] = 10 * np.random.rand(n, n) + 1j * np.random.rand(n, n)
- system.attach_lead(left_lead)
- system.attach_lead(right_lead)
- fsys = system.finalized()
- t = solve(fsys, 0, [1], [0]).data
- eig_should_be = np.linalg.eigvals(t * t.conjugate().transpose())
- n_eig = len(eig_should_be)
+def test_self_energy_reflection():
+ _test_sparse.test_self_energy_reflection(solve)
- fsys.leads[1] = LeadWithOnlySelfEnergy(fsys.leads[1])
- sol = solve(fsys, 0, [1], [0])
- ttdagnew = sol._a_ttdagger_a_inv(1, 0)
- eig_are = np.linalg.eigvals(ttdagnew)
- t_should_be = np.sum(eig_are)
- noise_should_be = np.sum(eig_are * (1 - eig_are))
- assert_almost_equal(eig_are.imag, 0)
- assert_almost_equal(np.sort(eig_are.real)[-n_eig :],
- np.sort(eig_should_be.real))
- assert_almost_equal(t_should_be, sol.transmission(1, 0))
- assert_almost_equal(noise_should_be, sol.noise(1, 0))
-
- fsys.leads[0] = LeadWithOnlySelfEnergy(fsys.leads[0])
- sol = solve(fsys, 0, [1], [0])
- ttdagnew = sol._a_ttdagger_a_inv(1, 0)
- eig_are = np.linalg.eigvals(ttdagnew)
- t_should_be = np.sum(eig_are)
- noise_should_be = np.sum(eig_are * (1 - eig_are))
- assert_almost_equal(eig_are.imag, 0)
- assert_almost_equal(np.sort(eig_are.real)[-n_eig :],
- np.sort(eig_should_be.real))
- assert_almost_equal(t_should_be, sol.transmission(1, 0))
- assert_almost_equal(noise_should_be, sol.noise(1, 0))
def test_very_singular_leads():
- sys = kwant.Builder()
- gr = kwant.lattice.Chain()
- left_lead = kwant.Builder(kwant.TranslationalSymmetry([(-1,)]))
- right_lead = kwant.Builder(kwant.TranslationalSymmetry([(1,)]))
- sys.default_site_group = gr
- left_lead.default_site_group = right_lead.default_site_group = gr
- sys[(0,)] = left_lead[(0,)] = right_lead[(0,)] = np.identity(2)
- left_lead[(0,), (1,)] = np.zeros((2, 2))
- right_lead[(0,), (1,)] = np.identity(2)
- sys.attach_lead(left_lead)
- sys.attach_lead(right_lead)
- fsys = sys.finalized()
- result = solve(fsys)
- assert [i[2] for i in result[1]] == [0, 2]
+ _test_sparse.test_very_singular_leads(solve)
+
def test_umfpack_del():
- assert hasattr(kwant.solvers.sparse.umfpack.UmfpackContext, '__del__')
+ if kwant.solvers.sparse.uses_umfpack:
+ assert hasattr(kwant.solvers.sparse.umfpack.UmfpackContext,
+ '__del__')
+ else:
+ raise SkipTest
def test_ldos():
- sys = kwant.Builder()
- gr = kwant.lattice.Chain()
- lead = kwant.Builder(kwant.TranslationalSymmetry((gr.vec((1,)),)))
- sys.default_site_group = lead.default_site_group = gr
- sys[(0,)] = sys[(1,)] = lead[(0,)] = 0
- sys[(0,), (1,)] = lead[(0,), (1,)] = 1
- sys.attach_lead(lead)
- sys.attach_lead(lead.reversed())
- fsys = sys.finalized()
- assert_almost_equal(kwant.solvers.sparse.ldos(fsys, 0),
- np.array([1, 1]) / (2 * np.pi))
+ _test_sparse.test_ldos(ldos)
--
GitLab