Add kron dense dia

qutip/core/data/kron.pxd

@@ -10,6 +10,8 @@ cpdef CSR kron_csr(CSR left, CSR right)
 cpdef Dia kron_dia(Dia left, Dia right)
 cpdef CSR kron_dense_csr_csr(Dense left, CSR right)
 cpdef CSR kron_csr_dense_csr(CSR left, Dense right)
+cpdef Dia kron_dia_dense_dia(Dia left, Dense right)
+cpdef Dia kron_dense_dia_dia(Dense left, Dia right)
 
 cpdef Data kron_transpose_data(Data left, Data right)
 cpdef Dense kron_transpose_dense(Dense left, Dense right)

qutip/core/data/kron.pyx

@@ -16,6 +16,7 @@ import numpy
 __all__ = [
     'kron', 'kron_csr', 'kron_dense', 'kron_dia',
     'kron_csr_dense_csr', 'kron_dense_csr_csr',
+    'kron_dia_dense_dia', 'kron_dense_dia_dia',
     'kron_transpose', 'kron_transpose_dense', 'kron_transpose_data',
 ]
 
@@ -176,6 +177,22 @@ cpdef Dia kron_dia(Dia left, Dia right):
     return out
 
 
+cpdef Dia kron_dia_dense_dia(Dia left, Dense right):
+    # The dispatcher would use kron_dense, but the output is at least as sparse
+    # as the sparse input. Since the dispatcher does not have precise control
+    # on which function to use when the signature is missing. We add
+    # then like this.
+    return kron_dia(left, _to(Dia, right))
+
+
+cpdef Dia kron_dense_dia_dia(Dense left, Dia right):
+    # The dispatcher would use kron_dense, but the output is at least as sparse
+    # as the sparse input. Since the dispatcher does not have precise control
+    # on which function to use when the signature is missing. We add
+    # then like this.
+    return kron_dia(_to(Dia, left), right)
+
+
 from .dispatch import Dispatcher as _Dispatcher
 import inspect as _inspect
 
@@ -200,6 +217,8 @@ kron.add_specialisations([
     (Dia, Dia, Dia, kron_dia),
     (CSR, Dense, CSR, kron_csr_dense_csr),
     (Dense, CSR, CSR, kron_dense_csr_csr),
+    (Dia, Dense, Dia, kron_dia_dense_dia),
+    (Dense, Dia, Dia, kron_dense_dia_dia),
 ], _defer=True)
 
 

qutip/core/tensor.py

@@ -467,4 +467,4 @@ def expand_operator(oper, dims, targets, dtype=None):
     for i, ind in enumerate(rest_pos):
         new_order[ind] = rest_qubits[i]
     id_list = [identity(dims[i]) for i in rest_pos]
-    return tensor([oper] + id_list).permute(new_order)
+    return tensor([oper] + id_list).permute(new_order).to(dtype)

qutip/tests/core/data/test_mathematics.py

@@ -764,6 +764,8 @@ def op_numpy(self, left, right):
         pytest.param(data.kron_dia, Dia, Dia, Dia),
         pytest.param(data.kron_dense_csr_csr, Dense, CSR, CSR),
         pytest.param(data.kron_csr_dense_csr, CSR, Dense, CSR),
+        pytest.param(data.kron_dense_dia_dia, Dense, Dia, Dia),
+        pytest.param(data.kron_dia_dense_dia, Dia, Dense, Dia),
     ]