Merge branch 'qutip:master' into parity

doc/changes/2257.misc

@@ -0,0 +1 @@
+Improve fidelity doc-string

qutip/core/coefficient.py

@@ -245,6 +245,9 @@ class CompilationOptions(QutipOptions):
 
     clean_on_error : bool [True]
         When writing a cython file that cannot be imported, erase it.
+
+    build_dir: str [None]
+        cythonize's build_dir.
     """
     _link_flags = ""
     _compiler_flags = ""
@@ -275,6 +278,7 @@ class CompilationOptions(QutipOptions):
         "link_flags": _link_flags,
         "extra_import": "",
         "clean_on_error": True,
+        "build_dir": None,
     }
     _settings_name = "compile"
 
@@ -554,7 +558,9 @@ def compile_code(code, file_name, parsed, c_opt):
                 include_dirs=[np.get_include()],
                 language='c++'
             )
-            ext_modules = cythonize(coeff_file, force=True)
+            ext_modules = cythonize(
+                coeff_file, force=True, build_dir=c_opt['build_dir']
+            )
             setup(ext_modules=ext_modules)
         except Exception as e:
             if c_opt['clean_on_error']:

qutip/core/data/csr.pyx

@@ -643,29 +643,26 @@ cdef CSR from_coo_pointers(
 
 
 cpdef CSR from_dia(Dia matrix):
-    if matrix.num_diag == 0:
-        return zeros(*matrix.shape)
-    mat = matrix.as_scipy()
-    ordered = np.argsort(mat.offsets)
-    nnz = len(mat.offsets) * max(mat.shape)
-    ptrs = np.zeros(mat.shape[0]+1, dtype=idxint_dtype)
-    indices = np.zeros(nnz, dtype=idxint_dtype)
-    data = np.zeros(nnz, dtype=complex)
-
-    ptr = 0
-    for row in range(mat.shape[0]):
-        for idx in ordered:
-            off = mat.offsets[idx]
-            if off + row < 0:
+    cdef base.idxint col, diag, i, ptr=0
+    cdef base.idxint nrows=matrix.shape[0], ncols=matrix.shape[1]
+    cdef base.idxint nnz = matrix.num_diag * min(matrix.shape)
+    cdef double complex[:] data = np.zeros(nnz, dtype=complex)
+    cdef base.idxint[:] cols = np.zeros(nnz, dtype=idxint_dtype)
+    cdef base.idxint[:] rows = np.zeros(nnz, dtype=idxint_dtype)
+
+    for i in range(matrix.num_diag):
+        diag = matrix.offsets[i]
+
+        for col in range(ncols):
+            if col - diag < 0 or col - diag >= nrows:
                 continue
-            elif off + row >= mat.shape[1]:
-                break
-            indices[ptr] = off + row
-            data[ptr] = mat.data[idx, off + row]
+            data[ptr] = matrix.data[i * ncols + col]
+            rows[ptr] = col - diag
+            cols[ptr] = col
             ptr += 1
-        ptrs[row + 1] = ptr
 
-    return CSR((data, indices, ptrs), matrix.shape, copy=False)
+    return from_coo_pointers(&rows[0], &cols[0], &data[0], matrix.shape[0],
+                             matrix.shape[1], nnz)
 
 
 cdef inline base.idxint _diagonal_length(

qutip/core/data/solve.py

@@ -75,6 +75,9 @@ def solve_csr_dense(matrix: Union[CSR, Dia], target: Dense, method=None,
         raise ValueError("mkl is not available")
     elif method == "mkl_spsolve":
         solver = mkl_spsolve
+        # mkl does not support dia.
+        if isinstance(matrix, Dia):
+            matrix = _data.to("CSR", matrix)
     else:
         raise ValueError(f"Unknown sparse solver {method}.")
 

qutip/core/metrics.py

@@ -9,13 +9,10 @@
            'hellinger_dist', 'hilbert_dist', 'average_gate_fidelity',
            'process_fidelity', 'unitarity', 'dnorm']
 
-import warnings
-
 import numpy as np
 from scipy import linalg as la
 import scipy.sparse as sp
-from .superop_reps import (to_kraus, to_choi, _to_superpauli, to_super,
-                           kraus_to_choi)
+from .superop_reps import to_choi, _to_superpauli, to_super, kraus_to_choi
 from .superoperator import operator_to_vector, vector_to_operator
 from .operators import qeye, qeye_like
 from .states import ket2dm
@@ -31,7 +28,13 @@
 def fidelity(A, B):
     """
     Calculates the fidelity (pseudo-metric) between two density matrices.
-    See: Nielsen & Chuang, "Quantum Computation and Quantum Information"
+
+    Notes
+    -----
+    Uses the definition from Nielsen & Chuang, "Quantum Computation and Quantum
+    Information". It is the square root of the fidelity defined in
+    R. Jozsa, Journal of Modern Optics, 41:12, 2315 (1994), used in
+    :func:`qutip.core.metrics.process_fidelity`.
 
     Parameters
     ----------
@@ -196,7 +199,7 @@ def process_fidelity(oper, target=None):
         if isinstance(oper, list):  # oper is a list of Kraus operators
             return _process_fidelity_to_id([k * target.dag() for k in oper])
         elif oper.type == 'oper':
-            return _process_fidelity_to_id(oper*target.dag())
+            return _process_fidelity_to_id(oper * target.dag())
         elif oper.type == 'super':
             oper_super = to_super(oper)
             target_dag_super = to_super(target.dag())

qutip/solver/steadystate.py

@@ -126,9 +126,10 @@ def steadystate(A, c_ops=[], *, method='direct', solver=None, **kwargs):
         raise TypeError('Cannot calculate the steady state for a ' +
                         'non-dissipative system.')
     if not A.issuper:
-        A = liouvillian(A)
-    for op in c_ops:
-        A += lindblad_dissipator(op)
+        A = liouvillian(A, c_ops)
+    else:
+        for op in c_ops:
+            A += lindblad_dissipator(op)
 
     if "-" in method:
         # to support v4's "power-gmres" method
@@ -202,10 +203,14 @@ def _steadystate_direct(A, weight, **kw):
     N = A.shape[0]
     n = int(N**0.5)
     dtype = type(A.data)
+    if dtype == _data.Dia:
+        # Dia is bad at vector, the following matmul is 10x slower with Dia
+        # than CSR and Dia is missing optimization such as `use_wbm`.
+        dtype = _data.CSR
     weight_vec = _data.column_stack(_data.diag([weight] * n, 0, dtype=dtype))
-    weight_mat = _data.kron(
-        weight_vec.transpose(),
-        _data.one_element[dtype]((N, 1), (0, 0), 1)
+    weight_mat = _data.matmul(
+        _data.one_element[dtype]((N, 1), (0, 0), 1),
+        weight_vec.transpose()
     )
     L = _data.add(weight_mat, A.data)
     b = _data.one_element[dtype]((N, 1), (0, 0), weight)
@@ -215,14 +220,14 @@ def _steadystate_direct(A, weight, **kw):
         if isinstance(L, _data.CSR):
             L, b = _permute_wbm(L, b)
         else:
-            warn("Only CSR matrice can be permuted.", RuntimeWarning)
+            warn("Only CSR matrices can be permuted.", RuntimeWarning)
     use_rcm = False
     if kw.pop("use_rcm", False):
         if isinstance(L, _data.CSR):
             L, b, perm = _permute_rcm(L, b)
             use_rcm = True
         else:
-            warn("Only CSR matrice can be permuted.", RuntimeWarning)
+            warn("Only CSR matrices can be permuted.", RuntimeWarning)
     if kw.pop("use_precond", False):
         if isinstance(L, (_data.CSR, _data.Dia)):
             kw["M"] = _compute_precond(L, kw)
@@ -271,14 +276,14 @@ def _steadystate_power(A, **kw):
         if isinstance(L, _data.CSR):
             L, y = _permute_wbm(L, y)
         else:
-            warn("Only CSR matrice can be permuted.", RuntimeWarning)
+            warn("Only CSR matrices can be permuted.", RuntimeWarning)
     use_rcm = False
     if kw.pop("use_rcm", False):
         if isinstance(L, _data.CSR):
             L, y, perm = _permute_rcm(L, y)
             use_rcm = True
         else:
-            warn("Only CSR matrice can be permuted.", RuntimeWarning)
+            warn("Only CSR matrices can be permuted.", RuntimeWarning)
     if kw.pop("use_precond", False):
         if isinstance(L, (_data.CSR, _data.Dia)):
             kw["M"] = _compute_precond(L, kw)

qutip/tests/solver/test_steadystate.py

@@ -29,10 +29,12 @@
     pytest.param('iterative-bicgstab', {'atol': 1e-12, "tol": 1e-10},
                  id="iterative-bicgstab"),
 ])
-def test_qubit(method, kwargs):
+@pytest.mark.parametrize("dtype", ["dense", "dia", "csr"])
+@pytest.mark.filterwarnings("ignore:Only CSR matrices:RuntimeWarning")
+def test_qubit(method, kwargs, dtype):
     # thermal steadystate of a qubit: compare numerics with analytical formula
-    sz = qutip.sigmaz()
-    sm = qutip.destroy(2)
+    sz = qutip.sigmaz().to(dtype)
+    sm = qutip.destroy(2, dtype=dtype)
 
     H = 0.5 * 2 * np.pi * sz
     gamma1 = 0.05