make modes pass the test

kwant/physics/leads.py

@@ -119,6 +119,7 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
         # u, v are matrices with shape n x n_nonsing.
         u = u[:, :n_nonsing]
         s = s[:n_nonsing]
+        u_s = u * s
         # pad v with zeros if necessary
         v = np.zeros((n, n_nonsing), dtype=vh.dtype)
         v[:vh.shape[1]] = vh[:n_nonsing].T.conj()
@@ -149,7 +150,7 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
         if need_to_stabilize:
             # Matrices are complex or need self-energy-like term to be
             # stabilized.
-            temp = dot(u * s, u.T.conj()) + dot(v, v.T.conj())
+            temp = dot(u_s, u_s.T.conj()) + dot(v, v.T.conj())
             h = h_onslice + 1j * temp
 
             sol = kla.lu_factor(h)
@@ -165,7 +166,7 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
         # the projected one (v^dagger psi lambda^-1, s u^dagger psi).
 
         def extract_wf(psi, lmbdainv):
-            wf = - dot(u * s, psi[: n_nonsing] * lmbdainv) - \
+            wf = - dot(u_s, psi[: n_nonsing] * lmbdainv) - \
                  dot(v, psi[n_nonsing:])
             if need_to_stabilize:
                 wf += 1j * (dot(v, psi[: n_nonsing]) +
@@ -185,30 +186,27 @@ def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
 
         begin, end = slice(n_nonsing), slice(n_nonsing, None)
 
-        A[begin, begin] = -np.identity(n_nonsing)
-
-        B[end, end] = np.identity(n_nonsing)
-        u_s = u * s
-
+        A[end, begin] = np.identity(n_nonsing)
         temp = kla.lu_solve(sol, v)
         temp2 = dot(u_s.T.conj(), temp)
         if need_to_stabilize:
-            A[end, begin] = 1j * temp2
-        A[end, end] = -temp2
+            A[begin, begin] = -1j * temp2
+        A[begin, end] = temp2
         temp2 = dot(v.T.conj(), temp)
         if need_to_stabilize:
-            A[begin, begin] += 1j *temp2
-        A[begin, end] = -temp2
+            A[end, begin] -= 1j *temp2
+        A[end, end] = temp2
 
+        B[begin, end] = -np.identity(n_nonsing)
         temp = kla.lu_solve(sol, u_s)
         temp2 = dot(u_s.T.conj(), temp)
-        B[end, begin] = temp2
+        B[begin, begin] = -temp2
         if need_to_stabilize:
-            B[end, end] -= 1j * temp2
+            B[begin, end] += 1j * temp2
         temp2 = dot(v.T.conj(), temp)
-        B[begin, begin] = temp2
+        B[end, begin] = -temp2
         if need_to_stabilize:
-            B[begin, end] = -1j * temp2
+            B[end, end] = 1j * temp2
 
         v_out = v[:m]
 

kwant/physics/tests/test_leads.py

@@ -8,7 +8,7 @@
 
 from __future__ import division
 import numpy as np
-from itertools import product
+from itertools import product, izip
 from numpy.testing import assert_almost_equal
 from kwant.physics import leads
 import kwant
@@ -241,15 +241,38 @@ def test_modes_bearded_ribbon():
 
 def test_algorithm_equivalence():
     np.random.seed(400)
-    n = 5
+    n = 12
     h = np.random.randn(n, n) + 1j * np.random.randn(n, n)
     h += h.T.conj()
     t = np.random.randn(n, n) + 1j * np.random.randn(n, n)
-    results = [leads.modes(h, t, algorithm=algo)
-               for algo in product(*(3 * [(True, False)]))]
-    for i in results:
-        assert np.allclose(results[0].vecs, i.vecs)
-        vecslmbdainv = i.vecslmbdainv
-        if i.svd is not None:
-            vecslmbdainv = np.dot(i.svd, vecslmbdainv)
-        assert np.allclose(vecslmbdainv, results[0].vecslmbdainv)
+    u, s, vh = np.linalg.svd(t)
+    prop_vecs = []
+    evan_vecs = []
+    algos = [(True,)] + list(product(*([(False,)] + 2 * [(True, False)])))
+    for algo in algos:
+        result = leads.modes(h, t, algorithm=algo)
+
+        vecs, vecslmbdainv = result.vecs, result.vecslmbdainv
+
+        # Bring the calculated vectors to real space
+        if not algo[0]:
+            vecslmbdainv = np.dot(vh.T.conj(), vecslmbdainv)
+            vecs = np.dot(vh.T.conj(), vecs)
+            np.testing.assert_almost_equal(result.svd, vh.T.conj())
+        full_vecs = np.r_[vecslmbdainv, vecs]
+
+        prop_vecs.append(full_vecs[:, : 2 * result.nmodes])
+        evan_vecs.append(full_vecs[:, 2 * result.nmodes :])
+
+    msg = 'Algorithm {0} failed.'
+    for vecs, algo in izip(prop_vecs, algos):
+        # Propagating modes should have identical ordering, and only vary
+        # By a phase
+        np.testing.assert_allclose(np.abs(np.sum(vecs/prop_vecs[0],
+                                                 axis=0)), vecs.shape[0],
+                                   err_msg=msg.format(algo))
+
+    for vecs, algo in izip(evan_vecs, algos):
+        # Evanescent modes must span the same linear space.
+        assert np.linalg.matrix_rank(np.c_[vecs, evan_vecs[0]], tol=1e-12) == \
+               vecs.shape[1], msg.format(algo)