From 4dd8bb7d1fb60b935245581b30924051d38ced5e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?T=C3=B3mas=20Rosdahl?= <torosdahl@gmail.com>
Date: Tue, 19 Feb 2019 17:12:27 +0100
Subject: [PATCH] add sparse option and a heuristic for using it in symmetry
 analysis

---
 kwant/qsymm.py            | 16 ++++++++++++++--
 kwant/tests/test_qsymm.py | 19 ++++++++++++++++---
 2 files changed, 30 insertions(+), 5 deletions(-)

diff --git a/kwant/qsymm.py b/kwant/qsymm.py
index 5ae0f449..93f4bdc3 100644
--- a/kwant/qsymm.py
+++ b/kwant/qsymm.py
@@ -409,7 +409,8 @@ def _get_builder_symmetries(builder):
 
 
 def find_builder_symmetries(builder, momenta=None, params=None,
-                            spatial_symmetries=True, prettify=True):
+                            spatial_symmetries=True, prettify=True,
+                            sparse=None):
     """Finds the symmetries of a Kwant system using qsymm.
 
     Parameters
@@ -430,6 +431,12 @@ def find_builder_symmetries(builder, momenta=None, params=None,
         Whether to carry out sparsification of the continuous symmetry
         generators, in general an arbitrary linear combination of the
         symmetry generators is returned.
+    sparse : bool, or None (default None)
+        Whether to use sparse linear algebra in the calculation.
+        Can give large performance gain in large systems.
+        If None, uses sparse or dense computation depending on
+        the size of the Hamiltonian.
+
 
     Returns
     -------
@@ -443,6 +450,10 @@ def find_builder_symmetries(builder, momenta=None, params=None,
     ham = builder_to_model(builder, momenta=momenta,
                            real_space=False, params=params)
 
+    # Use dense or sparse computation depending on Hamiltonian size
+    if sparse is None:
+        sparse = list(ham.values())[0].shape[0] > 20
+
     dim = len(np.array(builder.symmetry.periods))
 
     if spatial_symmetries:
@@ -455,5 +466,6 @@ def find_builder_symmetries(builder, momenta=None, params=None,
             qsymm.PointGroupElement(np.eye(dim), True, True, None),   # P
             qsymm.PointGroupElement(np.eye(dim), False, True, None)]  # C
     sg, cg = qsymm.symmetries(ham, candidates, prettify=prettify,
-                              continuous_rotations=False)
+                              continuous_rotations=False,
+                              sparse_linalg=sparse)
     return list(sg) + list(cg)
diff --git a/kwant/tests/test_qsymm.py b/kwant/tests/test_qsymm.py
index 7ad21aa6..2a13b5e4 100644
--- a/kwant/tests/test_qsymm.py
+++ b/kwant/tests/test_qsymm.py
@@ -404,14 +404,27 @@ def test_find_builder_discrete_symmetries():
         bulk[lat(0, 0)] = h_ons
         bulk[kwant.builder.HoppingKind((1, 0), lat)] = h_hop
         bulk[kwant.builder.HoppingKind((0, 1), lat)] = h_hop
-        builder_symmetries = find_builder_symmetries(bulk, spatial_symmetries=True, prettify=True)
+
+        builder_symmetries_default = find_builder_symmetries(bulk, spatial_symmetries=True,
+                                                             prettify=True)
+        builder_symmetries_sparse = find_builder_symmetries(bulk, spatial_symmetries=True,
+                                                            prettify=True, sparse=True)
+        builder_symmetries_dense = find_builder_symmetries(bulk, spatial_symmetries=True,
+                                                            prettify=True, sparse=False)
+
+        assert len(builder_symmetries_default) == len(builder_symmetries_sparse)
+        assert len(builder_symmetries_default) == len(builder_symmetries_dense)
 
         # Equality of symmetries ignores unitary part
         fourfold_rotation = qsymm.PointGroupElement(np.array([[0, 1],[1, 0]]), False, False, None)
-        assert fourfold_rotation in builder_symmetries
+        assert fourfold_rotation in builder_symmetries_default
+        assert fourfold_rotation in builder_symmetries_sparse
+        assert fourfold_rotation in builder_symmetries_dense
         class_symmetries = class_dict[sym]
         for class_symmetry in class_symmetries:
-            assert sym_dict[class_symmetry] in builder_symmetries
+            assert sym_dict[class_symmetry] in builder_symmetries_default
+            assert sym_dict[class_symmetry] in builder_symmetries_sparse
+            assert sym_dict[class_symmetry] in builder_symmetries_dense
 
 
 def random_onsite_hop(n, rng=0):
-- 
GitLab