How to use the qutip.superoperator.spre function in qutip

tr_op = tensor([identity(n) for n in L.dims[0][0]])
            tr_op_vec = operator_to_vector(tr_op)
            
            Pop = sp.kron(rhoss_vec.data, tr_op_vec.data.T, format='csr')
            Iop = sp.eye(N*N, N*N, format='csr')
            Q = Iop - Pop
            
            for k,w in enumerate(wlist):
                
                if w != 0.0:    
                    L_temp = 1.0j*w*spre(tr_op) + L
                else: #At zero frequency some solvers fail for small systems.
                      #Adding a small finite frequency of order 1e-15
                      #helps prevent the solvers from throwing an exception.
                    L_temp =  1.0j*(1e-15)*spre(tr_op) + L
                    
                if not settings.has_mkl:
                    A = L_temp.data.tocsc()
                else:
                    A = L_temp.data.tocsr()
                    A.sort_indices()                      
                      
                rhoss_vec = mat2vec(rhoss.full()).ravel()               
                
                for j, Jj in enumerate(J_ops):
                    Qj = Q.dot( Jj.data.dot( rhoss_vec))
                    try:
                        if settings.has_mkl:
                            X_rho_vec_j = mkl_spsolve(A,Qj)                            
                        else:
                            X_rho_vec_j = sp.linalg.splu(A, permc_spec

sops += [c / np.sqrt(2), -1j / np.sqrt(2) * c]
        sso.m_ops = m_ops
        if not isinstance(sso.dW_factors, list):
            sso.dW_factors = [np.sqrt(2)] * len(sops)
        elif len(sso.dW_factors) == len(sso.sc_ops):
            dW_factors = []
            for fact in sso.dW_factors:
                dW_factors += [np.sqrt(2) * fact, np.sqrt(2) * fact]
            sso.dW_factors = dW_factors
        elif len(sso.dW_factors) != len(sops):
            raise Exception("The len of dW_factors is not the same as sc_ops")
    else:
        raise Exception("The method must be one of homodyne or heterodyne")

    LH = 1 - (sso.H * 1j * sso.dt)
    sso.pp = spre(sso.H) * 0
    sso.sops = []
    sso.preops = []
    sso.postops = []
    sso.preops2 = []
    sso.postops2 = []

    def _prespostdag(op):
        return spre(op) * spost(op.dag())

    for op in sso.c_ops:
        LH -= op._cdc() * sso.dt * 0.5
        sso.pp += op.apply(_prespostdag)._f_norm2() * sso.dt

    for i, op in enumerate(sops):
        LH -= op._cdc() * sso.dt * 0.5
        sso.sops += [(spre(op) + spost(op.dag())) * sso.dt]

for op in sso.c_ops:
        LH -= op._cdc() * sso.dt * 0.5
        sso.pp += op.apply(_prespostdag)._f_norm2() * sso.dt

    for i, op in enumerate(sops):
        LH -= op._cdc() * sso.dt * 0.5
        sso.sops += [(spre(op) + spost(op.dag())) * sso.dt]
        sso.preops += [spre(op)]
        sso.postops += [spost(op.dag())]
        for op2 in sops[i:]:
            sso.preops2 += [spre(op * op2)]
            sso.postops2 += [spost(op.dag() * op2.dag())]

    sso.ce_ops = [QobjEvo(spre(op)) for op in sso.e_ops]
    sso.cm_ops = [QobjEvo(spre(op)) for op in sso.m_ops]
    sso.preLH = spre(LH)
    sso.postLH = spost(LH.dag())

    sso.preLH.compile()
    sso.postLH.compile()
    sso.pp.compile()
    [op.compile() for op in sso.sops]
    [op.compile() for op in sso.preops]
    [op.compile() for op in sso.postops]
    [op.compile() for op in sso.preops2]
    [op.compile() for op in sso.postops2]
    [op.compile() for op in sso.cm_ops]
    [op.compile() for op in sso.ce_ops]

    sso.solver_obj = PmSMESolver
    sso.solver_name = "smesolve_" + sso.solver
    res = _sesolve_generic(sso, sso.options, sso.progress_bar)

Internal function for computing the pseudo inverse of an Liouvillian using
    dense matrix methods. See pseudo_inverse for details.
    """
    rho_vec = np.transpose(mat2vec(rhoss.full()))

    tr_mat = tensor([identity(n) for n in L.dims[0][0]])
    tr_vec = np.transpose(mat2vec(tr_mat.full()))
    N = np.prod(L.dims[0][0])
    I = np.identity(N * N)
    P = np.kron(np.transpose(rho_vec), tr_vec)
    Q = I - P

    if w is None:
        L = L
    else:
        L = 1.0j*w*spre(tr_mat)+L

    if pseudo_args['method'] == 'direct':
        try:
            LIQ = np.linalg.solve(L.full(), Q)
        except:
            LIQ = np.linalg.lstsq(L.full(), Q)[0]

        R = np.dot(Q, LIQ)

        return Qobj(R, dims=L.dims)

    elif pseudo_args['method'] == 'numpy':
        return Qobj(np.dot(Q, np.dot(np.linalg.pinv(L.full()), Q)),
                    dims=L.dims)

    elif pseudo_args['method'] == 'scipy':

N_c = self.N_cut
        N_m = self.N_exp
        Q = coup_op # Q as shorthand for coupling operator
        beta = 1.0/(self.boltzmann*self.temperature)

        # Ntot is the total number of ancillary elements in the hierarchy
        # Ntot = factorial(N_c + N_m) / (factorial(N_c)*factorial(N_m))
        # Turns out to be the same as nstates from state_number_enumerate
        N_he, he2idx, idx2he = enr_state_dictionaries([N_c + 1]*N_m , N_c)

        unit_helems = fast_identity(N_he)
        if self.bnd_cut_approx:
            # the Tanimura boundary cut off operator
            if stats:
                stats.add_message('options', 'boundary cutoff approx', ss_conf)
            op = -2*spre(Q)*spost(Q.dag()) + spre(Q.dag()*Q) + spost(Q.dag()*Q)

            approx_factr = ((2*lam0 / (beta*gam*hbar)) - 1j*lam0) / hbar
            for k in range(N_m):
                approx_factr -= (c[k] / nu[k])
            L_bnd = -approx_factr*op.data
            L_helems = zcsr_kron(unit_helems, L_bnd)
        else:
            L_helems = fast_csr_matrix(shape=(N_he*sup_dim, N_he*sup_dim))

        # Build the hierarchy element interaction matrix
        if stats: start_helem_constr = timeit.default_timer()

        unit_sup = spre(unit_sys).data
        spreQ = spre(Q).data
        spostQ = spost(Q).data
        commQ = (spre(Q) - spost(Q)).data

Returns
    -------
    chi : array
        QPT chi matrix

    """

    E_ops = []
    # loop over all index permutations
    for inds in _index_permutations([len(ops) for ops in op_basis_list]):
        # loop over all composite systems
        E_op_list = [op_basis_list[k][inds[k]] for k in range(len(
            op_basis_list))]
        E_ops.append(tensor(E_op_list))

    EE_ops = [spre(E1) * spost(E2.dag()) for E1 in E_ops for E2 in E_ops]

    M = hstack([mat2vec(EE.full()) for EE in EE_ops])

    Uvec = mat2vec(U.full())

    chi_vec = la.solve(M, Uvec)

    return vec2mat(chi_vec)

def _prespostdag(op):
        return spre(op) * spost(op.dag())

def H2L(self, t, rho, args):
        Ht = self.f(t, args)
        Lt = -1.0j * (spre(Ht) - spost(Ht)).data
        for op in self.c_ops:
            Lt += op(t).data
        return Lt

rhoss_vec = operator_to_vector(rhoss)

    tr_op = tensor([identity(n) for n in L.dims[0][0]])
    tr_op_vec = operator_to_vector(tr_op)

    P = zcsr_kron(rhoss_vec.data, tr_op_vec.data.T)
    I = sp.eye(N*N, N*N, format='csr')
    Q = I - P

    if w is None:
        L = 1.0j*(1e-15)*spre(tr_op) + L
    else:
        if w != 0.0:
            L = 1.0j*w*spre(tr_op) + L
        else:
            L = 1.0j*(1e-15)*spre(tr_op) + L

    if pseudo_args['use_rcm']:
        perm = sp.csgraph.reverse_cuthill_mckee(L.data)
        A = sp_permute(L.data, perm, perm)
        Q = sp_permute(Q, perm, perm)
    else:
        if ss_args['solver'] == 'scipy':
            A = L.data.tocsc()
            A.sort_indices()

    if pseudo_args['method'] == 'splu':
        if settings.has_mkl:
            A = L.data.tocsr()
            A.sort_indices()
            LIQ = mkl_spsolve(A, Q.toarray())
        else:

sso.sops = []
    sso.preops = []
    sso.postops = []
    sso.preops2 = []
    sso.postops2 = []

    def _prespostdag(op):
        return spre(op) * spost(op.dag())

    for op in sso.c_ops:
        LH -= op._cdc() * sso.dt * 0.5
        sso.pp += op.apply(_prespostdag)._f_norm2() * sso.dt

    for i, op in enumerate(sops):
        LH -= op._cdc() * sso.dt * 0.5
        sso.sops += [(spre(op) + spost(op.dag())) * sso.dt]
        sso.preops += [spre(op)]
        sso.postops += [spost(op.dag())]
        for op2 in sops[i:]:
            sso.preops2 += [spre(op * op2)]
            sso.postops2 += [spost(op.dag() * op2.dag())]

    sso.ce_ops = [QobjEvo(spre(op)) for op in sso.e_ops]
    sso.cm_ops = [QobjEvo(spre(op)) for op in sso.m_ops]
    sso.preLH = spre(LH)
    sso.postLH = spost(LH.dag())

    sso.preLH.compile()
    sso.postLH.compile()
    sso.pp.compile()
    [op.compile() for op in sso.sops]
    [op.compile() for op in sso.preops]