-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathrhs8134.pyx
55 lines (49 loc) · 2.27 KB
/
rhs8134.pyx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#!python
#cython: language_level=3
# This file is generated automatically by QuTiP.
# (C) 2011 and later, QuSTaR
import numpy as np
cimport numpy as np
cimport cython
np.import_array()
cdef extern from "numpy/arrayobject.h" nogil:
void PyDataMem_NEW_ZEROED(size_t size, size_t elsize)
void PyArray_ENABLEFLAGS(np.ndarray arr, int flags)
void PyDataMem_FREE(void * ptr)
from qutip.cy.interpolate cimport interp, zinterp
from qutip.cy.math cimport erf, zerf
cdef double pi = 3.14159265358979323
from qutip.cy.brtools cimport (dense_add_mult, ZHEEVR, dense_to_eigbasis,
vec_to_eigbasis, vec_to_fockbasis, skew_and_dwmin,
diag_liou_mult, spec_func, farray_alloc)
from qutip.cy.brtools cimport (cop_super_mult, br_term_mult)
include '/home/neill/anaconda3/lib/python3.7/site-packages/qutip-4.5.0.dev0-py3.7-linux-x86_64.egg/qutip/cy/complex_math.pxi'
cdef complex spectral0(double w, double t): return 2.0 * 0.05 / (pi * 10.0 * 2.5448096943108576e-14) if (w==0) else (2.0*0.05*10.0 *w /(pi*(w**2+10.0**2))) * ((1/(exp((w) * 2.5448096943108576e-14)-1))+1)
@cython.cdivision(True)
@cython.boundscheck(False)
@cython.wraparound(False)
def cy_td_ode_rhs(
double t,
complex[::1] vec,
complex[::1,:] H0,
complex[::1,:] A0,
unsigned int nrows):
cdef double complex * out = <complex *>PyDataMem_NEW_ZEROED(nrows**2,sizeof(complex))
cdef complex[::1, :] H = farray_alloc(nrows)
cdef complex[::1, :] evecs = farray_alloc(nrows)
cdef double * eigvals = <double *>PyDataMem_NEW_ZEROED(nrows,sizeof(double))
dense_add_mult(H, H0, 1)
ZHEEVR(H, eigvals, evecs, nrows)
PyDataMem_FREE(&H[0,0])
cdef double complex * eig_vec = vec_to_eigbasis(vec, evecs, nrows)
diag_liou_mult(eigvals, eig_vec, out, nrows)
cdef double[:,::1] skew = <double[:nrows,:nrows]><double *>PyDataMem_NEW_ZEROED(nrows**2,sizeof(double))
cdef double dw_min = skew_and_dwmin(eigvals, skew, nrows)
br_term_mult(t, A0, evecs, skew, dw_min, spectral0, eig_vec, out, nrows, 1, 0.1, 1e-12)
cdef np.ndarray[complex, ndim=1, mode='c'] arr_out = vec_to_fockbasis(out, evecs, nrows)
PyDataMem_FREE(&skew[0,0])
PyDataMem_FREE(&evecs[0,0])
PyDataMem_FREE(eigvals)
PyDataMem_FREE(eig_vec)
PyDataMem_FREE(out)
return arr_out