From 9146ed7de0831699ecbc154c3a6b1eb93fe53419 Mon Sep 17 00:00:00 2001
From: Anton Akhmerov <anton.akhmerov@gmail.com>
Date: Sat, 24 Aug 2013 13:49:59 +0200
Subject: [PATCH] split solve -> smatrx, greens_function, same for BlockResult

---
 doc/source/images/ab_ring.py.diff             |   2 +-
 doc/source/images/graphene.py.diff            |   2 +-
 doc/source/images/quantum_well.py.diff        |   2 +-
 doc/source/images/spin_orbit.py.diff          |   2 +-
 doc/source/reference/kwant.rst                |   2 +-
 doc/source/reference/kwant.solvers.rst        |  14 +-
 doc/source/tutorial/ab_ring.py                |   2 +-
 doc/source/tutorial/graphene.py               |   2 +-
 doc/source/tutorial/quantum_well.py           |   2 +-
 doc/source/tutorial/quantum_wire.py           |   4 +-
 doc/source/tutorial/quantum_wire_revisited.py |   2 +-
 doc/source/tutorial/spin_orbit.py             |   2 +-
 .../tutorial/superconductor_transport.py      |   2 +-
 doc/source/tutorial/tutorial1.rst             |  11 +-
 doc/source/tutorial/tutorial2.rst             |   4 +-
 doc/source/tutorial/tutorial5.rst             |   4 +-
 doc/source/whatsnew/1.0.rst                   |  19 +-
 examples/advanced_construction.py             |   2 +-
 examples/square.py                            |   2 +-
 kwant/__init__.py                             |   2 +-
 kwant/physics/noise.py                        |   6 +-
 kwant/physics/tests/test_noise.py             |   4 +-
 kwant/solvers/common.py                       | 334 ++++++++++++------
 kwant/solvers/default.py                      |   5 +-
 kwant/solvers/mumps.py                        |   6 +-
 kwant/solvers/sparse.py                       |   5 +-
 kwant/solvers/tests/_test_sparse.py           |   8 +-
 kwant/solvers/tests/test_mumps.py             |  22 +-
 kwant/solvers/tests/test_sparse.py            |  22 +-
 kwant/tests/test_comprehensive.py             |   4 +-
 kwant/tests/test_system.py                    |   7 +-
 31 files changed, 317 insertions(+), 190 deletions(-)

diff --git a/doc/source/images/ab_ring.py.diff b/doc/source/images/ab_ring.py.diff
index d4f2c825..ce73ce4a 100644
--- a/doc/source/images/ab_ring.py.diff
+++ b/doc/source/images/ab_ring.py.diff
@@ -80,7 +80,7 @@
      # compute conductance
  
 @@ -87,18 +133,31 @@
-         smatrix = kwant.solve(sys, energy, args=[flux])
+         smatrix = kwant.smatrix(sys, energy, args=[flux])
          data.append(smatrix.transmission(1, 0))
  
 -    pyplot.figure()
diff --git a/doc/source/images/graphene.py.diff b/doc/source/images/graphene.py.diff
index 59f9d43e..741576cd 100644
--- a/doc/source/images/graphene.py.diff
+++ b/doc/source/images/graphene.py.diff
@@ -9,7 +9,7 @@
  
  # Define the graphene lattice
 @@ -100,22 +101,40 @@
-         smatrix = kwant.solve(sys, energy)
+         smatrix = kwant.smatrix(sys, energy)
          data.append(smatrix.transmission(0, 1))
  
 -    pyplot.figure()
diff --git a/doc/source/images/quantum_well.py.diff b/doc/source/images/quantum_well.py.diff
index 9262fd01..c1bcd61b 100644
--- a/doc/source/images/quantum_well.py.diff
+++ b/doc/source/images/quantum_well.py.diff
@@ -9,7 +9,7 @@
  
  def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
 @@ -52,19 +53,25 @@
-         smatrix = kwant.solve(sys, energy, args=[-welldepth])
+         smatrix = kwant.smatrix(sys, energy, args=[-welldepth])
          data.append(smatrix.transmission(1, 0))
  
 -    pyplot.figure()
diff --git a/doc/source/images/spin_orbit.py.diff b/doc/source/images/spin_orbit.py.diff
index d16e4a85..5f737a9d 100644
--- a/doc/source/images/spin_orbit.py.diff
+++ b/doc/source/images/spin_orbit.py.diff
@@ -9,7 +9,7 @@
  # define Pauli-matrices for convenience
  sigma_0 = tinyarray.array([[1, 0], [0, 1]])
 @@ -67,19 +68,24 @@
-         smatrix = kwant.solve(sys, energy)
+         smatrix = kwant.smatrix(sys, energy)
          data.append(smatrix.transmission(1, 0))
  
 -    pyplot.figure()
diff --git a/doc/source/reference/kwant.rst b/doc/source/reference/kwant.rst
index 13481d3d..9c3b6a92 100644
--- a/doc/source/reference/kwant.rst
+++ b/doc/source/reference/kwant.rst
@@ -38,4 +38,4 @@ From ``kwant.solvers.default``
 ------------------------------
 .. autosummary::
 
-   solve
+   smatrix
diff --git a/doc/source/reference/kwant.solvers.rst b/doc/source/reference/kwant.solvers.rst
index 586b27e4..94a3369b 100644
--- a/doc/source/reference/kwant.solvers.rst
+++ b/doc/source/reference/kwant.solvers.rst
@@ -26,18 +26,26 @@ reasons to use another.  The following functions are provided.
 .. autosummary::
    :toctree: generated/
 
-   solve
+   smatrix
+   greens_function
    wave_function
    ldos
 
-``solve`` returns an object of the following type:
+``smatrix`` returns an object of the following type:
 
 .. module:: kwant.solvers
 
 .. autosummary::
    :toctree: generated/
 
-   kwant.solvers.common.BlockResult
+   kwant.solvers.common.SMatrix
+
+The analog of ``smatrix``, ``greens_function`` accordingly returns:
+
+.. autosummary::
+   :toctree: generated/
+
+   kwant.solvers.common.GreensFunction
 
 Being just a thin wrapper around other solvers, the default solver selectively
 imports their functionality.  To find out the origin of any function in this
diff --git a/doc/source/tutorial/ab_ring.py b/doc/source/tutorial/ab_ring.py
index 8b302964..69aeb9cc 100644
--- a/doc/source/tutorial/ab_ring.py
+++ b/doc/source/tutorial/ab_ring.py
@@ -94,7 +94,7 @@ def plot_conductance(sys, energy, fluxes):
     normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
     data = []
     for flux in fluxes:
-        smatrix = kwant.solve(sys, energy, args=[flux])
+        smatrix = kwant.smatrix(sys, energy, args=[flux])
         data.append(smatrix.transmission(1, 0))
 
     pyplot.figure()
diff --git a/doc/source/tutorial/graphene.py b/doc/source/tutorial/graphene.py
index abdb3521..2c652d0b 100644
--- a/doc/source/tutorial/graphene.py
+++ b/doc/source/tutorial/graphene.py
@@ -113,7 +113,7 @@ def plot_conductance(sys, energies):
     # Compute transmission as a function of energy
     data = []
     for energy in energies:
-        smatrix = kwant.solve(sys, energy)
+        smatrix = kwant.smatrix(sys, energy)
         data.append(smatrix.transmission(0, 1))
 
     pyplot.figure()
diff --git a/doc/source/tutorial/quantum_well.py b/doc/source/tutorial/quantum_well.py
index fb901bb6..50001f70 100644
--- a/doc/source/tutorial/quantum_well.py
+++ b/doc/source/tutorial/quantum_well.py
@@ -54,7 +54,7 @@ def plot_conductance(sys, energy, welldepths):
     # Compute conductance
     data = []
     for welldepth in welldepths:
-        smatrix = kwant.solve(sys, energy, args=[-welldepth])
+        smatrix = kwant.smatrix(sys, energy, args=[-welldepth])
         data.append(smatrix.transmission(1, 0))
 
     pyplot.figure()
diff --git a/doc/source/tutorial/quantum_wire.py b/doc/source/tutorial/quantum_wire.py
index b1067a78..35e1ac6a 100644
--- a/doc/source/tutorial/quantum_wire.py
+++ b/doc/source/tutorial/quantum_wire.py
@@ -97,8 +97,8 @@ data = []
 for ie in xrange(100):
     energy = ie * 0.01
 
-    # compute the scattering matrix at energy energy
-    smatrix = kwant.solve(sys, energy)
+    # compute the scattering matrix at a given energy
+    smatrix = kwant.smatrix(sys, energy)
 
     # compute the transmission probability from lead 0 to
     # lead 1
diff --git a/doc/source/tutorial/quantum_wire_revisited.py b/doc/source/tutorial/quantum_wire_revisited.py
index 84ac6afc..c5242568 100644
--- a/doc/source/tutorial/quantum_wire_revisited.py
+++ b/doc/source/tutorial/quantum_wire_revisited.py
@@ -54,7 +54,7 @@ def plot_conductance(sys, energies):
     # Compute conductance
     data = []
     for energy in energies:
-        smatrix = kwant.solve(sys, energy)
+        smatrix = kwant.smatrix(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
     pyplot.figure()
diff --git a/doc/source/tutorial/spin_orbit.py b/doc/source/tutorial/spin_orbit.py
index f3b9f2ca..8f3a35a6 100644
--- a/doc/source/tutorial/spin_orbit.py
+++ b/doc/source/tutorial/spin_orbit.py
@@ -72,7 +72,7 @@ def plot_conductance(sys, energies):
     # Compute conductance
     data = []
     for energy in energies:
-        smatrix = kwant.solve(sys, energy)
+        smatrix = kwant.smatrix(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
     pyplot.figure()
diff --git a/doc/source/tutorial/superconductor_transport.py b/doc/source/tutorial/superconductor_transport.py
index c9b1ad47..3f953875 100644
--- a/doc/source/tutorial/superconductor_transport.py
+++ b/doc/source/tutorial/superconductor_transport.py
@@ -87,7 +87,7 @@ def plot_conductance(sys, energies):
     # Compute conductance
     data = []
     for energy in energies:
-        smatrix = kwant.solve(sys, energy)
+        smatrix = kwant.smatrix(sys, energy)
         # Conductance is N - R_ee + R_he
         data.append(smatrix.submatrix(0, 0).shape[0] -
                     smatrix.transmission(0, 0) +
diff --git a/doc/source/tutorial/tutorial1.rst b/doc/source/tutorial/tutorial1.rst
index 86157680..71efbed5 100644
--- a/doc/source/tutorial/tutorial1.rst
+++ b/doc/source/tutorial/tutorial1.rst
@@ -142,11 +142,12 @@ its conductance as a function of energy:
     :start-after: #HIDDEN_BEGIN_buzn
     :end-before: #HIDDEN_END_buzn
 
-We use ``kwant.solve`` which is a short name for `kwant.solvers.default.solve`
-of the default solver module `kwant.solvers.default`.  ``kwant.solve`` computes
-the scattering matrix ``smatrix`` solving a sparse linear system.  ``smatrix``
-itself allows to directly compute the total transmission probability from lead
-0 to lead 1 as ``smatrix.transmission(1, 0)``.
+We use ``kwant.smatrix`` which is a short name for
+`kwant.solvers.default.smatrix` of the default solver module
+`kwant.solvers.default`.  ``kwant.smatrix`` computes the scattering matrix
+``smatrix`` solving a sparse linear system.  ``smatrix`` itself allows to
+directly compute the total transmission probability from lead 0 to lead 1 as
+``smatrix.transmission(1, 0)``.
 
 Finally we can use `matplotlib` to make a plot of the computed data
 (although writing to file and using an external viewer such as
diff --git a/doc/source/tutorial/tutorial2.rst b/doc/source/tutorial/tutorial2.rst
index ce779b88..4750451e 100644
--- a/doc/source/tutorial/tutorial2.rst
+++ b/doc/source/tutorial/tutorial2.rst
@@ -163,7 +163,7 @@ Kwant now allows us to pass a function as a value to
 For each lattice point, the corresponding site is then passed as the
 first argument to the function `onsite`. The values of any additional
 parameters, which can be used to alter the Hamiltonian matrix elements
-at a later stage, are specified later during the call to `solve`.
+at a later stage, are specified later during the call to `smatrix`.
 Note that we had to define `onsite`, as it is
 not possible to mix values and functions as in ``sys[...] = 4 * t +
 potential``.
@@ -179,7 +179,7 @@ Finally, we compute the transmission probability:
     :start-after: #HIDDEN_BEGIN_sqvr
     :end-before: #HIDDEN_END_sqvr
 
-``kwant.solve`` allows us to specify a list, `args`, that will be passed as
+``kwant.smatrix`` allows us to specify a list, `args`, that will be passed as
 additional arguments to the functions that provide the Hamiltonian matrix
 elements.  In this example we are able to solve the system for different depths
 of the potential well by passing the potential value. We obtain the result:
diff --git a/doc/source/tutorial/tutorial5.rst b/doc/source/tutorial/tutorial5.rst
index ff078636..63e891bc 100644
--- a/doc/source/tutorial/tutorial5.rst
+++ b/doc/source/tutorial/tutorial5.rst
@@ -168,9 +168,9 @@ above the gap. At the gap edge, we observe a resonant Andreev reflection.
     - It is in fact possible to separate electron and hole degrees of
       freedom in the scattering matrix, even if one uses matrices for
       these degrees of freedom. In the solve step,
-      `~kwant.solvers.common.SparseSolver.solve` returns an array containing
+      `~kwant.solvers.common.SparseSolver.smatrix` returns an array containing
       the transverse wave functions of the lead modes (in
-      `BlockResult.lead_info <kwant.solvers.common.BlockResult.lead_info>`.
+      `SMatrix.lead_info <kwant.solvers.common.SMatrix.lead_info>`.
       By inspecting the wave functions, electron and hole wave
       functions can be distinguished (they only have entries in either
       the electron part *or* the hole part. If you encounter modes
diff --git a/doc/source/whatsnew/1.0.rst b/doc/source/whatsnew/1.0.rst
index 36a09ebb..c9a3bc69 100644
--- a/doc/source/whatsnew/1.0.rst
+++ b/doc/source/whatsnew/1.0.rst
@@ -54,6 +54,13 @@ Some renames
 * ``wave_func`` has been renamed to `~kwant.solvers.default.wave_function`,
 * ``MonatomicLattice`` has been renamed to `~kwant.lattice.Monatomic`,
 * ``PolyatomicLattice`` has been renamed to `~kwant.lattice.Polyatomic`.
+* ``solve`` was split into two functions: `~kwant.solvers.default.smatrix`, and
+  `~kwant.solvers.default.greens_function`. The former calculates the
+  scattering matrix, the latter the retarded Green's function between the sites
+  adjacent to the leads. It is temporarily not possible to mix self-energy and
+  modes leads within the same system.
+* The object that contained the results, `BlockResult` was also split into
+  `~kwant.solvers.common.SMatrix` and `~kwant.solvers.common.GreensFunction`.
 
 Band structure plots
 --------------------
@@ -98,7 +105,7 @@ regular function parameters (e.g. global variables).
 
 Now the value functions may accept arbitrary arguments after the `Site`
 arguments.  These extra arguments can be specified when
-`~kwant.solvers.default.solve` is called by setting the arguments:
+`~kwant.solvers.default.smatrix` is called by setting the arguments:
 
 args
     A tuple of values to be passed as the positional arguments to the
@@ -114,10 +121,10 @@ the following prototype::
         ...
 
 then the values of `t`, `B` and `pot` for which to solve the system can be
-passed to `solve` like this::
+passed to `~kwant.solvers.default.smatrix` like this::
 
-    kwant.solve(sys, energy,
-                args=(2., 3., 4.))
+    kwant.smatrix(sys, energy,
+                  args=(2., 3., 4.))
 
 With many parameters it can be less error-prone to collect all of them into a
 single object and pass this object as the single argument.  Such a parameter
@@ -137,7 +144,7 @@ collection could be a dictionary, or a class instance, for example::
 
     params = SimpleNamespace(t=1, mu=2)
     for params.B in B_values:
-        kwant.solve(sys, energy, args=[params])
+        kwant.smatrix(sys, energy, args=[params])
 
 Arguments can be passed in an equivalent way to
 `~kwant.solvers.default.wave_function`,
@@ -168,7 +175,7 @@ space, as well as their velocities and momenta.  All these quantities were
 previously not directly available.  The second object contains the propagating
 and evanescent modes in the compressed format expected by the sparse solver
 (previously this was the sole output of `~kwant.physics.modes`).  Accordingly,
-the `lead_info` attribute of `~kwant.solvers.default.BlockResult` contains the
+the `lead_info` attribute of `~kwant.solvers.common.SMatrix` contains the
 real space information about the modes in the leads (a list of
 `~kwant.physics.PropagatingModes` objects).
 
diff --git a/examples/advanced_construction.py b/examples/advanced_construction.py
index f47d4100..2165b404 100644
--- a/examples/advanced_construction.py
+++ b/examples/advanced_construction.py
@@ -48,7 +48,7 @@ def make_system(R):
 
 def main():
     sys = make_system(100)
-    print kwant.solve(sys, 1.1, args=[0.1]).transmission(0, 1)
+    print kwant.smatrix(sys, 1.1, args=[0.1]).transmission(0, 1)
 
 
 if __name__ == '__main__':
diff --git a/examples/square.py b/examples/square.py
index 0d31155b..c7df5d54 100644
--- a/examples/square.py
+++ b/examples/square.py
@@ -99,7 +99,7 @@ class System(kwant.system.FiniteSystem):
 def main():
     sys = System((10, 5), 1)
     energies = [0.04 * i for i in xrange(100)]
-    data = [kwant.solve(sys, energy).transmission(1, 0)
+    data = [kwant.smatrix(sys, energy).transmission(1, 0)
             for energy in energies]
 
     pyplot.plot(energies, data)
diff --git a/kwant/__init__.py b/kwant/__init__.py
index 5fdc4636..0f11f86a 100644
--- a/kwant/__init__.py
+++ b/kwant/__init__.py
@@ -19,7 +19,7 @@ from .lattice import TranslationalSymmetry
 __all__.append('TranslationalSymmetry')
 
 # Make kwant.solvers.default.solve available as kwant.solve.
-solve = solvers.default.solve
+smatrix = solvers.default.smatrix
 __all__.append('solve')
 
 # Importing plotter might not work, but this does not have to be a problem --
diff --git a/kwant/physics/noise.py b/kwant/physics/noise.py
index 92b9a9cb..847e2910 100644
--- a/kwant/physics/noise.py
+++ b/kwant/physics/noise.py
@@ -16,7 +16,7 @@ def two_terminal_shotnoise(smatrix):
 
     Parameters
     ----------
-    smatrix : `~kwant.solvers.common.BlockResult` instance
+    smatrix : `~kwant.solvers.common.SMatrix` instance
         A two terminal scattering matrix.
 
     Returns
@@ -25,6 +25,10 @@ def two_terminal_shotnoise(smatrix):
         Shot noise measured in noise quanta `2 e^3 |V| / pi hbar`.
     """
 
+    if hasattr(smatrix.lead_info[0], 'shape'):
+        raise NotImplementedError("Noise expressions in terms of Green's "
+                                  "functions are not implemented.")
+
     if len(smatrix.lead_info) != 2:
         raise ValueError("Only works for two-terminal systems!")
 
diff --git a/kwant/physics/tests/test_noise.py b/kwant/physics/tests/test_noise.py
index adc71d78..18b2c3ae 100644
--- a/kwant/physics/tests/test_noise.py
+++ b/kwant/physics/tests/test_noise.py
@@ -36,7 +36,7 @@ def test_multiterminal_input():
 
     sys = twoterminal_system()
     sys.attach_lead(sys.leads[0].builder)
-    sol = kwant.solve(sys.finalized(), out_leads=[0], in_leads=[0])
+    sol = kwant.smatrix(sys.finalized(), out_leads=[0], in_leads=[0])
     assert_raises(ValueError, two_terminal_shotnoise, sol)
 
 
@@ -45,7 +45,7 @@ def test_twoterminal():
 
     fsys = twoterminal_system().finalized()
 
-    sol = kwant.solve(fsys)
+    sol = kwant.smatrix(fsys)
     t = sol.submatrix(1, 0)
     Tn = np.linalg.eigvalsh(np.dot(t, t.conj().T))
     noise_should_be = np.sum(Tn * (1 - Tn))
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index bd839057..c8b95ba5 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -6,7 +6,7 @@
 # the AUTHORS file at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-__all__ = ['SparseSolver', 'BlockResult']
+__all__ = ['SparseSolver', 'SMatrix', 'GreensFunction']
 
 from collections import namedtuple
 import abc
@@ -90,7 +90,7 @@ class SparseSolver(object):
         """
         pass
 
-    def _make_linear_sys(self, sys, in_leads, energy=0, force_realspace=False,
+    def _make_linear_sys(self, sys, in_leads, energy=0, realspace=False,
                          check_hermiticity=True, args=()):
         """Make a sparse linear system of equations defining a scattering
         problem.
@@ -104,7 +104,7 @@ class SparseSolver(object):
             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
+        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.
@@ -126,20 +126,18 @@ class SparseSolver(object):
             total number of degrees of freedom in the scattering region.
 
         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.PropagatingModes` with a corresponding format,
-            otherwise the lead self-energy matrix.
+            Contains one entry for each lead.  If `realspace=False`, this is an
+            instance of `~kwant.physics.PropagatingModes` with a corresponding
+            format, 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
+        All the leads should implement a method `modes` if `realspace=False`
+        and a method `selfenergy`.
 
+        The system of equations that is created will be described in detail
+        elsewhere.
         """
-
         splhsmat = getattr(sp, self.lhsformat + '_matrix')
         sprhsmat = getattr(sp, self.rhsformat + '_matrix')
 
@@ -170,7 +168,7 @@ class SparseSolver(object):
         lead_info = []
         for leadnum, interface in enumerate(sys.lead_interfaces):
             lead = sys.leads[leadnum]
-            if hasattr(lead, 'modes') and not force_realspace:
+            if not realspace:
                 modes, (u, ulinv, nprop, svd_v) = lead.modes(energy, args=args)
                 lead_info.append(modes)
 
@@ -275,11 +273,10 @@ class SparseSolver(object):
 
         return LinearSys(lhs, rhs, indices, num_orb), lead_info
 
-    def solve(self, sys, energy=0, out_leads=None, in_leads=None,
-              force_realspace=False, check_hermiticity=True,
-              args=()):
+    def smatrix(self, sys, energy=0, out_leads=None, in_leads=None,
+                check_hermiticity=True, args=()):
         """
-        Compute the scattering matrix or Green's function between leads.
+        Compute the scattering matrix of a system.
 
         Parameters
         ----------
@@ -296,10 +293,6 @@ class SparseSolver(object):
             Numbers of leads in which current or wave function is injected.
             None is interpreted as all leads. Default is ``None`` and means
             "all leads".
-        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.
         args : tuple, defaults to empty
@@ -307,34 +300,19 @@ class SparseSolver(object):
 
         Returns
         -------
-        output : `~kwant.solvers.common.BlockResult`
-            See the notes below and `~kwant.solvers.common.BlockResult`
+        output : `~kwant.solvers.common.SMatrix`
+            See the notes below and `~kwant.solvers.common.SMatrix`
             documentation.
 
         Notes
         -----
         This function can be used to calculate the conductance and other
         transport properties of a system.  See the documentation for its output
-        type, `~kwant.solvers.common.BlockResult`.
+        type, `~kwant.solvers.common.SMatrix`.
 
-        It returns an object encapsulating the Green's function elements
-        between the desired leads. For leads defined as a self-energy, the
-        result is just the real-space retarded Green's function between the
-        system sites interfacing the leads in `in_leads` and those interfacing
-        the leads in `out_leads`. If, as is usually the case, the leads are
-        defined as tight-binding systems, then the Green's function from
-        incoming to outgoing modes is returned (more commonly known as the
-        scattering matrix).  If some leads are defined via a self-energy, and
-        some as tight-binding systems, the result has Green's function's
-        elements between modes and sites.  The returned object also contains
-        information about the modes or self-energies of the leads.
-
-        If `force_realspace` is set to ``True`` all leads will be treated as if
-        they would be defined in terms of self energies.  The returned Green's
-        function will be thus the retarded Green's function between sites in
-        real space, just like if the tradidional RGF algorithm would have been
-        used.  Enabling this option is more computationally expensive and can
-        be less stable.
+        The returned object contains the scattering matrix elements from the
+        `in_leads` to the `out_leads` as well as information about the lead
+        modes.
 
         Both `in_leads` and `out_leads` must be sorted and may only contain
         unique entries.
@@ -353,8 +331,7 @@ class SparseSolver(object):
         if len(in_leads) == 0 or len(out_leads) == 0:
             raise ValueError("No output is requested.")
 
-        linsys, lead_info = self._make_linear_sys(sys, in_leads, energy,
-                                                  force_realspace,
+        linsys, lead_info = self._make_linear_sys(sys, in_leads, energy, False,
                                                   check_hermiticity, args)
 
         kept_vars = np.concatenate([vars for i, vars in
@@ -365,7 +342,7 @@ class SparseSolver(object):
         len_rhs = sum(i.shape[1] for i in linsys.rhs)
         len_kv = len(kept_vars)
         if not(len_rhs and len_kv):
-            return BlockResult(np.zeros((len_kv, len_rhs)),
+            return SMatrix(np.zeros((len_kv, len_rhs)),
                                lead_info, out_leads, in_leads)
 
         # See comment about zero-shaped sparse matrices at the top of common.py.
@@ -374,8 +351,90 @@ class SparseSolver(object):
         flhs = self._factorized(linsys.lhs)
         data = self._solve_linear_sys(flhs, rhs, kept_vars)
 
-        return BlockResult(data, lead_info, out_leads, in_leads)
+        return SMatrix(data, lead_info, out_leads, in_leads)
+
+    def greens_function(self, sys, energy=0, out_leads=None, in_leads=None,
+                        check_hermiticity=True, args=()):
+        """
+        Compute the retarded Green's function of the system between its leads.
+
+        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 : sequence of integers or ``None``
+            Numbers of leads where current or wave function is extracted.  None
+            is interpreted as all leads. Default is ``None`` and means "all
+            leads".
+        in_leads : sequence of integers or ``None``
+            Numbers of leads in which current or wave function is injected.
+            None is interpreted as all leads. Default is ``None`` and means
+            "all leads".
+        check_hermiticity : ``bool``
+            Check if the Hamiltonian matrices are Hermitian.
+        args : tuple, defaults to empty
+            Positional arguments to pass to the ``hamiltonian`` method.
+
+        Returns
+        -------
+        output : `~kwant.solvers.common.GreensFunction`
+            See the notes below and `~kwant.solvers.common.GreensFunction`
+            documentation.
+
+        Notes
+        -----
+        This function can be used to calculate the conductance and other
+        transport properties of a system.  It is often slower and less stable
+        than the scattering matrix-based calculation executed by
+        `~kwant.smatrix`, and is currently provided mostly for testing
+        purposes and compatibility with RGF code.
+
+        It returns an object encapsulating the Green's function elements
+        between the system sites interfacing the leads in `in_leads` and those
+        interfacing the leads in `out_leads`.  The returned object also
+        contains a list with self-energies of the leads.
+
+        Both `in_leads` and `out_leads` must be sorted and may only contain
+        unique entries.
+        """
+
+        n = len(sys.lead_interfaces)
+        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, in_leads, energy, True,
+                                                  check_hermiticity, args)
+
+        kept_vars = np.concatenate([vars for i, vars in
+                                    enumerate(linsys.indices) if i in
+                                    out_leads])
+
+        # Do not perform factorization if no calculation is to be done.
+        len_rhs = sum(i.shape[1] for i in linsys.rhs)
+        len_kv = len(kept_vars)
+        if not(len_rhs and len_kv):
+            return GreensFunction(np.zeros((len_kv, len_rhs)),
+                               lead_info, out_leads, in_leads)
+
+        # See comment about zero-shaped sparse matrices at the top of common.py.
+        rhs = sp.bmat([[i for i in linsys.rhs if i.shape[1]]],
+                      format=self.rhsformat)
+        flhs = self._factorized(linsys.lhs)
+        data = self._solve_linear_sys(flhs, rhs, kept_vars)
 
+        return GreensFunction(data, lead_info, out_leads, in_leads)
 
     def ldos(self, fsys, energy=0, args=()):
         """
@@ -481,52 +540,16 @@ class WaveFunction(object):
 
 class BlockResult(object):
     """
-    Solution of a transport problem, subblock of retarded Green's function
-    or scattering matrix.
+    Container for a linear system solution with variable grouping.
 
-    Transport properties can be easily accessed using the
-    `~BlockResult.transmission` method (don't be fooled by the name,
-    it can also compute reflection, which is just transmission from one
-    lead back into the same lead.)
-
-    `BlockResult` however also allows for a more direct access to the result:
-    The data stored in `BlockResult` is either a real-space Green's
-    function (e.g. if ``force_realspace=True`` in
-    `~kwant.solvers.default.solve`) or a scattering matrix with respect to
-    lead modes. The details of this data can be directly accessed through
-    the instance variables `data` and `lead_info`. Subblocks of data
-    corresponding to particular leads are conveniently obtained by
-    `~BlockResult.submatrix`.
-
-    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.
-    out_leads, in_leads : list of integers
-        indices of the leads where current is extracted (out) or injected
-        (in). Only those are listed for which BlockResult contains the
-        calculated result.
+    This class is not intended to be used directly.
     """
-
-    def __init__(self, data, lead_info, out_leads, in_leads):
+    def __init__(self, data, lead_info, out_leads, in_leads, sizes):
         self.data = data
         self.lead_info = lead_info
         self.out_leads = out_leads
         self.in_leads = in_leads
-
-        sizes = []
-        for i in self.lead_info:
-            if isinstance(i, physics.PropagatingModes):
-                sizes.append(len(i.momenta) // 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)
@@ -539,16 +562,15 @@ class BlockResult(object):
         return self.out_block_coords(lead_out), self.in_block_coords(lead_in)
 
     def out_block_coords(self, lead_out):
-        """Return a slice corresponding to the rows in the block corresponding
-        to lead_out
+        """Return a slice with the rows in the block corresponding to lead_out.
         """
         lead_out = self.out_leads.index(lead_out)
         return slice(self._out_offsets[lead_out],
                      self._out_offsets[lead_out + 1])
 
     def in_block_coords(self, lead_in):
-        """Return a slice corresponding to the columns in the block
-        corresponding to lead_in
+        """
+        Return a slice with the columns in the block corresponding to lead_in.
         """
         lead_in = self.in_leads.index(lead_in)
         return slice(self._in_offsets[lead_in],
@@ -558,44 +580,124 @@ class BlockResult(object):
         """Return the matrix elements from lead_in to lead_out."""
         return self.data[self.block_coords(lead_out, lead_in)]
 
+
+class SMatrix(BlockResult):
+    """A scattering matrix.
+
+    Transport properties can be easily accessed using the
+    `~SMatrix.transmission` method (don't be fooled by the name,
+    it can also compute reflection, which is just transmission from one
+    lead back into the same lead.)
+
+    `SMatrix` however also allows for a more direct access to the result: The
+    data stored in `SMatrix` is a scattering matrix with respect to lead modes
+    and these modes themselves. The details of this data can be directly
+    accessed through the instance variables `data` and `lead_info`. Subblocks
+    of data corresponding to particular leads are conveniently obtained by
+    `~SMatrix.submatrix`.
+
+    Instance Variables
+    ------------------
+    data : NumPy array
+        a matrix containing all the requested matrix elements of the scattering
+        matrix.
+    lead_info : list of data
+        a list containing `kwant.physics.PropagatingModes` for each lead.
+    out_leads, in_leads : list of integers
+        indices of the leads where current is extracted (out) or injected
+        (in). Only those are listed for which SMatrix contains the
+        calculated result.
+    """
+
+    def __init__(self, data, lead_info, out_leads, in_leads):
+        sizes = [len(i.momenta) // 2 for i in lead_info]
+        super(SMatrix, self).__init__(data, lead_info, out_leads, in_leads,
+                                      sizes)
+
+    def transmission(self, lead_out, lead_in):
+        """Return transmission from lead_in to lead_out."""
+        return np.linalg.norm(self.submatrix(lead_out, lead_in)) ** 2
+
+    def __repr__(self):
+        return "SMatrix(data=%r, lead_info=%r, " \
+            "out_leads=%r, in_leads=%r)" % (self.data, self.lead_info,
+                                            self.out_leads, self.in_leads)
+
+
+class GreensFunction(BlockResult):
+    """
+    Retarded Green's function.
+
+    Transport properties can be easily accessed using the
+    `~GreensFunction.transmission` method (don't be fooled by the name, it can
+    also compute reflection, which is just transmission from one lead back into
+    the same lead).
+
+    `GreensFunction` however also allows for a more direct access to the
+    result: The data stored in `GreensFunction` is the real-space Green's
+    function. The details of this data can be directly accessed through the
+    instance variables `data` and `lead_info`. Subblocks of data corresponding
+    to particular leads are conveniently obtained by
+    `~GreensFunction.submatrix`.
+
+    Instance Variables
+    ------------------
+    data : NumPy array
+        a matrix containing all the requested matrix elements of Green's
+        function.
+    lead_info : list of matrices
+        a list with self-energies of each lead.
+    out_leads, in_leads : list of integers
+        indices of the leads where current is extracted (out) or injected
+        (in). Only those are listed for which SMatrix contains the
+        calculated result.
+    """
+
+    def __init__(self, data, lead_info, out_leads, in_leads):
+        sizes = [i.shape[0] for i in lead_info]
+        super(GreensFunction, self).__init__(data, lead_info, out_leads,
+                                             in_leads, sizes)
+
     def _a_ttdagger_a_inv(self, lead_out, lead_in):
+        """Return t * t^dagger in a certain basis."""
         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, physics.PropagatingModes):
-                # Lead is a "self energy lead": multiply gf2 with a gamma
-                # matrix.
-                factors.append(1j * (possible_se - possible_se.conj().T))
+            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], physics.PropagatingModes) and \
-           isinstance(self.lead_info[lead_in], physics.PropagatingModes):
-            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
+        gf = self.submatrix(lead_out, lead_in)
+        factors = []
+        for lead, gf2 in ((lead_out, gf), (lead_in, gf.conj().T)):
+            self_en = self.lead_info[lead]
+            factors.append(1j * (self_en - self_en.conj().T))
+            factors.append(gf2)
+        attdagainv = reduce(np.dot, factors)
+
+        result = np.trace(attdagainv).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
 
     def __repr__(self):
-        return "BlockResult(data=%r, lead_info=%r, " \
+        return "GreensFunction(data=%r, lead_info=%r, " \
             "out_leads=%r, in_leads=%r)" % (self.data, self.lead_info,
                                             self.out_leads, self.in_leads)
diff --git a/kwant/solvers/default.py b/kwant/solvers/default.py
index e92fdbca..56260d72 100644
--- a/kwant/solvers/default.py
+++ b/kwant/solvers/default.py
@@ -6,7 +6,7 @@
 # the AUTHORS file at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-__all__ = ['solve', 'ldos', 'wave_function']
+__all__ = ['smatrx', 'ldos', 'wave_function', 'greens_function']
 
 # MUMPS usually works best.  Use SciPy as fallback.
 try:
@@ -16,6 +16,7 @@ except ImportError:
 
 hidden_instance = smodule.Solver()
 
-solve = hidden_instance.solve
+smatrix = hidden_instance.smatrix
 ldos = hidden_instance.ldos
 wave_function = hidden_instance.wave_function
+greens_function = hidden_instance.greens_function
diff --git a/kwant/solvers/mumps.py b/kwant/solvers/mumps.py
index d99e615d..ab491e4e 100644
--- a/kwant/solvers/mumps.py
+++ b/kwant/solvers/mumps.py
@@ -6,7 +6,8 @@
 # the AUTHORS file at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-__all__ = ['solve', 'ldos', 'wave_function', 'options', 'Solver']
+__all__ = ['smatrix', 'ldos', 'wave_function', 'greens_function', 'options',
+           'Solver']
 
 import numpy as np
 from . import common
@@ -122,7 +123,8 @@ class Solver(common.SparseSolver):
 
 default_solver = Solver()
 
-solve = default_solver.solve
+smatrix = default_solver.smatrix
+greens_function = default_solver.greens_function
 ldos = default_solver.ldos
 wave_function = default_solver.wave_function
 options = default_solver.options
diff --git a/kwant/solvers/sparse.py b/kwant/solvers/sparse.py
index 3c7f157d..b2712f38 100644
--- a/kwant/solvers/sparse.py
+++ b/kwant/solvers/sparse.py
@@ -6,7 +6,7 @@
 # the AUTHORS file at the top-level directory of this distribution and at
 # http://kwant-project.org/authors.
 
-__all__ = ['solve', 'ldos', 'wave_function', 'Solver']
+__all__ = ['smatrix', 'greens_function', 'ldos', 'wave_function', 'Solver']
 
 import warnings
 import numpy as np
@@ -115,6 +115,7 @@ class Solver(common.SparseSolver):
 
 default_solver = Solver()
 
-solve = default_solver.solve
+smatrix = default_solver.smatrix
+greens_function = default_solver.greens_function
 ldos = default_solver.ldos
 wave_function = default_solver.wave_function
diff --git a/kwant/solvers/tests/_test_sparse.py b/kwant/solvers/tests/_test_sparse.py
index 9f8b2d83..1a3e49af 100644
--- a/kwant/solvers/tests/_test_sparse.py
+++ b/kwant/solvers/tests/_test_sparse.py
@@ -281,7 +281,7 @@ def test_tricky_singular_hopping(solve):
 
 # Test equivalence between self-energy and scattering matrix representations.
 # Also check that transmission works.
-def test_selfenergy(solve):
+def test_selfenergy(greens_function, smatrix):
     np.random.seed(4)
     system = kwant.Builder()
     left_lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
@@ -299,12 +299,12 @@ def test_selfenergy(solve):
     system.attach_lead(right_lead)
     fsys = system.finalized()
 
-    t = solve(fsys, 0, [1], [0]).data
+    t = smatrix(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])
+    sol = greens_function(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)
@@ -314,7 +314,7 @@ def test_selfenergy(solve):
     assert_almost_equal(t_should_be, sol.transmission(1, 0))
 
     fsys.leads[0] = LeadWithOnlySelfEnergy(fsys.leads[0])
-    sol = solve(fsys, 0, [1], [0])
+    sol = greens_function(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)
diff --git a/kwant/solvers/tests/test_mumps.py b/kwant/solvers/tests/test_mumps.py
index 3c78eb7c..9802e1fe 100644
--- a/kwant/solvers/tests/test_mumps.py
+++ b/kwant/solvers/tests/test_mumps.py
@@ -10,7 +10,7 @@
 from numpy.testing.decorators import skipif
 try:
     from kwant.solvers.mumps import \
-        solve, ldos, wave_function, options, reset_options
+        smatrix, greens_function, ldos, wave_function, options, reset_options
     import _test_sparse
     _no_mumps = False
 except ImportError:
@@ -30,7 +30,7 @@ def test_output():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_output(solve)
+        _test_sparse.test_output(smatrix)
 
 
 @skipif(_no_mumps)
@@ -38,7 +38,7 @@ def test_one_lead():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_one_lead(solve)
+        _test_sparse.test_one_lead(smatrix)
 
 
 @skipif(_no_mumps)
@@ -46,7 +46,7 @@ def test_smatrix_shape():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_smatrix_shape(solve)
+        _test_sparse.test_smatrix_shape(smatrix)
 
 
 @skipif(_no_mumps)
@@ -54,14 +54,14 @@ def test_two_equal_leads():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_two_equal_leads(solve)
+        _test_sparse.test_two_equal_leads(smatrix)
 
 @skipif(_no_mumps)
 def test_graph_system():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_graph_system(solve)
+        _test_sparse.test_graph_system(smatrix)
 
 
 @skipif(_no_mumps)
@@ -69,7 +69,7 @@ def test_singular_graph_system():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_singular_graph_system(solve)
+        _test_sparse.test_singular_graph_system(smatrix)
 
 
 @skipif(_no_mumps)
@@ -77,7 +77,7 @@ def test_tricky_singular_hopping():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_tricky_singular_hopping(solve)
+        _test_sparse.test_tricky_singular_hopping(smatrix)
 
 
 @skipif(_no_mumps)
@@ -85,7 +85,7 @@ def test_selfenergy():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_selfenergy(solve)
+        _test_sparse.test_selfenergy(greens_function, smatrix)
 
 
 @skipif(_no_mumps)
@@ -93,7 +93,7 @@ def test_selfenergy_reflection():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_selfenergy_reflection(solve)
+        _test_sparse.test_selfenergy_reflection(greens_function)
 
 
 @skipif(_no_mumps)
@@ -101,7 +101,7 @@ def test_very_singular_leads():
     for opts in opt_list:
         reset_options()
         options(**opts)
-        _test_sparse.test_very_singular_leads(solve)
+        _test_sparse.test_very_singular_leads(smatrix)
 
 
 @skipif(_no_mumps)
diff --git a/kwant/solvers/tests/test_sparse.py b/kwant/solvers/tests/test_sparse.py
index 8b0a2375..399a1a9c 100644
--- a/kwant/solvers/tests/test_sparse.py
+++ b/kwant/solvers/tests/test_sparse.py
@@ -7,48 +7,48 @@
 # http://kwant-project.org/authors.
 
 from nose.plugins.skip import Skip, SkipTest
-from  kwant.solvers.sparse import solve, ldos, wave_function
+from  kwant.solvers.sparse import smatrix, greens_function, ldos, wave_function
 import kwant.solvers.sparse
 import _test_sparse
 
 def test_output():
-    _test_sparse.test_output(solve)
+    _test_sparse.test_output(smatrix)
 
 
 def test_one_lead():
-    _test_sparse.test_one_lead(solve)
+    _test_sparse.test_one_lead(smatrix)
 
 
 def test_smatrix_shape():
-    _test_sparse.test_smatrix_shape(solve)
+    _test_sparse.test_smatrix_shape(smatrix)
 
 
 def test_two_equal_leads():
-    _test_sparse.test_two_equal_leads(solve)
+    _test_sparse.test_two_equal_leads(smatrix)
 
 
 def test_graph_system():
-    _test_sparse.test_graph_system(solve)
+    _test_sparse.test_graph_system(smatrix)
 
 
 def test_singular_graph_system():
-    _test_sparse.test_singular_graph_system(solve)
+    _test_sparse.test_singular_graph_system(smatrix)
 
 
 def test_tricky_singular_hopping():
-    _test_sparse.test_tricky_singular_hopping(solve)
+    _test_sparse.test_tricky_singular_hopping(smatrix)
 
 
 def test_selfenergy():
-    _test_sparse.test_selfenergy(solve)
+    _test_sparse.test_selfenergy(greens_function, smatrix)
 
 
 def test_selfenergy_reflection():
-    _test_sparse.test_selfenergy_reflection(solve)
+    _test_sparse.test_selfenergy_reflection(greens_function)
 
 
 def test_very_singular_leads():
-    _test_sparse.test_very_singular_leads(solve)
+    _test_sparse.test_very_singular_leads(smatrix)
 
 
 def test_umfpack_del():
diff --git a/kwant/tests/test_comprehensive.py b/kwant/tests/test_comprehensive.py
index 0bb59a05..4b0086f7 100644
--- a/kwant/tests/test_comprehensive.py
+++ b/kwant/tests/test_comprehensive.py
@@ -41,7 +41,7 @@ def test_qhe(W=16, L=8):
     # reciprocal_phis = numpy.linspace(0.1, 7, 200)
     # conductances = []
     # for phi in 1 / reciprocal_phis:
-    #     smatrix = kwant.solve(sys, 1.0, args=[phi, ""])
+    #     smatrix = kwant.smatrix(sys, 1.0, args=[phi, ""])
     #     conductances.append(smatrix.transmission(1, 0))
     # pyplot.plot(reciprocal_phis, conductances)
     # pyplot.show()
@@ -51,5 +51,5 @@ def test_qhe(W=16, L=8):
                                        ((5.2, 5.5), 3, 1e-1)]:
         for r_phi in r_phis:
             args = (1.0 / r_phi, "")
-            T_actual = kwant.solve(sys, 1.0, args=args).transmission(1, 0)
+            T_actual = kwant.smatrix(sys, 1.0, args=args).transmission(1, 0)
             assert abs(T_nominal - T_actual) < max_err
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index 898cdf69..8359b9ab 100644
--- a/kwant/tests/test_system.py
+++ b/kwant/tests/test_system.py
@@ -63,11 +63,12 @@ def test_hamiltonian_submatrix():
     lead[gr(0), gr(1)] = np.random.randn(2, 2)
     sys.attach_lead(lead)
     sys2 = sys.finalized()
-    smatrix = kwant.solve(sys2, .1).data
+    smatrix = kwant.smatrix(sys2, .1).data
     sys3 = sys2.precalculate(.1, calculate_selfenergy=False)
-    smatrix2 = kwant.solve(sys3, .1).data
+    smatrix2 = kwant.smatrix(sys3, .1).data
     np.testing.assert_almost_equal(smatrix, smatrix2)
-    assert_raises(ValueError, kwant.solve, sys3, 0.2, None, None, True)
+    assert_raises(ValueError, kwant.solvers.default.greens_function, sys3, 0.2,
+                  None, None)
 
     # Test for shape errors.
     sys[gr(0), gr(2)] = np.array([[1, 2]])
-- 
GitLab