本文整理汇总了Python中qutip.qobj.Qobj.dims方法的典型用法代码示例。如果您正苦于以下问题:Python Qobj.dims方法的具体用法?Python Qobj.dims怎么用?Python Qobj.dims使用的例子?那么恭喜您, 这里精选的方法代码示例或许可以为您提供帮助。您也可以进一步了解该方法所在类qutip.qobj.Qobj
的用法示例。
在下文中一共展示了Qobj.dims方法的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: _steadystate_power
# 需要导入模块: 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
示例2: operator_to_vector
# 需要导入模块: 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
示例3: vector_to_operator
# 需要导入模块: 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
示例4: vector_to_operator
# 需要导入模块: 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
示例5: _steadystate_power
# 需要导入模块: 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
示例6: rand_unitary_haar
# 需要导入模块: 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
示例7: spre
# 需要导入模块: 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
示例8: spost
# 需要导入模块: 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
示例9: liouvillian
# 需要导入模块: 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
示例10: rand_super_bcsz
# 需要导入模块: 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)
示例11: tensor
# 需要导入模块: 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)
示例12: _steadystate_power
# 需要导入模块: 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:
#.........这里部分代码省略.........
示例13: _steadystate_power
# 需要导入模块: 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
示例14: _steadystate_power
# 需要导入模块: 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,)]
#.........这里部分代码省略.........
示例15: tensor
# 需要导入模块: 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