Python Qobj.dims方法代码示例

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def _steadystate_power(L, maxiter=10, tol=1e-6, itertol=1e-5,
                       verbose=False):
    """
    Inverse power method for steady state solving.
    """
    if verbose:
        print('Starting iterative power method Solver...')
    use_solver(assumeSortedIndices=True)
    rhoss = Qobj()
    sflag = issuper(L)
    if sflag:
        rhoss.dims = L.dims[0]
        rhoss.shape = [prod(rhoss.dims[0]), prod(rhoss.dims[1])]
    else:
        rhoss.dims = [L.dims[0], 1]
        rhoss.shape = [prod(rhoss.dims[0]), 1]
    n = prod(rhoss.shape)
    L = L.data.tocsc() - (tol ** 2) * sp.eye(n, n, format='csc')
    L.sort_indices()
    v = mat2vec(rand_dm(rhoss.shape[0], 0.5 / rhoss.shape[0] + 0.5).full())
    if verbose:
        start_time = time.time()
    it = 0
    while (la.norm(L * v, np.inf) > tol) and (it < maxiter):
        v = spsolve(L, v)
        v = v / la.norm(v, np.inf)
        it += 1
    if it >= maxiter:
        raise Exception('Failed to find steady state after ' +
                        str(maxiter) + ' iterations')
    # normalise according to type of problem
    if sflag:
        trow = sp.eye(rhoss.shape[0], rhoss.shape[0], format='coo')
        trow = sp_reshape(trow, (1, n))
        data = v / sum(trow.dot(v))
    else:
        data = data / la.norm(v)

    data = sp.csr_matrix(vec2mat(data))
    rhoss.data = 0.5 * (data + data.conj().T)
    rhoss.isherm = True
    if verbose:
        print('Power solver time: ', time.time() - start_time)
    if qset.auto_tidyup:
        return rhoss.tidyup()
    else:
        return rhoss

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def operator_to_vector(op):
    """
    Create a vector representation of a quantum operator given
    the matrix representation.
    """
    q = Qobj()
    q.dims = [op.dims, [1]]
    q.data = sp_reshape(op.data.T, (np.prod(op.shape), 1))
    return q

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def vector_to_operator(op):
    """
    Create a matrix representation given a quantum operator in
    vector form.
    """
    q = Qobj()
    q.dims = op.dims[0]
    q.data = sp_reshape(op.data.T, q.shape).T
    return q

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def vector_to_operator(op):
    """
    Create a matrix representation given a quantum operator in
    vector form.
    """
    q = Qobj()
    q.dims = op.dims[0]
    n = int(np.sqrt(op.shape[0]))
    q.data = sp_reshape(op.data.T, (n, n)).T
    return q

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def _steadystate_power(L, ss_args):
    """
    Inverse power method for steady state solving.
    """
    if settings.debug:
        print('Starting iterative power method Solver...')
    tol=ss_args['tol']
    maxiter=ss_args['maxiter']
    use_solver(assumeSortedIndices=True)
    rhoss = Qobj()
    sflag = issuper(L)
    if sflag:
        rhoss.dims = L.dims[0]
    else:
        rhoss.dims = [L.dims[0], 1]
    n = prod(rhoss.shape)
    L = L.data.tocsc() - (tol ** 2) * sp.eye(n, n, format='csc')
    L.sort_indices()
    v = mat2vec(rand_dm(rhoss.shape[0], 0.5 / rhoss.shape[0] + 0.5).full())

    it = 0
    while (la.norm(L * v, np.inf) > tol) and (it < maxiter):
        v = spsolve(L, v)
        v = v / la.norm(v, np.inf)
        it += 1
    if it >= maxiter:
        raise Exception('Failed to find steady state after ' +
                        str(maxiter) + ' iterations')
    # normalise according to type of problem
    if sflag:
        trow = sp.eye(rhoss.shape[0], rhoss.shape[0], format='coo')
        trow = sp_reshape(trow, (1, n))
        data = v / sum(trow.dot(v))
    else:
        data = data / la.norm(v)

    data = sp.csr_matrix(vec2mat(data))
    rhoss.data = 0.5 * (data + data.conj().T)
    rhoss.isherm = True

    return rhoss

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def rand_unitary_haar(N=2, dims=None):
    """
    Returns a Haar random unitary matrix of dimension
    ``dim``, using the algorithm of [Mez07]_.

    Parameters
    ----------
    N : int
        Dimension of the unitary to be returned.
    dims : list of lists of int, or None
        Dimensions of quantum object.  Used for specifying
        tensor structure. Default is dims=[[N],[N]].

    Returns
    -------
    U : Qobj
        Unitary of dims ``[[dim], [dim]]`` drawn from the Haar
        measure.
    """
    if dims is not None:
        _check_dims(dims, N, N)
    else:
        dims = [[N], [N]]

    # Mez01 STEP 1: Generate an N × N matrix Z of complex standard
    #               normal random variates.
    Z = randnz((N, N))

    # Mez01 STEP 2: Find a QR decomposition Z = Q · R.
    Q, R = la.qr(Z)

    # Mez01 STEP 3: Create a diagonal matrix Lambda by rescaling
    #               the diagonal elements of R.
    Lambda = np.diag(R).copy()
    Lambda /= np.abs(Lambda)

    # Mez01 STEP 4: Note that R' := Λ¯¹ · R has real and
    #               strictly positive elements, such that
    #               Q' = Q · Λ is Haar distributed.
    # NOTE: Λ is a diagonal matrix, represented as a vector
    #       of the diagonal entries. Thus, the matrix dot product
    #       is represented nicely by the NumPy broadcasting of
    #       the *scalar* multiplication. In particular,
    #       Q · Λ = Q_ij Λ_jk = Q_ij δ_jk λ_k = Q_ij λ_j.
    #       As NumPy arrays, Q has shape (N, N) and
    #       Lambda has shape (N, ), such that the broadcasting
    #       represents precisely Q_ij λ_j.
    U = Qobj(Q * Lambda)
    U.dims = dims
    return U

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def spre(A):
    """Superoperator formed from pre-multiplication by operator A.

    Parameters
    ----------
    A : qobj
        Quantum operator for pre-multiplication.

    Returns
    --------
    super :qobj
        Superoperator formed from input quantum object.
    """
    if not isinstance(A, Qobj):
        raise TypeError('Input is not a quantum object')

    if not A.isoper:
        raise TypeError('Input is not a quantum operator')

    S = Qobj(isherm=A.isherm, superrep='super')
    S.dims = [[A.dims[0], A.dims[1]], [A.dims[0], A.dims[1]]]
    S.data = sp.kron(sp.identity(np.prod(A.shape[1])), A.data, format='csr')
    return S

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def spost(A):
    """Superoperator formed from post-multiplication by operator A

    Parameters
    ----------
    A : qobj
        Quantum operator for post multiplication.

    Returns
    -------
    super : qobj
        Superoperator formed from input qauntum object.
    """
    if not isinstance(A, Qobj):
        raise TypeError('Input is not a quantum object')

    if not A.isoper:
        raise TypeError('Input is not a quantum operator')

    S = Qobj(isherm=A.isherm, superrep='super')
    S.dims = [[A.dims[0], A.dims[1]], [A.dims[0], A.dims[1]]]
    S.data = zcsr_kron(A.data.T, 
                fast_identity(np.prod(A.shape[0])))
    return S

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def liouvillian(H, c_ops=[], data_only=False, chi=None):
    """Assembles the Liouvillian superoperator from a Hamiltonian
    and a ``list`` of collapse operators. Like liouvillian, but with an
    experimental implementation which avoids creating extra Qobj instances,
    which can be advantageous for large systems.

    Parameters
    ----------
    H : qobj
        System Hamiltonian.

    c_ops : array_like
        A ``list`` or ``array`` of collapse operators.

    Returns
    -------
    L : qobj
        Liouvillian superoperator.

    """

    if chi and len(chi) != len(c_ops):
        raise ValueError('chi must be a list with same length as c_ops')

    if H is not None:
        if H.isoper:
            op_dims = H.dims
            op_shape = H.shape
        elif H.issuper:
            op_dims = H.dims[0]
            op_shape = [np.prod(op_dims[0]), np.prod(op_dims[0])]
        else:
            raise TypeError("Invalid type for Hamiltonian.")
    else:
        # no hamiltonian given, pick system size from a collapse operator
        if isinstance(c_ops, list) and len(c_ops) > 0:
            c = c_ops[0]
            if c.isoper:
                op_dims = c.dims
                op_shape = c.shape
            elif c.issuper:
                op_dims = c.dims[0]
                op_shape = [np.prod(op_dims[0]), np.prod(op_dims[0])]
            else:
                raise TypeError("Invalid type for collapse operator.")
        else:
            raise TypeError("Either H or c_ops must be given.")

    sop_dims = [[op_dims[0], op_dims[0]], [op_dims[1], op_dims[1]]]
    sop_shape = [np.prod(op_dims), np.prod(op_dims)]

    spI = sp.identity(op_shape[0], format='csr')

    if H:
        if H.isoper:
            data = -1j * (sp.kron(spI, H.data, format='csr')
                          - sp.kron(H.data.T, spI, format='csr'))
        else:
            data = H.data
    else:
        data = sp.csr_matrix((sop_shape[0], sop_shape[1]), dtype=complex)

    for idx, c_op in enumerate(c_ops):
        if c_op.issuper:
            data = data + c_op.data
        else:
            cd = c_op.data.T.conj()
            c = c_op.data
            if chi:
                data = data + np.exp(1j * chi[idx]) * sp.kron(cd.T, c,
                                                              format='csr')
            else:
                data = data + sp.kron(cd.T, c, format='csr')
            cdc = cd * c
            data = data - 0.5 * sp.kron(spI, cdc, format='csr')
            data = data - 0.5 * sp.kron(cdc.T, spI, format='csr')

    if data_only:
        return data
    else:
        L = Qobj()
        L.dims = sop_dims
        L.data = data
        L.isherm = False
        L.superrep = 'super'
        return L

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def rand_super_bcsz(N=2, enforce_tp=True, rank=None, dims=None):
    """
    Returns a random superoperator drawn from the Bruzda
    et al ensemble for CPTP maps [BCSZ08]_. Note that due to
    finite numerical precision, for ranks less than full-rank,
    zero eigenvalues may become slightly negative, such that the
    returned operator is not actually completely positive.


    Parameters
    ----------
    N : int
        Square root of the dimension of the superoperator to be returned.
    enforce_tp : bool
        If True, the trace-preserving condition of [BCSZ08]_ is enforced;
        otherwise only complete positivity is enforced.
    rank : int or None
        Rank of the sampled superoperator. If None, a full-rank
        superoperator is generated.
    dims : list
        Dimensions of quantum object.  Used for specifying
        tensor structure. Default is dims=[[[N],[N]], [[N],[N]]].

    Returns
    -------
    rho : Qobj
        A superoperator acting on vectorized dim × dim density operators,
        sampled from the BCSZ distribution.
    """
    if dims is not None:
        # TODO: check!
        pass
    else:
        dims = [[[N], [N]], [[N], [N]]]

    if rank is None:
        rank = N ** 2
    if rank > N ** 2:
        raise ValueError("Rank cannot exceed superoperator dimension.")

    # We use mainly dense matrices here for speed in low
    # dimensions. In the future, it would likely be better to switch off
    # between sparse and dense matrices as the dimension grows.

    # We start with a Ginibre uniform matrix X of the appropriate rank,
    # and use it to construct a positive semidefinite matrix X X⁺.
    X = randnz((N ** 2, rank), norm="ginibre")

    # Precompute X X⁺, as we'll need it in two different places.
    XXdag = np.dot(X, X.T.conj())

    if enforce_tp:
        # We do the partial trace over the first index by using dense reshape
        # operations, so that we can avoid bouncing to a sparse representation
        # and back.
        Y = np.einsum("ijik->jk", XXdag.reshape((N, N, N, N)))

        # Now we have the matrix 𝟙 ⊗ Y^{-1/2}, which we can find by doing
        # the square root and the inverse separately. As a possible improvement,
        # iterative methods exist to find inverse square root matrices directly,
        # as this is important in statistics.
        Z = np.kron(np.eye(N), sqrtm(la.inv(Y)))

        # Finally, we dot everything together and pack it into a Qobj,
        # marking the dimensions as that of a type=super (that is,
        # with left and right compound indices, each representing
        # left and right indices on the underlying Hilbert space).
        D = Qobj(np.dot(Z, np.dot(XXdag, Z)))
    else:
        D = N * Qobj(XXdag / np.trace(XXdag))

    D.dims = [
        # Left dims
        [[N], [N]],
        # Right dims
        [[N], [N]],
    ]

    # Since [BCSZ08] gives a row-stacking Choi matrix, but QuTiP
    # expects a column-stacking Choi matrix, we must permute the indices.
    D = D.permute([[1], [0]])

    D.dims = dims

    # Mark that we've made a Choi matrix.
    D.superrep = "choi"

    return sr.to_super(D)

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def tensor(*args):
    """Calculates the tensor product of input operators. 
    
    Parameters
    ----------
    args : array_like 
        ``list`` or ``array`` of quantum objects for tensor product.
        
    Returns
    --------
    obj : qobj
        A composite quantum object.
    
    Examples
    --------    
    >>> tensor([sigmax(), sigmax()])
    Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isHerm = True
    Qobj data = 
    [[ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]
     [ 0.+0.j  0.+0.j  1.+0.j  0.+0.j]
     [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]
     [ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]]

    """
    if not args:
        raise TypeError("Requires at least one input argument")
    num_args=len(args)
    step=0
    for n in range(num_args):
        if isinstance(args[n],Qobj):
            qos=args[n]
            if step==0:
                dat=qos.data
                dim=qos.dims
                shp=qos.shape
                step=1
            else:
                dat=sp.kron(dat,qos.data, format='csr') #sparse Kronecker product
                dim=[dim[0]+qos.dims[0],dim[1]+qos.dims[1]] #append dimensions of Qobjs
                shp=[dat.shape[0],dat.shape[1]] #new shape of matrix
                
        elif isinstance(args[n],(list,ndarray)):#checks if input is list/array of Qobjs
            qos=args[n]
            items=len(qos) #number of inputs
            if not all([isinstance(k,Qobj) for k in qos]): #raise error if one of the inputs is not a quantum object
                raise TypeError("One of inputs is not a quantum object")
            if items==1:# if only one Qobj, do nothing
                if step==0: 
                    dat=qos[0].data
                    dim=qos[0].dims
                    shp=qos[0].shape
                    step=1
                else:
                    dat=sp.kron(dat,qos[0].data, format='csr') #sparse Kronecker product
                    dim=[dim[0]+qos[0].dims[0],dim[1]+qos[0].dims[1]] #append dimensions of Qobjs
                    shp=[dat.shape[0],dat.shape[1]] #new shape of matrix
            elif items!=1:
                if step==0:
                    dat=qos[0].data
                    dim=qos[0].dims
                    shp=qos[0].shape
                    step=1
                for k in range(items-1): #cycle over all items
                    dat=sp.kron(dat,qos[k+1].data, format='csr') #sparse Kronecker product
                    dim=[dim[0]+qos[k+1].dims[0],dim[1]+qos[k+1].dims[1]] #append dimensions of Qobjs
                    shp=[dat.shape[0],dat.shape[1]] #new shape of matrix
    out=Qobj()
    out.data=dat
    out.dims=dim
    out.shape=shp
    if qutip.settings.auto_tidyup:
        return Qobj(out).tidyup() #returns tidy Qobj
    else:
        return Qobj(out)

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def _steadystate_power(L, ss_args):
    """
    Inverse power method for steady state solving.
    """
    ss_args['info'].pop('weight', None)
    if settings.debug:
        logger.debug('Starting iterative inverse-power method solver.')
    tol = ss_args['tol']
    mtol = ss_args['mtol']
    if mtol is None:
        mtol = max(0.1*tol, 1e-15)
    maxiter = ss_args['maxiter']

    use_solver(assumeSortedIndices=True)
    rhoss = Qobj()
    sflag = issuper(L)
    if sflag:
        rhoss.dims = L.dims[0]
    else:
        rhoss.dims = [L.dims[0], 1]
    n = L.shape[0]
    # Build Liouvillian
    if ss_args['solver'] == 'mkl' and ss_args['method'] == 'power':
        has_mkl = 1
    else:
        has_mkl = 0
    L, perm, perm2, rev_perm, ss_args = _steadystate_power_liouvillian(L,
                                                                       ss_args,
                                                                       has_mkl)
    orig_nnz = L.nnz
    # start with all ones as RHS
    v = np.ones(n, dtype=complex)
    if ss_args['use_rcm']:
        v = v[np.ix_(perm2,)]

    # Do preconditioning
    if ss_args['solver'] == 'scipy':
        if ss_args['M'] is None and ss_args['use_precond'] and \
                ss_args['method'] in ['power-gmres',
                                      'power-lgmres',
                                      'power-bicgstab']:
            ss_args['M'], ss_args = _iterative_precondition(L, int(np.sqrt(n)),
                                                            ss_args)
            if ss_args['M'] is None:
                warnings.warn("Preconditioning failed. Continuing without.",
                              UserWarning)

    ss_iters = {'iter': 0}

    def _iter_count(r):
        ss_iters['iter'] += 1
        return

    _power_start = time.time()
    # Get LU factors
    if ss_args['method'] == 'power':
        if ss_args['solver'] == 'mkl':
            lu = mkl_splu(L, max_iter_refine=ss_args['max_iter_refine'],
                          scaling_vectors=ss_args['scaling_vectors'],
                          weighted_matching=ss_args['weighted_matching'])
        else:
            lu = splu(L, permc_spec=ss_args['permc_spec'],
                      diag_pivot_thresh=ss_args['diag_pivot_thresh'],
                      options=dict(ILU_MILU=ss_args['ILU_MILU']))

            if settings.debug and _scipy_check:
                L_nnz = lu.L.nnz
                U_nnz = lu.U.nnz
                logger.debug('L NNZ: %i ; U NNZ: %i' % (L_nnz, U_nnz))
                logger.debug('Fill factor: %f' % ((L_nnz+U_nnz)/orig_nnz))

    it = 0
    # FIXME: These atol keyword except checks can be removed once scipy 1.1
    # is a minimum requirement
    while (la.norm(L * v, np.inf) > tol) and (it < maxiter):
        check = 0
        if ss_args['method'] == 'power':
            v = lu.solve(v)
        elif ss_args['method'] == 'power-gmres':
            try:
                v, check = gmres(L, v, tol=mtol, atol=ss_args['matol'],
                                 M=ss_args['M'], x0=ss_args['x0'],
                                 restart=ss_args['restart'],
                                 maxiter=ss_args['maxiter'],
                                 callback=_iter_count)
            except TypeError as e:
                if "unexpected keyword argument 'atol'" in str(e):
                    v, check = gmres(L, v, tol=mtol,
                                     M=ss_args['M'], x0=ss_args['x0'],
                                     restart=ss_args['restart'],
                                     maxiter=ss_args['maxiter'],
                                     callback=_iter_count)

        elif ss_args['method'] == 'power-lgmres':
            try:
                v, check = lgmres(L, v, tol=mtol, atol=ss_args['matol'],
                                  M=ss_args['M'], x0=ss_args['x0'],
                                  maxiter=ss_args['maxiter'],
                                  callback=_iter_count)
            except TypeError as e:
#.........这里部分代码省略.........

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def _steadystate_power(L, ss_args):
    """
    Inverse power method for steady state solving.
    """
    ss_args['info'].pop('weight', None)
    if settings.debug:
        print('Starting iterative inverse-power method solver...')
    tol = ss_args['tol']
    maxiter = ss_args['maxiter']

    use_solver(assumeSortedIndices=True)
    rhoss = Qobj()
    sflag = issuper(L)
    if sflag:
        rhoss.dims = L.dims[0]
    else:
        rhoss.dims = [L.dims[0], 1]
    n = prod(rhoss.shape)
    L = L.data.tocsc() - (1e-15) * sp.eye(n, n, format='csc')
    L.sort_indices()
    orig_nnz = L.nnz

    # start with all ones as RHS
    v = np.ones(n,dtype=complex)

    if ss_args['use_rcm']:
        if settings.debug:
            old_band = sp_bandwidth(L)[0]
            print('Original bandwidth:', old_band)
        perm = reverse_cuthill_mckee(L)
        rev_perm = np.argsort(perm)
        L = sp_permute(L, perm, perm, 'csc')
        v = v[np.ix_(perm,)]
        if settings.debug:
            new_band = sp_bandwidth(L)[0]
            print('RCM bandwidth:', new_band)
            print('Bandwidth reduction factor:', round(old_band/new_band, 2))

    # Get LU factors
    lu = splu(L, permc_spec=ss_args['permc_spec'],
              diag_pivot_thresh=ss_args['diag_pivot_thresh'],
              options=dict(ILU_MILU=ss_args['ILU_MILU']))

    if settings.debug and _scipy_check:
        L_nnz = lu.L.nnz
        U_nnz = lu.U.nnz
        print('L NNZ:', L_nnz, ';', 'U NNZ:', U_nnz)
        print('Fill factor:', (L_nnz+U_nnz)/orig_nnz)

    _power_start = time.time()
    it = 0
    while (la.norm(L * v, np.inf) > tol) and (it < maxiter):
        v = lu.solve(v)
        v = v / la.norm(v, np.inf)
        it += 1
    if it >= maxiter:
        raise Exception('Failed to find steady state after ' +
                        str(maxiter) + ' iterations')
    
    _power_end = time.time()
    ss_args['info']['solution_time'] = _power_end-_power_start
    ss_args['info']['iterations'] = it
    if settings.debug:
        print('Number of iterations:', it)

    if ss_args['use_rcm']:
        v = v[np.ix_(rev_perm,)]

    # normalise according to type of problem
    if sflag:
        trow = sp.eye(rhoss.shape[0], rhoss.shape[0], format='coo')
        trow = sp_reshape(trow, (1, n))
        data = v / sum(trow.dot(v))
    else:
        data = data / la.norm(v)

    data = sp.csr_matrix(vec2mat(data))
    rhoss.data = 0.5 * (data + data.conj().T)
    rhoss.isherm = True
    if ss_args['return_info']:
        return rhoss, ss_args['info']
    else:
        return rhoss

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def _steadystate_power(L, ss_args):
    """
    Inverse power method for steady state solving.
    """
    ss_args['info'].pop('weight', None)
    if settings.debug:
        logger.debug('Starting iterative inverse-power method solver.')
    tol = ss_args['tol']
    maxiter = ss_args['maxiter']

    use_solver(assumeSortedIndices=True)
    rhoss = Qobj()
    sflag = issuper(L)
    if sflag:
        rhoss.dims = L.dims[0]
    else:
        rhoss.dims = [L.dims[0], 1]
    n = L.shape[0]
    # Build Liouvillian
    if settings.has_mkl and ss_args['method'] == 'power':
        has_mkl = 1
    else:
        has_mkl = 0
    L, perm, perm2, rev_perm, ss_args = _steadystate_power_liouvillian(L, 
                                                ss_args, has_mkl)
    orig_nnz = L.nnz
    # start with all ones as RHS
    v = np.ones(n, dtype=complex)
    if ss_args['use_rcm']:
        v = v[np.ix_(perm2,)]
    
    # Do preconditioning
    if ss_args['M'] is None and ss_args['use_precond'] and \
            ss_args['method'] in ['power-gmres', 
                                'power-lgmres', 'power-bicgstab']:
        ss_args['M'], ss_args = _iterative_precondition(L, int(np.sqrt(n)), ss_args)
        if ss_args['M'] is None:
            warnings.warn("Preconditioning failed. Continuing without.",
                          UserWarning)
    
    ss_iters = {'iter': 0}

    def _iter_count(r):
        ss_iters['iter'] += 1
        return
    
    _power_start = time.time()
    # Get LU factors
    if ss_args['method'] == 'power':
        if settings.has_mkl:
            lu = mkl_splu(L)
        else: 
            lu = splu(L, permc_spec=ss_args['permc_spec'],
              diag_pivot_thresh=ss_args['diag_pivot_thresh'],
              options=dict(ILU_MILU=ss_args['ILU_MILU']))

            if settings.debug and _scipy_check:
                L_nnz = lu.L.nnz
                U_nnz = lu.U.nnz
                logger.debug('L NNZ: %i ; U NNZ: %i' % (L_nnz, U_nnz))
                logger.debug('Fill factor: %f' % ((L_nnz+U_nnz)/orig_nnz))

    it = 0
    _tol = max(ss_args['tol']/10, 1e-15) # Should make this user accessible
    while (la.norm(L * v, np.inf) > tol) and (it < maxiter):
        
        if ss_args['method'] == 'power':
            v = lu.solve(v)
        elif ss_args['method'] == 'power-gmres':
            v, check = gmres(L, v, tol=_tol, M=ss_args['M'],
                                x0=ss_args['x0'], restart=ss_args['restart'],
                                maxiter=ss_args['maxiter'], callback=_iter_count)
        elif ss_args['method'] == 'power-lgmres':
            v, check = lgmres(L, v, tol=_tol, M=ss_args['M'],
                              x0=ss_args['x0'], maxiter=ss_args['maxiter'],
                              callback=_iter_count)
        elif ss_args['method'] == 'power-bicgstab':
            v, check = bicgstab(L, v, tol=_tol, M=ss_args['M'],
                                x0=ss_args['x0'],
                                maxiter=ss_args['maxiter'], callback=_iter_count)
        else:
            raise Exception("Invalid iterative solver method.")
            
        v = v / la.norm(v, np.inf)
        it += 1
    if ss_args['method'] == 'power' and settings.has_mkl:
        lu.delete()
    if it >= maxiter:
        raise Exception('Failed to find steady state after ' +
                        str(maxiter) + ' iterations')
    _power_end = time.time()
    ss_args['info']['solution_time'] = _power_end-_power_start
    ss_args['info']['iterations'] = it
    if ss_args['return_info']:
        ss_args['info']['residual_norm'] = la.norm(L*v)
    if settings.debug:
        logger.debug('Number of iterations: %i' % it)

    if ss_args['use_rcm']:
        v = v[np.ix_(rev_perm,)]
#.........这里部分代码省略.........

# 需要导入模块: from qutip.qobj import Qobj [as 别名]
# 或者: from qutip.qobj.Qobj import dims [as 别名]
def tensor(*args):
    """Calculates the tensor product of input operators.

    Parameters
    ----------
    args : array_like
        ``list`` or ``array`` of quantum objects for tensor product.

    Returns
    -------
    obj : qobj
        A composite quantum object.

    Examples
    --------
    >>> tensor([sigmax(), sigmax()])
    Quantum object: dims = [[2, 2], [2, 2]], \
shape = [4, 4], type = oper, isHerm = True
    Qobj data =
    [[ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]
     [ 0.+0.j  0.+0.j  1.+0.j  0.+0.j]
     [ 0.+0.j  1.+0.j  0.+0.j  0.+0.j]
     [ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]]
    """

    if not args:
        raise TypeError("Requires at least one input argument")

    if len(args) == 1 and isinstance(args[0], (list, np.ndarray)):
        # this is the case when tensor is called on the form:
        # tensor([q1, q2, q3, ...])
        qlist = args[0]

    elif len(args) == 1 and isinstance(args[0], Qobj):
        # tensor is called with a single Qobj as an argument, do nothing
        return args[0]

    else:
        # this is the case when tensor is called on the form:
        # tensor(q1, q2, q3, ...)
        qlist = args

    if not all([isinstance(q, Qobj) for q in qlist]):
        # raise error if one of the inputs is not a quantum object
        raise TypeError("One of inputs is not a quantum object")

    out = Qobj()

    if qlist[0].issuper:
        out.superrep = qlist[0].superrep
        if not all([q.superrep == out.superrep for q in qlist]):
            raise TypeError("In tensor products of superroperators, all must" +
                            "have the same representation")

    out.isherm = True
    for n, q in enumerate(qlist):
        if n == 0:
            out.data = q.data
            out.dims = q.dims
        else:
            out.data = sp.kron(out.data, q.data, format='csr')
            out.dims = [out.dims[0] + q.dims[0], out.dims[1] + q.dims[1]]

        out.isherm = out.isherm and q.isherm

    if not out.isherm:
        out._isherm = None

    return out.tidyup() if qutip.settings.auto_tidyup else out