quspin.tools.misc.matvec

quspin.tools.misc.matvec(array, other, overwrite_out=False, out=None, a=1.0)[source]

Calculates omp-parallelized matrix vector products.

Let \(A\) be a matrix (array), \(v\) be a vector (other), and \(a\) (a) be a scalar. This function provides an omp-parallelized implementation of

\[x += aAv \qquad \mathrm{or} \qquad x = aAv.\]
Parameters:
arrayarray_like object (e.g. numpy.ndarray, scipy.sparse.csr_matrix, scipy.sparse.csc_matrix, scipy.sparse.dia_matrix)

Sparse or dense array to take the dot product with.

otherarray_like

array which contains the vector to take the dot product with.

ascalar, optional

value to scale the vector with after the product with array is taken: \(x += a A v\) or \(x = a A v\).

outarray_like

output array to put the results of the calculation in.

overwrite_outbool, optional

If set to True, the function overwrites the values in out with the result (cf. \(x = a A v\)). Otherwise the result is added to the values in out (in-pace addition, cf. \(x += a A v\)).

Returns:
numpy.ndarray

result of the matrix-vector product \(a A v\).

  • if out is not None and overwrite_out = True, the function returns out with the original data overwritten, otherwise if overwrite_out = False the result is added to out.

  • if out is None, the result is stored in a new array which is returned by the function.

Notes

  • for QuSpin builds which support OpenMP, this function will be multithreaded.

  • using out=v will result in incorrect results.

  • matvec determines the correct omp-parallelized matrix-vector product, depending on the type of the input array (csr, csc, dia, other [e.g., dense]), every time matvec is called. To avoid this overhead, see quspin.tools.misc.get_matvec_function().

Examples

The example shows how to use the get_matvec_function() (line 43) and the matvec() function (lines 47-81) in a user-defined ODE which solves the Lindblad equation for a single qubit (see also Example 17).

  1#
  2from quspin.operators import hamiltonian, commutator, anti_commutator
  3from quspin.basis import spin_basis_1d  # Hilbert space spin basis_1d
  4from quspin.tools.evolution import evolve
  5from quspin.tools.misc import get_matvec_function
  6import numpy as np
  7from six import iteritems  # loop over elements of dictionary
  8import matplotlib.pyplot as plt  # plotting library
  9
 10#
 11###### model parameters
 12#
 13L = 1  # one qubit
 14delta = 1.0  # detuning
 15Omega_0 = np.sqrt(2)  # bare Rabi frequency
 16Omega_Rabi = np.sqrt(Omega_0**2 + delta**2)  # Rabi frequency
 17gamma = 0.5 * np.sqrt(3)  # decay rate
 18#
 19##### create Hamiltonian to evolve unitarily
 20# basis
 21basis = spin_basis_1d(
 22    L, pauli=-1
 23)  # uses convention "+" = [[0,1],[0,0]] (mind the spin factor 1/2)
 24# site-coupling lists
 25hx_list = [[Omega_0, i] for i in range(L)]
 26hz_list = [[delta, i] for i in range(L)]
 27# static opstr list
 28static_H = [["x", hx_list], ["z", hz_list]]
 29# hamiltonian
 30H = hamiltonian(static_H, [], basis=basis, dtype=np.float64)
 31print("Hamiltonian:\n", H.toarray())
 32#
 33##### create Lindbladian
 34# site-coupling lists
 35L_list = [[1.0, i] for i in range(L)]
 36# static opstr list
 37static_L = [["+", L_list]]
 38# Lindblad operator
 39L = hamiltonian(static_L, [], basis=basis, dtype=np.complex128, check_herm=False)
 40print("Lindblad operator:\n", L.toarray())
 41# pre-compute operators for efficiency
 42L_dagger = L.getH()
 43L_daggerL = L_dagger * L
 44#
 45#### determine the corresponding matvec routines ####
 46#
 47# different matvec functions are required since we use both matrices and their transposed, cf. Lindblad_EOM_v3
 48matvec_csr = get_matvec_function(H.static)  # csr matvec function
 49matvec_csc = get_matvec_function(H.static.T)  # csc matvec function
 50
 51
 52#
 53# fast function (not as memory efficient)
 54def Lindblad_EOM(time, rho, rho_out, rho_aux):
 55    """
 56    This function solves the complex-valued time-dependent GPE:
 57    $$ \dot\rho(t) = -i[H,\rho(t)] + 2\gamma\left( L\rho L^\dagger - \frac{1}{2}\{L^\dagger L, \rho \} \right) $$
 58    """
 59    rho = rho.reshape((H.Ns, H.Ns))  # reshape vector from ODE solver input
 60    ### Hamiltonian part
 61    # commutator term (unitary
 62    # rho_out = H.static.dot(rho))
 63    matvec_csr(H.static, rho, out=rho_out, a=+1.0, overwrite_out=True)
 64    # rho_out -= (H.static.T.dot(rho.T)).T // RHS~rho.dot(H)
 65    matvec_csc(H.static.T, rho.T, out=rho_out.T, a=-1.0, overwrite_out=False)
 66    #
 67    for func, Hd in iteritems(H._dynamic):
 68        ft = func(time)
 69        # rho_out += ft*Hd.dot(rho)
 70        matvec_csr(Hd, rho, out=rho_out, a=+ft, overwrite_out=False)
 71        # rho_out -= ft*(Hd.T.dot(rho.T)).T
 72        matvec_csc(Hd.T, rho.T, out=rho_out.T, a=-ft, overwrite_out=False)
 73    # multiply by -i
 74    rho_out *= -1.0j
 75    #
 76    ### Lindbladian part (static only)
 77    ## 1st Lindblad term (nonunitary): 2\gamma * L*rho*L^\dagger
 78    # rho_aux = 2\gamma*L.dot(rho)
 79    matvec_csr(L.static, rho, out=rho_aux, a=+2.0 * gamma, overwrite_out=True)
 80    # rho_out += (L.static.conj().dot(rho_aux.T)).T // RHS ~ rho_aux.dot(L_dagger)
 81    matvec_csr(L.static.conj(), rho_aux.T, out=rho_out.T, a=+1.0, overwrite_out=False)
 82    #
 83    ## anticommutator (2nd Lindblad) term (nonunitary): -\gamma \{L^\dagger * L, \eho}
 84    # rho_out += gamma*L_daggerL.static.dot(rho)
 85    matvec_csr(L_daggerL.static, rho, out=rho_out, a=-gamma, overwrite_out=False)
 86    # rho_out += gamma*(L_daggerL.static.T.dot(rho.T)).T // RHS~rho.dot(L_daggerL)
 87    matvec_csc(L_daggerL.static.T, rho.T, out=rho_out.T, a=-gamma, overwrite_out=False)
 88    #
 89    return rho_out.ravel()  # ODE solver accepts vectors only
 90
 91
 92#
 93# define auxiliary arguments
 94EOM_args = (
 95    np.zeros(
 96        (H.Ns, H.Ns), dtype=np.complex128, order="C"
 97    ),  # auxiliary variable rho_out
 98    np.zeros((H.Ns, H.Ns), dtype=np.complex128, order="C"),
 99)  # auxiliary variable rho_aux
100#
101##### time-evolve state according to Lindlad equation
102# define real time vector
103t_max = 6.0
104time = np.linspace(0.0, t_max, 101)
105# define initial state
106rho0 = np.array([[0.5, 0.5j], [-0.5j, 0.5]], dtype=np.complex128)
107#
108# evolution
109rho_t = evolve(
110    rho0,
111    time[0],
112    time,
113    Lindblad_EOM,
114    f_params=EOM_args,
115    iterate=True,
116    atol=1e-12,
117    rtol=1e-12,
118)
119#
120# compute state evolution
121population_down = np.zeros(time.shape, dtype=np.float64)
122for i, rho_flattened in enumerate(rho_t):
123    rho = rho_flattened.reshape(H.Ns, H.Ns)
124    population_down[i] = rho_flattened[1, 1].real
125    print(
126        "time={0:.2f}, population of down state = {1:0.8f}".format(
127            time[i], population_down[i]
128        )
129    )