allow to choose modes algorithm, add a systematic (failing) test for modes

kwant/physics/leads.py

@@ -21,7 +21,7 @@ __all__ = ['selfenergy', 'modes', 'ModesTuple', 'selfenergy_from_modes']
 
 Linsys = namedtuple('Linsys', ['eigenproblem', 'v', 'extract', 'project'])
 
-def setup_linsys(h_onslice, h_hop, tol=1e6):
+def setup_linsys(h_onslice, h_hop, tol=1e6, algorithm=None):
     """
     Make an eigenvalue problem for eigenvectors of translation operator.
 
@@ -34,6 +34,18 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
     tol : float
         Numbers are considered zero when they are smaller than `tol` times
         the machine precision.
+    algorithm : tuple or 3 booleans or None
+        Which steps of the eigenvalue prolem stabilization to perform.
+        The first value selects whether to work in the basis of the
+        hopping svd, or lattice basis. If the real space basis is chosen, the
+        following two options do not apply.
+        The second value selects whether to add an anti-Hermitian term to the
+        slice Hamiltonian before inverting. Finally the third value selects
+        whether to reduce a generalized eigenvalue problem to the regular one.
+        The default value, `None`, results in kwant selecting the algorithm
+        that is the most efficient without sacrificing stability. Manual
+        selection may result in either slower performance, or large numerical
+        errors, and is mostly required for testing purposes.
 
     Returns
     -------
@@ -69,7 +81,8 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
     assert m == vh.shape[1], "Corrupt output of svd."
     n_nonsing = np.sum(s > eps * s[0])
 
-    if n_nonsing == n:
+    if (n_nonsing == n and algorithm is None) or (algorithm is not None and
+                                                  algorithm[0]):
         # The hopping matrix is well-conditioned and can be safely inverted.
         # Hence the regular transfer matrix may be used.
         hop_inv = la.inv(h_hop)
@@ -94,6 +107,11 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
         matrices = (A, None)
         v_out = None
     else:
+        if algorithm is not None:
+            need_to_stabilize, divide = algorithm[1:]
+        else:
+            need_to_stabilize = None
+
         # The hopping matrix has eigenvalues close to 0 - those
         # need to be eliminated.
 
@@ -115,21 +133,23 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
         # always if the original Hamiltonian is complex, and check for
         # invertibility first if it is real
 
-        need_to_stabilize = True
-
-        if issubclass(np.common_type(h_onslice, h_hop), np.floating):
+        h = h_onslice
+        sol = kla.lu_factor(h)
+        if issubclass(np.common_type(h_onslice, h_hop), np.floating) \
+           and need_to_stabilize is None:
             # Check if stabilization is needed.
-            h = h_onslice
-            sol = kla.lu_factor(h)
             rcond = kla.rcond_from_lu(sol, npl.norm(h, 1))
 
             if rcond > eps:
                 need_to_stabilize = False
 
+        if need_to_stabilize is None:
+            need_to_stabilize = True
+
         if need_to_stabilize:
             # Matrices are complex or need self-energy-like term to be
             # stabilized.
-            temp = dot(u, s.reshape(-1, 1) * u.T.conj()) + dot(v, v.T.conj())
+            temp = dot(u * s, u.T.conj()) + dot(v, v.T.conj())
             h = h_onslice + 1j * temp
 
             sol = kla.lu_factor(h)
@@ -163,30 +183,34 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
         A = np.zeros((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
         B = np.zeros((2 * n_nonsing, 2 * n_nonsing), np.common_type(h, h_hop))
 
-        A[:n_nonsing, :n_nonsing] = -np.identity(n_nonsing)
+        begin, end = slice(n_nonsing), slice(n_nonsing, None)
+
+        A[begin, begin] = -np.identity(n_nonsing)
 
-        B[n_nonsing:, n_nonsing:] = np.identity(n_nonsing)
+        B[end, end] = np.identity(n_nonsing)
         u_s = u * s
 
         temp = kla.lu_solve(sol, v)
         temp2 = dot(u_s.T.conj(), temp)
         if need_to_stabilize:
-            A[n_nonsing:, :n_nonsing] = 1j * temp2
-        A[n_nonsing:, n_nonsing:] = -temp2
+            A[end, begin] = 1j * temp2
+        A[end, end] = -temp2
         temp2 = dot(v.T.conj(), temp)
         if need_to_stabilize:
-            A[:n_nonsing, :n_nonsing] += 1j *temp2
-        A[:n_nonsing, n_nonsing:] = -temp2
+            A[begin, begin] += 1j *temp2
+        A[begin, end] = -temp2
 
         temp = kla.lu_solve(sol, u_s)
         temp2 = dot(u_s.T.conj(), temp)
-        B[n_nonsing:, :n_nonsing] = temp2
+        B[end, begin] = temp2
         if need_to_stabilize:
-            B[n_nonsing:, n_nonsing:] -= 1j * temp2
+            B[end, end] -= 1j * temp2
         temp2 = dot(v.T.conj(), temp)
-        B[:n_nonsing, :n_nonsing] = temp2
+        B[begin, begin] = temp2
         if need_to_stabilize:
-            B[:n_nonsing, n_nonsing:] = -1j * temp2
+            B[begin, end] = -1j * temp2
+
+        v_out = v[:m]
 
         # Solving a generalized eigenproblem is about twice as expensive
         # as solving a regular eigenvalue problem.
@@ -197,16 +221,16 @@ def setup_linsys(h_onslice, h_hop, tol=1e6):
         # the matrix B can be safely inverted.
 
         lu_b = kla.lu_factor(B)
-        rcond = kla.rcond_from_lu(lu_b, npl.norm(B, 1))
+        if algorithm is None:
+            rcond = kla.rcond_from_lu(lu_b, npl.norm(B, 1))
+            # A more stringent condition is used here since errors can
+            # accumulate from here to the eigenvalue calculation later.
+            divide = rcond > eps * tol
 
-        # A more stringent condition is used here since errors can accumulate
-        # from here to the eigenvalue calculation later.
-        v_out = v[:m]
-        if rcond > eps * tol:
+        if divide:
             matrices = (kla.lu_solve(lu_b, A), None)
         else:
             matrices = (A, B)
-
     return Linsys(matrices, v_out, extract_wf, project_wf)
 
 
@@ -398,7 +422,7 @@ d.update(__doc__=modes_docstring)
 ModesTuple = type('ModesTuple', _Modes.__bases__, d)
 del _Modes, d, modes_docstring
 
-def modes(h_onslice, h_hop, tol=1e6):
+def modes(h_onslice, h_hop, tol=1e6, algorithm=None):
     """
     Compute the eigendecomposition of a translation operator of a lead.
 
@@ -412,6 +436,18 @@ def modes(h_onslice, h_hop, tol=1e6):
     tol : float
         Numbers and differences are considered zero when they are smaller
         than `tol` times the machine precision.
+    algorithm : tuple or 3 booleans or None
+        Which steps of the eigenvalue prolem stabilization to perform. The
+        default value, `None`, results in kwant selecting the algorithm that is
+        the most efficient without sacrificing stability. Manual selection may
+        result in either slower performance, or large numerical errors, and is
+        mostly required for testing purposes.  The first value selects whether
+        to work in the basis of the hopping svd, or lattice basis. If the real
+        space basis is chosen, the following two options do not apply.  The
+        second value selects whether to add an anti-Hermitian term to the slice
+        Hamiltonian before inverting. Finally the third value selects whether
+        to reduce a generalized eigenvalue problem to the regular one.
+
 
     Returns
     -------
@@ -447,6 +483,7 @@ def modes(h_onslice, h_hop, tol=1e6):
 
     This function uses the most stable and efficient algorithm for calculating
     the mode decomposition. Its details will be published elsewhere.
+
     """
     m = h_hop.shape[1]
 
@@ -462,7 +499,8 @@ def modes(h_onslice, h_hop, tol=1e6):
         return ModesTuple(np.empty((0, 0)), np.empty((0, 0)), 0, v)
 
     # Defer most of the calculation to helper routines.
-    matrices, v, extract, project = setup_linsys(h_onslice, h_hop, tol)
+    matrices, v, extract, project = setup_linsys(h_onslice, h_hop, tol,
+                                                 algorithm)
     ev, evanselect, propselect, vec_gen, ord_schur =\
          unified_eigenproblem(*(matrices + (tol,)))
 

kwant/physics/tests/test_leads.py

@@ -8,6 +8,7 @@
 
 from __future__ import division
 import numpy as np
+from itertools import product
 from numpy.testing import assert_almost_equal
 from kwant.physics import leads
 import kwant
@@ -54,6 +55,7 @@ def test_regular_fully_degenerate():
 
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
+
 def test_regular_degenerate_with_crossing():
     """This is a testcase with invertible hopping matrices,
     and degenerate k-values with a crossing such that one
@@ -84,6 +86,7 @@ def test_regular_degenerate_with_crossing():
 
     assert_almost_equal(g, modes_se(h_onslice, hop))
 
+
 def test_singular():
     """This testcase features a rectangular (and hence singular)
      hopping matrix without degeneracies.
@@ -107,9 +110,9 @@ def test_singular():
     h_onslice[w:, w:] = h_onslice_s
     g = leads.square_selfenergy(w, t, e)
 
-    print np.round(g, 5) / np.round(modes_se(h_onslice, h_hop), 5)
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
+
 def test_singular_but_square():
     """This testcase features a singular, square hopping matrices
     without degeneracies.
@@ -136,6 +139,7 @@ def test_singular_but_square():
     g[:w, :w] = leads.square_selfenergy(w, t, e)
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
+
 def test_singular_fully_degenerate():
     """This testcase features a rectangular (and hence singular)
      hopping matrix with complete degeneracy.
@@ -169,6 +173,7 @@ def test_singular_fully_degenerate():
 
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
+
 def test_singular_degenerate_with_crossing():
     """This testcase features a rectangular (and hence singular)
      hopping matrix with degeneracy k-values including a crossing
@@ -203,6 +208,7 @@ def test_singular_degenerate_with_crossing():
 
     assert_almost_equal(g, modes_se(h_onslice, h_hop))
 
+
 def test_singular_h_and_t():
     h = 0.1 * np.identity(6)
     t = np.eye(6, 6, 4)
@@ -211,6 +217,7 @@ def test_singular_h_and_t():
     sigma_should_be[4, 4] = sigma_should_be[5, 5] = -10
     assert_almost_equal(sigma, sigma_should_be)
 
+
 def test_modes():
     h, t = .3, .7
     vecs, vecslinv, nrprop, svd = leads.modes(np.array([[h]]), np.array([[t]]))
@@ -219,7 +226,8 @@ def test_modes():
     np.testing.assert_almost_equal((vecs[0] *  vecslinv[0].conj()).imag,
                                    [0.5, -0.5])
 
-def test_modes_bearder_ribbon():
+
+def test_modes_bearded_ribbon():
     # Check if bearded graphene ribbons work.
     lat = kwant.lattice.honeycomb()
     sys = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
@@ -229,3 +237,19 @@ def test_modes_bearder_ribbon():
     sys = sys.finalized()
     h, t = sys.slice_hamiltonian(), sys.inter_slice_hopping()
     assert leads.modes(h, t).nmodes == 8  # Calculated by plotting dispersion.
+
+
+def test_algorithm_equivalence():
+    np.random.seed(400)
+    n = 5
+    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)

kwant/solvers/tests/_test_sparse.py

@@ -55,7 +55,6 @@ def test_output(solve):
     s2, modes2 = result2.data, result2.lead_info
     assert s2.shape == (modes2[1][2], modes2[0][2])
     assert_almost_equal(s1, s2)
-    print np.dot(s.T.conj(), s)
     assert_almost_equal(np.dot(s.T.conj(), s),
                         np.identity(s.shape[0]))
     assert_raises(ValueError, solve, fsys, 0, [])