add sparse option and a heuristic for using it in symmetry analysis

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)

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):