Merge pull request #2008 from Ericgig/FillDense

Add missing Dense specialisations

qutip/core/data/constant.py

@@ -30,8 +30,8 @@
     (which is their representation of 0), and dense matrices will still be
     filled.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     rows, cols : int
         The number of rows and columns in the output matrix.
     """
@@ -58,8 +58,8 @@
     `scale` can be given, where all the diagonal elements will be that instead
     of 1.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     dimension : int
         The dimension of the square output identity matrix.
     scale : complex, optional

qutip/core/data/dispatch.pyx

@@ -347,8 +347,8 @@ cdef class Dispatcher:
         """
         Create a new data layer dispatching operator.
 
-        Arguments
-        ---------
+        Parameters
+        ----------
         signature_source : callable or inspect.Signature
             An object from which the call signature of operation can be
             determined.  You can pass any callable defined in Python space, and
@@ -430,8 +430,8 @@ cdef class Dispatcher:
         recent version; you can use this to override currently known
         specialisations if desired.
 
-        Arguments
-        ---------
+        Parameters
+        ----------
         specialisations : iterable of tuples
             An iterable where each element specifies a new specialisation for
             this operation.  Each element of the iterable should be a tuple,

qutip/core/data/inner.pyx

@@ -5,12 +5,14 @@ cdef extern from "<complex>" namespace "std" nogil:
     double complex conj(double complex x)
 
 from qutip.core.data.base cimport idxint, Data
-from qutip.core.data cimport csr
+from qutip.core.data cimport csr, dense
 from qutip.core.data.csr cimport CSR
+from qutip.core.data.dense cimport Dense
+from qutip.core.data.matmul cimport matmul_dense
 
 __all__ = [
-    'inner', 'inner_csr',
-    'inner_op', 'inner_op_csr',
+    'inner', 'inner_csr', 'inner_dense',
+    'inner_op', 'inner_op_csr', 'inner_op_dense',
 ]
 
 
@@ -88,6 +90,33 @@ cpdef double complex inner_csr(CSR left, CSR right, bint scalar_is_ket=False) no
         return _inner_csr_bra_ket(left, right)
     return _inner_csr_ket_ket(left, right)
 
+cpdef double complex inner_dense(Dense left, Dense right, bint scalar_is_ket=False) nogil except *:
+    """
+    Compute the complex inner product <left|right>.  The shape of `left` is
+    used to determine if it has been supplied as a ket or a bra.  The result of
+    this function will be identical if passed `left` or `adjoint(left)`.
+
+    The parameter `scalar_is_ket` is only intended for the case where `left`
+    and `right` are both of shape (1, 1).  In this case, `left` will be assumed
+    to be a ket unless `scalar_is_ket` is False.  This parameter is ignored at
+    all other times.
+    """
+    _check_shape_inner(left, right)
+    if left.shape[0] == left.shape[1] == right.shape[1] == 1:
+        return (
+                conj(left.data[0]) * right.data[0] if scalar_is_ket
+                else left.data[0] * right.data[0]
+        )
+    cdef double complex out = 0
+    cdef size_t i
+    if left.shape[0] == 1:
+        for i in range(right.shape[0]):
+            out += left.data[i] * right.data[i]
+    else:
+        for i in range(right.shape[0]):
+            out += conj(left.data[i]) * right.data[i]
+    return out
+
 
 cdef double complex _inner_op_csr_bra_ket(CSR left, CSR op, CSR right) nogil:
     cdef size_t ptr_l, ptr_op, ptr_r, row, col
@@ -143,6 +172,21 @@ cpdef double complex inner_op_csr(CSR left, CSR op, CSR right,
         return _inner_op_csr_bra_ket(left, op, right)
     return _inner_op_csr_ket_ket(left, op, right)
 
+cpdef double complex inner_op_dense(Dense left, Dense op, Dense right,
+                                  bint scalar_is_ket=False) except *:
+    """
+    Compute the complex inner product <left|op|right>.  The shape of `left` is
+    used to determine if it has been supplied as a ket or a bra.  The result of
+    this function will be identical if passed `left` or `adjoint(left)`.
+
+    The parameter `scalar_is_ket` is only intended for the case where `left`
+    and `right` are both of shape (1, 1).  In this case, `left` will be assumed
+    to be a ket unless `scalar_is_ket` is False.  This parameter is ignored at
+    all other times.
+    """
+    _check_shape_inner_op(left, op, right)
+    return inner_dense(left, matmul_dense(op, right), scalar_is_ket)
+
 
 from .dispatch import Dispatcher as _Dispatcher
 import inspect as _inspect
@@ -172,8 +216,8 @@ inner.__doc__ =\
     to be a ket unless `scalar_is_ket` is False.  This parameter is ignored at
     all other times.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     left : Data
         The left operand as either a bra or a ket matrix.
 
@@ -187,6 +231,7 @@ inner.__doc__ =\
     """
 inner.add_specialisations([
     (CSR, CSR, inner_csr),
+    (Dense, Dense, inner_dense),
 ], _defer=True)
 
 inner_op = _Dispatcher(
@@ -216,8 +261,8 @@ inner_op.__doc__ =\
     to be a ket unless `scalar_is_ket` is False.  This parameter is ignored at
     all other times.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     left : Data
         The left operand as either a bra or a ket matrix.
 
@@ -235,6 +280,7 @@ inner_op.__doc__ =\
     """
 inner_op.add_specialisations([
     (CSR, CSR, CSR, inner_op_csr),
+    (Dense, Dense, Dense, inner_op_dense),
 ], _defer=True)
 
 del _inspect, _Dispatcher

qutip/core/data/matmul.pyx

@@ -255,6 +255,23 @@ cpdef Dense matmul_dense(Dense left, Dense right, double complex scale=1, Dense
     cdef double complex *b
     cdef char transa, transb
     cdef int m, n, k=left.shape[1], lda, ldb
+    if right.shape[1] == 1:
+        # Matrix Vector product
+        a, b = left.data, right.data
+        if left.fortran:
+            lda = left.shape[0]
+            transa = b'n'
+            m = left.shape[0]
+            n = left.shape[1]
+        else:
+            lda = left.shape[1]
+            transa = b't'
+            m = left.shape[1]
+            n = left.shape[0]
+        ldb = 1
+        blas.zgemv(&transa, &m , &n, &scale, a, &lda, b, &ldb,
+                   &out_scale, out.data, &ldb)
+        return out
     # We use the BLAS routine zgemm for every single call and pretend that
     # we're always supplying it with Fortran-ordered matrices, but to achieve
     # what we want, we use the property of matrix multiplication that
@@ -397,8 +414,8 @@ matmul.__doc__ =\
     where `scale` is (optionally) a scalar, and `left` and `right` are
     matrices.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     left : Data
         The left operand as either a bra or a ket matrix.
 

qutip/core/data/permute.pyx

@@ -14,8 +14,7 @@ import numpy as np
 cimport numpy as cnp
 
 from qutip.core.data.base cimport idxint, idxint_DTYPE
-from qutip.core.data cimport csr
-from qutip.core.data.csr cimport CSR
+from qutip.core.data cimport csr, dense, CSR, Dense
 
 from qutip.core.data.base import idxint_dtype
 
@@ -201,6 +200,16 @@ cpdef CSR indices_csr(CSR matrix, object row_perm=None, object col_perm=None):
                              np.asarray(row_perm, dtype=idxint_dtype),
                              np.asarray(col_perm, dtype=idxint_dtype))
 
+cpdef Dense indices_dense(Dense matrix, object row_perm=None, object col_perm=None):
+    if row_perm is None and col_perm is None:
+        return matrix.copy()
+    array = matrix.as_ndarray()
+    if row_perm is not None:
+        array = array[np.argsort(row_perm), :]
+    if col_perm is not None:
+        array = array[:, np.argsort(col_perm)]
+    return Dense(array)
+
 
 cdef CSR _dimensions_csr_columns(CSR matrix, _Indexer index):
     if matrix.shape[0] != 1:
@@ -293,6 +302,18 @@ cpdef CSR dimensions_csr(CSR matrix, object dimensions, object order):
         return _indices_csr_full(matrix, permutation, permutation)
     return _dimensions_csr_sparse(matrix, index)
 
+@cython.cdivision(True)
+cpdef Dense dimensions_dense(Dense matrix, object dimensions, object order):
+    cdef _Indexer index = _Indexer(np.asarray(dimensions, dtype=idxint_dtype),
+                                   np.asarray(order, dtype=idxint_dtype))
+    cdef idxint[:] permutation = index.all()
+    row_perm, col_perm = None, None
+    if matrix.shape[0] != 1:
+        row_perm = permutation
+    if matrix.shape[1] != 1:
+        col_perm = permutation
+    return indices_dense(matrix, row_perm, col_perm)
+
 
 from .dispatch import Dispatcher as _Dispatcher
 import inspect as _inspect
@@ -323,8 +344,8 @@ dimensions.__doc__ =\
     In other words, the inputs to `kron` are reordered so that input `n` moves
     to position `order[n]`.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     matrix : Data
         Input matrix to reorder.  This can either be a square matrix
         representing an operator, or a bra- or ket-like vector.
@@ -341,6 +362,7 @@ dimensions.__doc__ =\
     """
 dimensions.add_specialisations([
     (CSR, CSR, dimensions_csr),
+    (Dense, Dense, dimensions_dense),
 ], _defer=True)
 
 indices = _Dispatcher(
@@ -363,8 +385,8 @@ indices.__doc__ =\
     of quantum states; if you want to "reorder" the tensor-product structure of
     a system, you want `permute.dimensions` instead.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     matrix : Data
         The input matrix.
 
@@ -377,6 +399,7 @@ indices.__doc__ =\
     """
 indices.add_specialisations([
     (CSR, CSR, indices_csr),
+    (Dense, Dense, indices_dense),
 ], _defer=True)
 
 del _inspect, _Dispatcher

qutip/core/data/pow.pyx

@@ -3,12 +3,14 @@
 
 cimport cython
 
-from qutip.core.data cimport csr
+from qutip.core.data cimport csr, dense
 from qutip.core.data.csr cimport CSR
+from qutip.core.data.dense cimport Dense
 from qutip.core.data.matmul cimport matmul_csr
+import numpy as np
 
 __all__ = [
-    'pow', 'pow_csr',
+    'pow', 'pow_csr', 'pow_dense',
 ]
 
 
@@ -39,6 +41,16 @@ cpdef CSR pow_csr(CSR matrix, unsigned long long n):
     return out
 
 
+cpdef Dense pow_dense(Dense matrix, unsigned long long n):
+    if matrix.shape[0] != matrix.shape[1]:
+        raise ValueError("matrix power only works with square matrices")
+    if n == 0:
+        return dense.identity(matrix.shape[0])
+    if n == 1:
+        return matrix.copy()
+    return Dense(np.linalg.matrix_power(matrix.as_ndarray(), n))
+
+
 from .dispatch import Dispatcher as _Dispatcher
 import inspect as _inspect
 
@@ -58,8 +70,8 @@ pow.__doc__ =\
     must be an integer >= 0.  `A ** 0` is defined to be the identity matrix of
     the same shape.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     matrix : Data
         Input matrix to take the power of.
 
@@ -68,6 +80,7 @@ pow.__doc__ =\
     """
 pow.add_specialisations([
     (CSR, CSR, pow_csr),
+    (Dense, Dense, pow_dense),
 ], _defer=True)
 
 del _inspect, _Dispatcher

qutip/core/data/project.pyx

@@ -131,8 +131,8 @@ project.__doc__ =\
     Get the projector of a state with itself.  Mathematically, if passed an
     object `|a>` or `<a|`, then return the matrix `|a><a|`.
 
-    Arguments
-    ---------
+    Parameters
+    ----------
     state : Data
         The input state bra- or ket-like vector.
     """