本文整理汇总了Python中qutip.qobj.issuper函数的典型用法代码示例。如果您正苦于以下问题:Python issuper函数的具体用法?Python issuper怎么用?Python issuper使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了issuper函数的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: _mesolve_const
def _mesolve_const(H, rho0, tlist, c_op_list, e_ops, args, opt,
progress_bar):
"""
Evolve the density matrix using an ODE solver, for constant hamiltonian
and collapse operators.
"""
if debug:
print(inspect.stack()[0][3])
#
# check initial state
#
if isket(rho0):
# if initial state is a ket and no collapse operator where given,
# fall back on the unitary schrodinger equation solver
if len(c_op_list) == 0 and isoper(H):
return _sesolve_const(H, rho0, tlist, e_ops, args, opt,
progress_bar)
# Got a wave function as initial state: convert to density matrix.
rho0 = ket2dm(rho0)
#
# construct liouvillian
#
if opt.tidy:
H = H.tidyup(opt.atol)
L = liouvillian(H, c_op_list)
#
# setup integrator
#
initial_vector = mat2vec(rho0.full()).ravel('F')
if issuper(rho0):
r = scipy.integrate.ode(_ode_super_func)
r.set_f_params(L.data)
else:
if opt.use_openmp and L.data.nnz >= qset.openmp_thresh:
r = scipy.integrate.ode(cy_ode_rhs_openmp)
r.set_f_params(L.data.data, L.data.indices, L.data.indptr,
opt.openmp_threads)
else:
r = scipy.integrate.ode(cy_ode_rhs)
r.set_f_params(L.data.data, L.data.indices, L.data.indptr)
# r = scipy.integrate.ode(_ode_rho_test)
# r.set_f_params(L.data)
r.set_integrator('zvode', method=opt.method, order=opt.order,
atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
first_step=opt.first_step, min_step=opt.min_step,
max_step=opt.max_step)
r.set_initial_value(initial_vector, tlist[0])
#
# call generic ODE code
#
return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)
示例2: _steadystate_setup
def _steadystate_setup(A, c_op_list):
"""Build Liouvillian (if necessary) and check input.
"""
if isoper(A):
if len(c_op_list) > 0:
return liouvillian(A, c_op_list)
raise TypeError('Cannot calculate the steady state for a ' +
'non-dissipative system ' +
'(no collapse operators given)')
elif issuper(A):
return A
else:
raise TypeError('Solving for steady states requires ' +
'Liouvillian (super) operators')
示例3: _steadystate_power
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
示例4: _steadystate_power
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
示例5: _spectrum_pi
def _spectrum_pi(H, wlist, c_ops, a_op, b_op, use_pinv=False):
"""
Internal function for calculating the spectrum of the correlation function
:math:`\left<A(\\tau)B(0)\\right>`.
"""
L = H if issuper(H) else liouvillian(H, c_ops)
tr_mat = tensor([qeye(n) for n in L.dims[0][0]])
N = np.prod(L.dims[0][0])
A = L.full()
b = spre(b_op).full()
a = spre(a_op).full()
tr_vec = np.transpose(mat2vec(tr_mat.full()))
rho_ss = steadystate(L)
rho = np.transpose(mat2vec(rho_ss.full()))
I = np.identity(N * N)
P = np.kron(np.transpose(rho), tr_vec)
Q = I - P
spectrum = np.zeros(len(wlist))
for idx, w in enumerate(wlist):
if use_pinv:
MMR = np.linalg.pinv(-1.0j * w * I + A)
else:
MMR = np.dot(Q, np.linalg.solve(-1.0j * w * I + A, Q))
s = np.dot(tr_vec,
np.dot(a, np.dot(MMR, np.dot(b, np.transpose(rho)))))
spectrum[idx] = -2 * np.real(s[0, 0])
return spectrum
示例6: ode2es
def ode2es(L, rho0):
"""Creates an exponential series that describes the time evolution for the
initial density matrix (or state vector) `rho0`, given the Liouvillian
(or Hamiltonian) `L`.
Parameters
----------
L : qobj
Liouvillian of the system.
rho0 : qobj
Initial state vector or density matrix.
Returns
-------
eseries : :class:`qutip.eseries`
``eseries`` represention of the system dynamics.
"""
if issuper(L):
# check initial state
if isket(rho0):
# Got a wave function as initial state: convert to density matrix.
rho0 = rho0 * rho0.dag()
w, v = L.eigenstates()
v = np.hstack([ket.full() for ket in v])
# w[i] = eigenvalue i
# v[:,i] = eigenvector i
rlen = np.prod(rho0.shape)
r0 = mat2vec(rho0.full())
v0 = la.solve(v, r0)
vv = v * sp.spdiags(v0.T, 0, rlen, rlen)
out = None
for i in range(rlen):
qo = Qobj(vec2mat(vv[:, i]), dims=rho0.dims, shape=rho0.shape)
if out:
out += eseries(qo, w[i])
else:
out = eseries(qo, w[i])
elif isoper(L):
if not isket(rho0):
raise TypeError("Second argument must be a ket if first" + "is a Hamiltonian.")
w, v = L.eigenstates()
v = np.hstack([ket.full() for ket in v])
# w[i] = eigenvalue i
# v[:,i] = eigenvector i
rlen = np.prod(rho0.shape)
r0 = rho0.full()
v0 = la.solve(v, r0)
vv = v * sp.spdiags(v0.T, 0, rlen, rlen)
out = None
for i in range(rlen):
qo = Qobj(np.matrix(vv[:, i]).T, dims=rho0.dims, shape=rho0.shape)
if out:
out += eseries(qo, -1.0j * w[i])
else:
out = eseries(qo, -1.0j * w[i])
else:
raise TypeError("First argument must be a Hamiltonian or Liouvillian.")
return estidy(out)
示例7: steadystate
def steadystate(A, c_op_list=[], method='direct', sparse=True, use_rcm=True,
sym=False, use_precond=True, M=None, drop_tol=1e-3,
fill_factor=12, diag_pivot_thresh=None, maxiter=1000, tol=1e-5,
verbose=False):
"""Calculates the steady state for quantum evolution subject to the
supplied Hamiltonian or Liouvillian operator and (if given a Hamiltonian) a
list of collapse operators.
If the user passes a Hamiltonian then it, along with the list of collapse
operators, will be converted into a Liouvillian operator in Lindblad form.
Parameters
----------
A : qobj
A Hamiltonian or Liouvillian operator.
c_op_list : list
A list of collapse operators.
method : str {'direct', 'iterative', 'iterative-bicg', 'lu', 'svd', 'power'}
Method for solving the underlying linear equation. Direct solver
'direct' (default), iterative GMRES method 'iterative',
iterative method BICGSTAB 'iterative-bicg', LU decomposition 'lu',
SVD 'svd' (dense), or inverse-power method 'power'.
sparse : bool, default=True
Solve for the steady state using sparse algorithms. If set to False,
the underlying Liouvillian operator will be converted into a dense
matrix. Use only for 'smaller' systems.
maxiter : int, optional
Maximum number of iterations to perform if using an iterative method
such as 'iterative' (default=1000), or 'power' (default=10).
tol : float, optional, default=1e-5
Tolerance used for terminating solver solution when using iterative
solvers.
use_precond : bool optional, default = True
ITERATIVE ONLY. Use an incomplete sparse LU decomposition as a
preconditioner for the 'iterative' LGMRES and BICG solvers.
Speeds up convergence time by orders of magnitude in many cases.
M : {sparse matrix, dense matrix, LinearOperator}
Preconditioner for A. The preconditioner should approximate the inverse of A.
Effective preconditioning dramatically improves the rate of convergence,
for iterative methods. Does not affect other solvers.
fill_factor : float, default=12
ITERATIVE ONLY. Specifies the fill ratio upper bound (>=1) of the iLU
preconditioner. Lower values save memory at the cost of longer
execution times and a possible singular factorization.
drop_tol : float, default=1e-3
ITERATIVE ONLY. Sets the threshold for the magnitude of preconditioner
elements that should be dropped. Can be reduced for a courser factorization
at the cost of an increased number of iterations, and a possible singular
factorization.
diag_pivot_thresh : float, default=None
ITERATIVE ONLY. Sets the threshold between [0,1] for which diagonal
elements are considered acceptable pivot points when using a
preconditioner. A value of zero forces the pivot to be the diagonal
element.
verbose : bool default=False
Flag for printing out detailed information on the steady state solver.
Returns
-------
dm : qobj
Steady state density matrix.
Notes
-----
The SVD method works only for dense operators (i.e. small systems).
"""
n_op = len(c_op_list)
if isoper(A):
if n_op == 0:
raise TypeError('Cannot calculate the steady state for a ' +
'non-dissipative system ' +
'(no collapse operators given)')
else:
A = liouvillian_fast(A, c_op_list)
if not issuper(A):
raise TypeError('Solving for steady states requires ' +
'Liouvillian (super) operators')
if method == 'direct':
if sparse:
return _steadystate_direct_sparse(A, verbose=verbose)
else:
return _steadystate_direct_dense(A, verbose=verbose)
#.........这里部分代码省略.........
示例8: mesolve
#.........这里部分代码省略.........
if isinstance(c_ops, Qobj):
c_ops = [c_ops]
if isinstance(e_ops, Qobj):
e_ops = [e_ops]
if isinstance(e_ops, dict):
e_ops_dict = e_ops
e_ops = [e for e in e_ops.values()]
else:
e_ops_dict = None
# convert array based time-dependence to string format
H, c_ops, args = _td_wrap_array_str(H, c_ops, args, tlist)
# check for type (if any) of time-dependent inputs
_, n_func, n_str = _td_format_check(H, c_ops)
if options is None:
options = Options()
if (not options.rhs_reuse) or (not config.tdfunc):
# reset config collapse and time-dependence flags to default values
config.reset()
res = None
#
# dispatch the appropriate solver
#
if (
(c_ops and len(c_ops) > 0)
or (not isket(rho0))
or (isinstance(H, Qobj) and issuper(H))
or (isinstance(H, list) and isinstance(H[0], Qobj) and issuper(H[0]))
):
#
# we have collapse operators
#
#
# find out if we are dealing with all-constant hamiltonian and
# collapse operators or if we have at least one time-dependent
# operator. Then delegate to appropriate solver...
#
if isinstance(H, Qobj):
# constant hamiltonian
if n_func == 0 and n_str == 0:
# constant collapse operators
res = _mesolve_const(H, rho0, tlist, c_ops, e_ops, args, options, progress_bar)
elif n_str > 0:
# constant hamiltonian but time-dependent collapse
# operators in list string format
res = _mesolve_list_str_td([H], rho0, tlist, c_ops, e_ops, args, options, progress_bar)
elif n_func > 0:
# constant hamiltonian but time-dependent collapse
# operators in list function format
res = _mesolve_list_func_td([H], rho0, tlist, c_ops, e_ops, args, options, progress_bar)
elif isinstance(H, (types.FunctionType, types.BuiltinFunctionType, partial)):
# function-callback style time-dependence: must have constant
# collapse operators
if n_str > 0: # or n_func > 0:
raise TypeError(
示例9: steadystate
#.........这里部分代码省略.........
explicitly specified.
use_precond : bool optional, default = True
ITERATIVE ONLY. Use an incomplete sparse LU decomposition as a
preconditioner for the 'iterative' GMRES and BICG solvers.
Speeds up convergence time by orders of magnitude in many cases.
M : {sparse matrix, dense matrix, LinearOperator}, optional
Preconditioner for A. The preconditioner should approximate the inverse
of A. Effective preconditioning dramatically improves the rate of
convergence, for iterative methods only . If no preconditioner is
given and ``use_precond=True``, then one is generated automatically.
fill_factor : float, optional, default=10
ITERATIVE ONLY. Specifies the fill ratio upper bound (>=1) of the iLU
preconditioner. Lower values save memory at the cost of longer
execution times and a possible singular factorization.
drop_tol : float, optional, default=1e-3
ITERATIVE ONLY. Sets the threshold for the magnitude of preconditioner
elements that should be dropped. Can be reduced for a courser
factorization at the cost of an increased number of iterations, and a
possible singular factorization.
diag_pivot_thresh : float, optional, default=None
ITERATIVE ONLY. Sets the threshold between [0,1] for which diagonal
elements are considered acceptable pivot points when using a
preconditioner. A value of zero forces the pivot to be the diagonal
element.
ILU_MILU : str, optional, default='smilu_2'
Selects the incomplete LU decomposition method algoithm used in
creating the preconditoner. Should only be used by advanced users.
Returns
-------
dm : qobj
Steady state density matrix.
Notes
-----
The SVD method works only for dense operators (i.e. small systems).
"""
ss_args = _default_steadystate_args()
for key in kwargs.keys():
if key in ss_args.keys():
ss_args[key] = kwargs[key]
else:
raise Exception(
"Invalid keyword argument '"+key+"' passed to steadystate.")
# Set column perm to NATURAL if using RCM and not specified by user
if ss_args['use_rcm'] and ('permc_spec' not in kwargs.keys()):
ss_args['permc_spec'] = 'NATURAL'
# Set use_wbm=True if using iterative solver with preconditioner and
# not explicitly set to False by user
if (ss_args['method'] in ['iterative-lgmres', 'iterative-gmres']) \
and ('use_wbm' not in kwargs.keys()):
ss_args['use_wbm'] = True
n_op = len(c_op_list)
if isoper(A):
if n_op == 0:
raise TypeError('Cannot calculate the steady state for a ' +
'non-dissipative system ' +
'(no collapse operators given)')
else:
A = liouvillian(A, c_op_list)
if not issuper(A):
raise TypeError('Solving for steady states requires ' +
'Liouvillian (super) operators')
# Set weight parameter to avg abs val in L if not set explicitly
if 'weight' not in kwargs.keys():
ss_args['weight'] = np.mean(np.abs(A.data.data.max()))
if ss_args['method'] == 'direct':
if ss_args['sparse']:
return _steadystate_direct_sparse(A, ss_args)
else:
return _steadystate_direct_dense(A)
elif ss_args['method'] == 'eigen':
return _steadystate_eigen(A, ss_args)
elif ss_args['method'] in ['iterative-gmres', 'iterative-lgmres']:
return _steadystate_iterative(A, ss_args)
elif ss_args['method'] == 'svd':
return _steadystate_svd_dense(A, ss_args)
elif ss_args['method'] == 'power':
return _steadystate_power(A, ss_args)
else:
raise ValueError('Invalid method argument for steadystate.')
示例10: _mesolve_list_func_td
def _mesolve_list_func_td(H_list, rho0, tlist, c_list, e_ops, args, opt, progress_bar):
"""
Internal function for solving the master equation. See mesolve for usage.
"""
if debug:
print(inspect.stack()[0][3])
#
# check initial state
#
if isket(rho0):
rho0 = rho0 * rho0.dag()
#
# construct liouvillian in list-function format
#
L_list = []
if opt.rhs_with_state:
constant_func = lambda x, y, z: 1.0
else:
constant_func = lambda x, y: 1.0
# add all hamitonian terms to the lagrangian list
for h_spec in H_list:
if isinstance(h_spec, Qobj):
h = h_spec
h_coeff = constant_func
elif isinstance(h_spec, list) and isinstance(h_spec[0], Qobj):
h = h_spec[0]
h_coeff = h_spec[1]
else:
raise TypeError("Incorrect specification of time-dependent " + "Hamiltonian (expected callback function)")
if isoper(h):
L_list.append([(-1j * (spre(h) - spost(h))).data, h_coeff, False])
elif issuper(h):
L_list.append([h.data, h_coeff, False])
else:
raise TypeError(
"Incorrect specification of time-dependent " + "Hamiltonian (expected operator or superoperator)"
)
# add all collapse operators to the liouvillian list
for c_spec in c_list:
if isinstance(c_spec, Qobj):
c = c_spec
c_coeff = constant_func
c_square = False
elif isinstance(c_spec, list) and isinstance(c_spec[0], Qobj):
c = c_spec[0]
c_coeff = c_spec[1]
c_square = True
else:
raise TypeError(
"Incorrect specification of time-dependent " + "collapse operators (expected callback function)"
)
if isoper(c):
L_list.append([liouvillian(None, [c], data_only=True), c_coeff, c_square])
elif issuper(c):
L_list.append([c.data, c_coeff, c_square])
else:
raise TypeError(
"Incorrect specification of time-dependent "
+ "collapse operators (expected operator or "
+ "superoperator)"
)
#
# setup integrator
#
initial_vector = mat2vec(rho0.full()).ravel()
if opt.rhs_with_state:
r = scipy.integrate.ode(drho_list_td_with_state)
else:
r = scipy.integrate.ode(drho_list_td)
r.set_integrator(
"zvode",
method=opt.method,
order=opt.order,
atol=opt.atol,
rtol=opt.rtol,
nsteps=opt.nsteps,
first_step=opt.first_step,
min_step=opt.min_step,
max_step=opt.max_step,
)
r.set_initial_value(initial_vector, tlist[0])
r.set_f_params(L_list, args)
#.........这里部分代码省略.........
示例11: _mesolve_list_str_td
def _mesolve_list_str_td(H_list, rho0, tlist, c_list, e_ops, args, opt, progress_bar):
"""
Internal function for solving the master equation. See mesolve for usage.
"""
if debug:
print(inspect.stack()[0][3])
#
# check initial state: must be a density matrix
#
if isket(rho0):
rho0 = rho0 * rho0.dag()
#
# construct liouvillian
#
Lconst = 0
Ldata = []
Linds = []
Lptrs = []
Lcoeff = []
# loop over all hamiltonian terms, convert to superoperator form and
# add the data of sparse matrix representation to
for h_spec in H_list:
if isinstance(h_spec, Qobj):
h = h_spec
if isoper(h):
Lconst += -1j * (spre(h) - spost(h))
elif issuper(h):
Lconst += h
else:
raise TypeError(
"Incorrect specification of time-dependent "
+ "Hamiltonian (expected operator or "
+ "superoperator)"
)
elif isinstance(h_spec, list):
h = h_spec[0]
h_coeff = h_spec[1]
if isoper(h):
L = -1j * (spre(h) - spost(h))
elif issuper(h):
L = h
else:
raise TypeError(
"Incorrect specification of time-dependent "
+ "Hamiltonian (expected operator or "
+ "superoperator)"
)
Ldata.append(L.data.data)
Linds.append(L.data.indices)
Lptrs.append(L.data.indptr)
Lcoeff.append(h_coeff)
else:
raise TypeError("Incorrect specification of time-dependent " + "Hamiltonian (expected string format)")
# loop over all collapse operators
for c_spec in c_list:
if isinstance(c_spec, Qobj):
c = c_spec
if isoper(c):
cdc = c.dag() * c
Lconst += spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)
elif issuper(c):
Lconst += c
else:
raise TypeError(
"Incorrect specification of time-dependent "
+ "Liouvillian (expected operator or "
+ "superoperator)"
)
elif isinstance(c_spec, list):
c = c_spec[0]
c_coeff = c_spec[1]
if isoper(c):
cdc = c.dag() * c
L = spre(c) * spost(c.dag()) - 0.5 * spre(cdc) - 0.5 * spost(cdc)
c_coeff = "(" + c_coeff + ")**2"
elif issuper(c):
L = c
else:
raise TypeError(
"Incorrect specification of time-dependent "
+ "Liouvillian (expected operator or "
+ "superoperator)"
)
#.........这里部分代码省略.........
示例12: _steadystate_power
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: _mesolve_func_td
def _mesolve_func_td(L_func, rho0, tlist, c_op_list, e_ops, args, opt,
progress_bar):
"""
Evolve the density matrix using an ODE solver with time dependent
Hamiltonian.
"""
if debug:
print(inspect.stack()[0][3])
#
# check initial state
#
if isket(rho0):
rho0 = ket2dm(rho0)
#
# construct liouvillian
#
new_args = None
if len(c_op_list) > 0:
L_data = liouvillian(None, c_op_list).data
else:
n, m = rho0.shape
if issuper(rho0):
L_data = sp.csr_matrix((n, m), dtype=complex)
else:
L_data = sp.csr_matrix((n ** 2, m ** 2), dtype=complex)
if type(args) is dict:
new_args = {}
for key in args:
if isinstance(args[key], Qobj):
if isoper(args[key]):
new_args[key] = (
-1j * (spre(args[key]) - spost(args[key])))
else:
new_args[key] = args[key]
else:
new_args[key] = args[key]
elif type(args) is list or type(args) is tuple:
new_args = []
for arg in args:
if isinstance(arg, Qobj):
if isoper(arg):
new_args.append((-1j * (spre(arg) - spost(arg))).data)
else:
new_args.append(arg.data)
else:
new_args.append(arg)
if type(args) is tuple:
new_args = tuple(new_args)
else:
if isinstance(args, Qobj):
if isoper(args):
new_args = (-1j * (spre(args) - spost(args)))
else:
new_args = args
else:
new_args = args
#
# setup integrator
#
initial_vector = mat2vec(rho0.full()).ravel('F')
if issuper(rho0):
if not opt.rhs_with_state:
r = scipy.integrate.ode(_ode_super_func_td)
else:
r = scipy.integrate.ode(_ode_super_func_td_with_state)
else:
if not opt.rhs_with_state:
r = scipy.integrate.ode(cy_ode_rho_func_td)
else:
r = scipy.integrate.ode(_ode_rho_func_td_with_state)
r.set_integrator('zvode', method=opt.method, order=opt.order,
atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
first_step=opt.first_step, min_step=opt.min_step,
max_step=opt.max_step)
r.set_initial_value(initial_vector, tlist[0])
r.set_f_params(L_data, L_func, new_args)
#
# call generic ODE code
#
return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)
示例14: _steadystate_power
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
示例15: steadystate
def steadystate(A, c_op_list=[], **kwargs):
"""Calculates the steady state for quantum evolution subject to the
supplied Hamiltonian or Liouvillian operator and (if given a Hamiltonian) a
list of collapse operators.
If the user passes a Hamiltonian then it, along with the list of collapse
operators, will be converted into a Liouvillian operator in Lindblad form.
Parameters
----------
A : qobj
A Hamiltonian or Liouvillian operator.
c_op_list : list
A list of collapse operators.
method : str {'direct', 'iterative', 'iterative-bicg', 'svd', 'power'}
Method for solving the underlying linear equation. Direct solver
'direct' (default), iterative GMRES method 'iterative',
iterative method BICGSTAB 'iterative-bicg', SVD 'svd' (dense),
or inverse-power method 'power'.
sparse : bool, {True, False}
Solve for the steady state using sparse algorithms. If set to False,
the underlying Liouvillian operator will be converted into a dense
matrix. Use only for 'smaller' systems.
use_rcm : bool, {True, False}
Use reverse Cuthill-Mckee reordering to minimize fill-in in the
LU factorization of the Liouvillian.
use_umfpack : bool {False, True}
Use umfpack solver instead of SuperLU. For SciPy 0.14+, this option
requires installing scikits.umfpack.
maxiter : int, optional
Maximum number of iterations to perform if using an iterative method
such as 'iterative' (default=1000), or 'power' (default=10).
tol : float, optional, default=1e-5
Tolerance used for terminating solver solution when using iterative
solvers.
use_precond : bool optional, default = True
ITERATIVE ONLY. Use an incomplete sparse LU decomposition as a
preconditioner for the 'iterative' GMRES and BICG solvers.
Speeds up convergence time by orders of magnitude in many cases.
M : {sparse matrix, dense matrix, LinearOperator}, optional
Preconditioner for A. The preconditioner should approximate the inverse
of A. Effective preconditioning dramatically improves the rate of
convergence, for iterative methods. Does not affect other solvers.
fill_factor : float, default=12
ITERATIVE ONLY. Specifies the fill ratio upper bound (>=1) of the iLU
preconditioner. Lower values save memory at the cost of longer
execution times and a possible singular factorization.
drop_tol : float, default=1e-3
ITERATIVE ONLY. Sets the threshold for the magnitude of preconditioner
elements that should be dropped. Can be reduced for a courser
factorization at the cost of an increased number of iterations, and a
possible singular factorization.
diag_pivot_thresh : float, default=None
ITERATIVE ONLY. Sets the threshold between [0,1] for which diagonal
elements are considered acceptable pivot points when using a
preconditioner. A value of zero forces the pivot to be the diagonal
element.
Returns
-------
dm : qobj
Steady state density matrix.
Notes
-----
The SVD method works only for dense operators (i.e. small systems).
"""
ss_args = _default_steadystate_args()
for key in kwargs.keys():
if key in ss_args.keys():
ss_args[key]=kwargs[key]
else:
raise Exception('Invalid keyword argument passed to steadystate.')
n_op = len(c_op_list)
if isoper(A):
if n_op == 0:
raise TypeError('Cannot calculate the steady state for a ' +
'non-dissipative system ' +
'(no collapse operators given)')
else:
A = liouvillian(A, c_op_list)
if not issuper(A):
#.........这里部分代码省略.........