def test_02_2_qft_bounds(self):
        """
        control.pulseoptim: QFT gate with linear initial pulses (bounds)
        assert that amplitudes remain in bounds
        """
        Sx = sigmax()
        Sy = sigmay()
        Sz = sigmaz()
        Si = 0.5*identity(2)

        H_d = 0.5*(tensor(Sx, Sx) + tensor(Sy, Sy) + tensor(Sz, Sz))
        H_c = [tensor(Sx, Si), tensor(Sy, Si), tensor(Si, Sx), tensor(Si, Sy)]
        U_0 = identity(4)
        # Target for the gate evolution - Quantum Fourier Transform gate
        U_targ = qft.qft(2)

        n_ts = 10
        evo_time = 10

        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-9,
                        amp_lbound=-1.0, amp_ubound=1.0,
                        init_pulse_type='LIN',
                        gen_stats=True)
        assert_((result.final_amps >= -1.0).all() and
                    (result.final_amps <= 1.0).all(),
                    msg="Amplitude bounds exceeded for QFT")

def Hamiltonian_toymodel(theta, phi):

## Defining Constants
	omega = 1.4 * 10**6
	m_b = 5.7883818066*10**-5
	par=m_b*2.417990504024*10**11 # s-1/mT
	a1 = np.array([-.0989, -.0989, 1.7569]) * par # /mT
	a2 = np.array([0.0,0.0,1.0812]) * par

	H_a_zeeman_nt = omega * (Sx * mt.sin(theta) *  mt.cos(phi) + Sy *  mt.sin(theta) * mt. sin(phi) + Sz * mt.cos(theta))
	Ha_zeeman  = qt.tensor(H_a_zeeman_nt, qt.identity(3))

	H_a_hfi = a1[0] *  qt.tensor(Sx, I2_x) + a1[1] * qt.tensor(Sy, I2_y) + a1[2] * qt.tensor(Sz, I2_z)

	H_a = Ha_zeeman + H_a_hfi
	H_a_eigen, H_a_eigenstates = H_a.eigenstates()


	H_b_zeeman_nt = omega * (Sx * mt.sin(theta) *  mt.cos(phi) + Sy *  mt.sin(theta) * mt. sin(phi) + Sz * mt.cos(theta))
	Hb_zeeman  = qt.tensor(H_b_zeeman_nt, qt.identity(3))

	H_b_hfi = a2[0] *  qt.tensor(Sx, I2_x) + a2[1] * qt.tensor(Sy, I2_y) + a2[2] * qt.tensor(Sz, I2_z)

	H_b = Hb_zeeman + H_b_hfi
	H_b_eigen, H_b_eigenstates = H_b.eigenstates()

	return (H_a_eigen, H_a_eigenstates, H_b_eigen, H_b_eigenstates)

def testJCZeroTemperature():
    """
    brmesolve: Jaynes-Cummings model, zero temperature
    """

    N = 10
    a = tensor(destroy(N), identity(2))
    sm = tensor(identity(N), destroy(2))
    psi0 = ket2dm(tensor(basis(N, 1), basis(2, 0)))
    a_ops = [(a + a.dag())]
    e_ops = [a.dag() * a, sm.dag() * sm]

    w0 = 1.0 * 2 * np.pi
    g = 0.05 * 2 * np.pi
    kappa = 0.05
    times = np.linspace(0, 2 * 2 * np.pi / g, 1000)

    c_ops = [np.sqrt(kappa) * a]
    H = w0 * a.dag() * a + w0 * sm.dag() * sm + \
        g * (a + a.dag()) * (sm + sm.dag())

    res_me = mesolve(H, psi0, times, c_ops, e_ops)
    res_brme = brmesolve(H, psi0, times, a_ops, e_ops,
                         spectra_cb=[lambda w: kappa * (w >= 0)])

    for idx, e in enumerate(e_ops):
        diff = abs(res_me.expect[idx] - res_brme.expect[idx]).max()
        assert_(diff < 5e-2)  # accept 5% error

def test_enr_destory_full():
    "Excitation-number-restricted state-space: full state space"
    a1, a2 = enr_destroy([4, 4], 4**2)
    b1, b2 = tensor(destroy(4), identity(4)), tensor(identity(4), destroy(4))

    assert_(a1 == b1)
    assert_(a2 == b2)

 def test_state_to_state(self):
     """
     control.pulseoptim: state-to-state transfer 
     linear initial pulse used
     assert that goal is achieved
     """       
     # 2 qubits with Ising interaction
     # some arbitary coupling constants
     alpha = [0.9, 0.7]
     beta  = [0.8, 0.9]
     Sx = sigmax()
     Sz = sigmaz()
     H_d = (alpha[0]*tensor(Sx,identity(2)) + 
           alpha[1]*tensor(identity(2),Sx) +
           beta[0]*tensor(Sz,identity(2)) +
           beta[1]*tensor(identity(2),Sz))
     H_c = [tensor(Sz,Sz)]
     
     q1_0 = q2_0 = Qobj([[1], [0]])
     q1_T = q2_T = Qobj([[0], [1]])
     
     psi_0 = tensor(q1_0, q2_0)
     psi_T = tensor(q1_T, q2_T)
     
     n_ts = 10
     evo_time = 18
     
     # Run the optimisation
     result = cpo.optimize_pulse_unitary(H_d, H_c, psi_0, psi_T, 
                     n_ts, evo_time, 
                     fid_err_targ=1e-10, 
                     init_pulse_type='LIN', 
                     gen_stats=True)
     assert_(result.goal_achieved, msg="State-to-state goal not achieved. "
                 "Terminated due to: {}, with infidelity: {}".format(
                 result.termination_reason, result.fid_err))
     assert_almost_equal(result.fid_err, 0.0, decimal=10, 
                         err_msg="Hadamard infidelity too high")
                         
     #Try with Qobj propagation
     result = cpo.optimize_pulse_unitary(H_d, H_c, psi_0, psi_T, 
                     n_ts, evo_time, 
                     fid_err_targ=1e-10, 
                     init_pulse_type='LIN', 
                     dyn_params={'oper_dtype':Qobj},
                     gen_stats=True)
     assert_(result.goal_achieved, msg="State-to-state goal not achieved "
                 "(Qobj propagation)"
                 "Terminated due to: {}, with infidelity: {}".format(
                 result.termination_reason, result.fid_err))

    def test_01_2_unitary_hadamard_no_stats(self):
        """
        control.pulseoptim: Hadamard gate with linear initial pulses (no stats)
        assert that goal is achieved
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 10
        evo_time = 10

        # Run the optimisation
        #Try without stats
        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-10,
                        init_pulse_type='LIN',
                        gen_stats=False)
        assert_(result.goal_achieved, msg="Hadamard goal not achieved "
                                            "(no stats). "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result.termination_reason, result.fid_err))

    def test_01_4_unitary_hadamard_qobj(self):
        """
        control.pulseoptim: Hadamard gate with linear initial pulses (Qobj)
        assert that goal is achieved
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 10
        evo_time = 10

        # Run the optimisation
        #Try with Qobj propagation
        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-10,
                        init_pulse_type='LIN',
                        dyn_params={'oper_dtype':Qobj},
                        gen_stats=True)
        assert_(result.goal_achieved, msg="Hadamard goal not achieved "
                                            "(Qobj propagation). "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result.termination_reason, result.fid_err))

    def test_9_time_dependent_drift(self):
        """
        control.pulseoptim: Hadamard gate with fixed and time varying drift
        assert that goal is achieved for both and that different control
        pulses are produced (only) when they should be
        """
        # Hadamard
        H_0 = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 20
        evo_time = 10
        
        drift_amps_flat = np.ones([n_ts], dtype=float)
        dript_amps_step = [np.round(float(k)/n_ts) for k in range(n_ts)]
        
        # Run the optimisations
        result_fixed = cpo.optimize_pulse_unitary(H_0, H_c, U_0, U_targ, 
                        n_ts, evo_time, 
                        fid_err_targ=1e-10, 
                        init_pulse_type='LIN', 
                        gen_stats=True)
        assert_(result_fixed.goal_achieved, 
                    msg="Fixed drift goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result_fixed.termination_reason, result_fixed.fid_err))
                    
        H_d = [drift_amps_flat[k]*H_0 for k in range(n_ts)]
        result_flat = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ, 
                        n_ts, evo_time, 
                        fid_err_targ=1e-10, 
                        init_pulse_type='LIN', 
                        gen_stats=True)
        assert_(result_flat.goal_achieved, msg="Flat drift goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result_flat.termination_reason, result_flat.fid_err))
                    
        # Check fixed and flat produced the same pulse
        assert_almost_equal(result_fixed.final_amps, result_flat.final_amps, 
                            decimal=9, 
                            err_msg="Flat and fixed drift result in "
                                    "different control pules")
                            
        H_d = [dript_amps_step[k]*H_0 for k in range(n_ts)]
        result_step = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ, 
                        n_ts, evo_time, 
                        fid_err_targ=1e-10, 
                        init_pulse_type='LIN', 
                        gen_stats=True)
        assert_(result_step.goal_achieved, msg="Step drift goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result_step.termination_reason, result_step.fid_err))
                    
        # Check step and flat produced different results
        assert_(np.any(
            np.abs(result_flat.final_amps - result_step.final_amps) > 1e-3), 
                            msg="Flat and step drift result in "
                                    "the same control pules")

    def test_crab(self):
        """
        Optimise pulse for Hadamard gate using CRAB algorithm
        Apply guess and ramping pulse
        assert that goal is achieved and fidelity error is below threshold
        assert that starting amplitude is zero
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 12
        evo_time = 10
        
        # Run the optimisation
        result = cpo.opt_pulse_crab_unitary(H_d, H_c, U_0, U_targ, 
                n_ts, evo_time, 
                fid_err_targ=1e-5, 
                alg_params={'crab_pulse_params':{'randomize_coeffs':False, 
                                                 'randomize_freqs':False}},
                init_coeff_scaling=0.5,
                guess_pulse_type='GAUSSIAN', 
                guess_pulse_params={'variance':0.1*evo_time},
                guess_pulse_scaling=1.0, guess_pulse_offset=1.0,
                amp_lbound=None, amp_ubound=None,
                ramping_pulse_type='GAUSSIAN_EDGE', 
                ramping_pulse_params={'decay_time':evo_time/100.0},
                gen_stats=True)
        assert_(result.goal_achieved, msg="Hadamard goal not achieved")
        assert_almost_equal(result.fid_err, 0.0, decimal=3, 
                            err_msg="Hadamard infidelity too high")
        assert_almost_equal(result.final_amps[0, 0], 0.0, decimal=3, 
                            err_msg="lead in amplitude not zero")

def _pseudo_inverse_dense(L, rhoss, method='direct', **pseudo_args):
    """
    Internal function for computing the pseudo inverse of an Liouvillian using
    dense matrix methods. See pseudo_inverse for details.
    """
    if method == 'direct':
        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
        LIQ = np.linalg.solve(L.full(), Q)
        R = np.dot(Q, LIQ)

        return Qobj(R, dims=L.dims)

    elif method == 'numpy':
        return Qobj(np.linalg.pinv(L.full()), dims=L.dims)

    elif method == 'scipy':
        return Qobj(la.pinv(L.full()), dims=L.dims)

    elif method == 'scipy2':
        return Qobj(la.pinv2(L.full()), dims=L.dims)

    else:
        raise ValueError("Unsupported method '%s'. Use 'direct' or 'numpy'" %
                         method)

 def count_waves(n_ts, evo_time, ptype, freq=None, num_waves=None):
     
     # Any dyn config will do 
     #Hadamard
     H_d = sigmaz()
     H_c = [sigmax()]
     U_0 = identity(2)
     U_targ = hadamard_transform(1)
     
     pulse_params = {}
     if freq is not None:
         pulse_params['freq'] = freq
     if num_waves is not None:
         pulse_params['num_waves'] = num_waves
     
     optim = cpo.create_pulse_optimizer(H_d, H_c, U_0, U_targ, 
                                 n_ts, evo_time, 
                                 dyn_type='UNIT', 
                                 init_pulse_type=ptype,
                                 init_pulse_params=pulse_params,
                                 gen_stats=False)
     pgen = optim.pulse_generator
     pulse = pgen.gen_pulse()
     
     # count number of waves
     zero_cross = pulse[0:-2]*pulse[1:-1] < 0
     
     return (sum(zero_cross) + 1) / 2

    def test_01_1_unitary_hadamard(self):
        """
        control.pulseoptim: Hadamard gate with linear initial pulses
        assert that goal is achieved and fidelity error is below threshold
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 10
        evo_time = 10

        # Run the optimisation
        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-10,
                        init_pulse_type='LIN',
                        gen_stats=True)
        assert_(result.goal_achieved, msg="Hadamard goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result.termination_reason, result.fid_err))
        assert_almost_equal(result.fid_err, 0.0, decimal=10,
                            err_msg="Hadamard infidelity too high")

    def test_01_6_unitary_hadamard_grad(self):
        """
        control.pulseoptim: Hadamard gate gradient check
        assert that gradient approx and exact gradient match in tolerance
        """
        # Hadamard
        H_d = sigmaz()
        H_c = [sigmax()]
        U_0 = identity(2)
        U_targ = hadamard_transform(1)

        n_ts = 10
        evo_time = 10

        # Create the optim objects
        optim = cpo.create_pulse_optimizer(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-10,
                        dyn_type='UNIT',
                        init_pulse_type='LIN',
                        gen_stats=True)
        dyn = optim.dynamics

        init_amps = optim.pulse_generator.gen_pulse().reshape([-1, 1])
        dyn.initialize_controls(init_amps)

        # Check the exact gradient
        func = optim.fid_err_func_wrapper
        grad = optim.fid_err_grad_wrapper
        x0 = dyn.ctrl_amps.flatten()
        grad_diff = check_grad(func, grad, x0)
        assert_almost_equal(grad_diff, 0.0, decimal=6,
                            err_msg="Unitary gradient outside tolerance")

def genstate(s):
    """helper function to obtain the correct initialization"""
    newstr = s.replace("+","0")
    newstr = newstr.replace("-","1")
    print(newstr)
    state = qt.ket(newstr)
    if s[0] == "+" or s[0]=="-":
        operator=qt.hadamard_transform()
    else:
        operator=qt.identity(2)
    for i in range(1,len(s),1):
        if s[i] == "+" or s[i]=="-":
            #apply hadamard
            operator = qt.tensor([operator,qt.hadamard_transform()])
        else:
            operator = qt.tensor([operator,qt.identity(2)])
    return operator * state

    def test_02_1_qft(self):
        """
        control.pulseoptim: QFT gate with linear initial pulses
        assert that goal is achieved and fidelity error is below threshold
        """
        Sx = sigmax()
        Sy = sigmay()
        Sz = sigmaz()
        Si = 0.5*identity(2)

        H_d = 0.5*(tensor(Sx, Sx) + tensor(Sy, Sy) + tensor(Sz, Sz))
        H_c = [tensor(Sx, Si), tensor(Sy, Si), tensor(Si, Sx), tensor(Si, Sy)]
        U_0 = identity(4)
        # Target for the gate evolution - Quantum Fourier Transform gate
        U_targ = qft.qft(2)

        n_ts = 10
        evo_time = 10

        result = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-9,
                        init_pulse_type='LIN',
                        gen_stats=True)

        assert_(result.goal_achieved, msg="QFT goal not achieved. "
                    "Terminated due to: {}, with infidelity: {}".format(
                    result.termination_reason, result.fid_err))
        assert_almost_equal(result.fid_err, 0.0, decimal=7,
                            err_msg="QFT infidelity too high")

        # check bounds
        result2 = cpo.optimize_pulse_unitary(H_d, H_c, U_0, U_targ,
                        n_ts, evo_time,
                        fid_err_targ=1e-9,
                        amp_lbound=-1.0, amp_ubound=1.0,
                        init_pulse_type='LIN',
                        gen_stats=True)
        assert_((result2.final_amps >= -1.0).all() and
                    (result2.final_amps <= 1.0).all(),
                    msg="Amplitude bounds exceeded for QFT")

    def jc_integrate(self, N, wc, wa, g, kappa, gamma,
                     pump, psi0, use_rwa, tlist):

        # Hamiltonian
        a = tensor(destroy(N), identity(2))
        sm = tensor(identity(N), destroy(2))

        if use_rwa:
            # use the rotating wave approxiation
            H = wc * a.dag() * a + wa * sm.dag() * sm + g * (
                a.dag() * sm + a * sm.dag())
        else:
            H = wc * a.dag() * a + wa * sm.dag() * sm + g * (
                a.dag() + a) * (sm + sm.dag())

        # collapse operators
        c_op_list = []

        n_th_a = 0.0  # zero temperature

        rate = kappa * (1 + n_th_a)
        c_op_list.append(np.sqrt(rate) * a)

        rate = kappa * n_th_a
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * a.dag())

        rate = gamma
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sm)

        rate = pump
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sm.dag())

        # evolve and calculate expectation values
        output = mesolve(
            H, psi0, tlist, c_op_list, [a.dag() * a, sm.dag() * sm])
        expt_list = output.expect[0], output.expect[1]
        return expt_list[0], expt_list[1]

def Singlet_yield(theta, rate) :
	(Ha_eigen_theta,H_a_eigenstates, Hb_eigen_theta, H_b_eigenstates) = Hamiltonian_toymodel(theta,0)
	Sx_aq 	 = qt.tensor(Sx, qt.identity(3)).transform(H_a_eigenstates)
	Sy_aq 	 = qt.tensor(Sy, qt.identity(3)).transform(H_a_eigenstates)
	Sz_aq    = qt.tensor(Sz, qt.identity(3)).transform(H_a_eigenstates)
	Sx_bq 	 = qt.tensor(Sx, qt.identity(3)).transform(H_b_eigenstates)
	Sy_bq 	 = qt.tensor(Sy, qt.identity(3)).transform(H_b_eigenstates)
	Sz_bq    = qt.tensor(Sz, qt.identity(3)).transform(H_b_eigenstates)
	Sx_a  	 = Sx_aq.data
	Sy_a	 = Sy_aq.data
	Sz_a     = Sz_aq.data
	Sx_b 	 = Sx_bq.data
	Sy_b 	 = Sy_bq.data
	Sz_b     = Sz_bq.data
	S_a 	 = np.array([Sx_a, Sy_a, Sz_a])
	S_b      = np.array([Sx_b, Sy_b, Sz_b])
	g_a1 = 0.0
	for p in range(0,3):
		for q in range(0,3):
			(rolpA,colpA) = S_a[p].nonzero()
			(rolqA,colqA) = S_a[q].nonzero()
			(rolpB,colpB) = S_b[p].nonzero()
			(rolqB,colqB) = S_b[q].nonzero()
			A = list(set(zip(rolpA, colpA)).intersection(zip(rolqA, colqA))) # Techniques to iterate over sparse matrices
			B = list(set(zip(rolpB, colpB)).intersection(zip(rolqB, colqB))) # Taking Intersction over non-zero values
			#print((10.0*p + q)/100)
			for n,m in A:
				#print (str(n/576.0*100.0))        		
				for r,s in B:
					wa_mn = Ha_eigen_theta[m] - Ha_eigen_theta[n]						# Frequencies as deined in timmel et. al. 1998
					wb_rs = Hb_eigen_theta[r] - Hb_eigen_theta[s]
					g_a1 +=  S_a[p][n, m]*S_a[q][m, n]*S_b[p][r,s]*S_b[q][s, r]* (rate**2/(rate**2 + (wa_mn - wb_rs)**2))				

	singletyield = (g_a1/9.0) + .25
	return singletyield

    def jc_steadystate(self, N, wc, wa, g, kappa, gamma,
                       pump, psi0, use_rwa, tlist):

        # Hamiltonian
        a = tensor(destroy(N), identity(2))
        sm = tensor(identity(N), destroy(2))

        if use_rwa:
            # use the rotating wave approxiation
            H = wc * a.dag(
            ) * a + wa * sm.dag() * sm + g * (a.dag() * sm + a * sm.dag())
        else:
            H = wc * a.dag() * a + wa * sm.dag() * sm + g * (
                a.dag() + a) * (sm + sm.dag())

        # collapse operators
        c_op_list = []

        n_th_a = 0.0  # zero temperature

        rate = kappa * (1 + n_th_a)
        c_op_list.append(np.sqrt(rate) * a)

        rate = kappa * n_th_a
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * a.dag())

        rate = gamma
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sm)

        rate = pump
        if rate > 0.0:
            c_op_list.append(np.sqrt(rate) * sm.dag())

        # find the steady state
        rho_ss = steadystate(H, c_op_list)

        return expect(a.dag() * a, rho_ss), expect(sm.dag() * sm, rho_ss)

def _pseudo_inverse_sparse(L, rhoss, method='splu', **pseudo_args):
    """
    Internal function for computing the pseudo inverse of an Liouvillian using
    sparse matrix methods. See pseudo_inverse for details.
    """

    N = np.prod(L.dims[0][0])

    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 = sp.kron(rhoss_vec.data, tr_op_vec.data.T, format='csr')
    I = sp.eye(N*N, N*N, format='csr')
    Q = I - P

    if pseudo_args['use_rcm']:
        perm = reverse_cuthill_mckee(L.data)
        A = sp_permute(L.data, perm, perm, 'csr')
        Q = sp_permute(Q, perm, perm, 'csr')
    else:
        if not settings.has_mkl:
            A = L.data.tocsc()
        A.sort_indices()
    
    if method == 'splu':
        if settings.has_mkl:
            LIQ = mkl_spsolve(A,Q.toarray())
        else:
            lu = sp.linalg.splu(A, permc_spec=pseudo_args['permc_spec'],
                            diag_pivot_thresh=pseudo_args['diag_pivot_thresh'],
                            options=dict(ILU_MILU=pseudo_args['ILU_MILU']))
            LIQ = lu.solve(Q.toarray())

    elif method == 'spilu':
        lu = sp.linalg.spilu(A, permc_spec=pseudo_args['permc_spec'],
                             fill_factor=pseudo_args['fill_factor'], 
                             drop_tol=pseudo_args['drop_tol'])
        LIQ = lu.solve(Q.toarray())

    else:
        raise ValueError("unsupported method '%s'" % method)

    R = sp.csr_matrix(Q * LIQ)

    if pseudo_args['use_rcm']:
        rev_perm = np.argsort(perm)
        R = sp_permute(R, rev_perm, rev_perm, 'csr')

    return Qobj(R, dims=L.dims)

def _pseudo_inverse_sparse(L, rhoss, method='splu', use_umfpack=False,
                           use_rcm=False):
    """
    Internal function for computing the pseudo inverse of an Liouvillian using
    sparse matrix methods. See pseudo_inverse for details.
    """

    N = np.prod(L.dims[0][0])

    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 = sp.kron(rhoss_vec.data, tr_op_vec.data.T, format='csc')
    I = sp.eye(N*N, N*N, format='csc')
    Q = I - P

    if use_rcm:
        perm = reverse_cuthill_mckee(L.data)
        A = sp_permute(L.data, perm, perm, 'csc').tocsc()
        Q = sp_permute(Q, perm, perm, 'csc')
        permc_spec = 'NATURAL'
    else:
        A = L.data.tocsc()
        A.sort_indices()
        permc_spec = 'COLAMD'

    if method == 'spsolve':
        sp.linalg.use_solver(assumeSortedIndices=True, useUmfpack=use_umfpack)
        LIQ = sp.linalg.spsolve(A, Q)

    elif method == 'splu':
        lu = sp.linalg.splu(A, permc_spec=permc_spec)
        LIQ = lu.solve(Q.toarray())

    elif method == 'spilu':
        lu = sp.linalg.spilu(A, permc_spec=permc_spec,
                             fill_factor=10, drop_tol=1e-8)
        LIQ = lu.solve(Q.toarray())

    else:
        raise ValueError("unsupported method '%s'" % method)

    R = sp.csc_matrix(Q * LIQ)

    if use_rcm:
        rev_perm = np.argsort(perm)
        R = sp_permute(R, rev_perm, rev_perm, 'csc')

    return Qobj(R, dims=L.dims)