QuTiP development notebook for testing steadystate solvers



In [1]:

    
%pylab inline









    



Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].
For more information, type 'help(pylab)'.



In [2]:

    
from qutip import *



In [3]:

    
import time

Setup



In [4]:

    
reps = 1



In [5]:

    
from IPython.display import HTML



In [6]:

    
def show_bm_mat(bm_mat, solvers):
    
    m = bm_mat[bm_mat > 0].min()
    
    html = "<table>"

    html += "<tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr>"
    
    for idx, (desc, func) in enumerate(solvers):
    
        if bm_mat[idx] == m:
            html += "<tr><td><b>%s</b></td><td><b>%.8f</b></td><td><b>%.2f</b></td></tr>" % \
                (desc, bm_mat[idx], bm_mat[idx]/m)
        else:
            html += "<tr><td>%s</td><td>%.8f</td><td>%.2f</td></tr>" % \
                (desc, bm_mat[idx], bm_mat[idx]/m)
            
    
    html += "</table>"

    return HTML(html)



In [7]:

    
def benchmark_steadystate_solvers(args, solvers, problem_func):

    bm_mat = zeros(len(solvers))

    H, c_ops = problem_func(args)
    
    for sol_idx, solver in enumerate(solvers):
        solver_name, solver_func = solver
        try:
            t1 = time.time()
            for r in range(reps):
                rhoss = solver_func(H, c_ops)
            t2 = time.time()
            bm_mat[sol_idx] = (t2 - t1)/reps
    
        except:
            bm_mat[sol_idx] = nan
            
    return bm_mat



In [8]:

    
solvers = [["power use_umfpack=True",          
            lambda H, c_ops: steadystate(H, c_ops, method='power', use_umfpack=True)],
           ["power use_umfpack=False",         
            lambda H, c_ops: steadystate(H, c_ops, method='power', use_umfpack=False)],
           ["direct sparse use_umfpack=True",  
            lambda H, c_ops: steadystate(H, c_ops, use_umfpack=True, sparse=True)],
           ["direct sparse use_umfpack=False", 
            lambda H, c_ops: steadystate(H, c_ops, use_umfpack=False, sparse=True)],
           ["iterative use_precond=True",      
            lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True)],
           ["iterative use_precond=False",     
            lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False)],
           ["iterative-bicg use_precond=True",  
            lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True)],
           ["iterative-bicg use_precond=False", 
            lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False)],
           ["direct dense use_umfpack=True",   
            lambda H, c_ops: steadystate(H, c_ops, use_umfpack=True, sparse=False)],
           ["direct dense use_umfpack=False",  
            lambda H, c_ops: steadystate(H, c_ops, use_umfpack=False, sparse=False)],
           ["svd_dense",                       
            lambda H, c_ops: steadystate(H, c_ops, method='svd')],
           ["lu",                              
            lambda H, c_ops: steadystate(H, c_ops, method='lu')],
          ]



In [9]:

    
large_solvers = [
    ["power use_umfpack=True",          
     lambda H, c_ops: steadystate(H, c_ops, method='power', use_umfpack=True)],
    ["power use_umfpack=False",         
     lambda H, c_ops: steadystate(H, c_ops, method='power', use_umfpack=False)],
    ["direct sparse use_umfpack=True",  
     lambda H, c_ops: steadystate(H, c_ops, use_umfpack=True, sparse=True)],
    ["direct sparse use_umfpack=False", 
     lambda H, c_ops: steadystate(H, c_ops, use_umfpack=False, sparse=True)],
    ["iterative use_precond=True",      
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True)],
    ["iterative-bicg use_precond=True",  
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True)],
   ]

Test problem 1



In [10]:

    
def bm_problem1(N):

    a = tensor(destroy(N), identity(2))
    b = tensor(identity(N), destroy(2))

    H = a.dag() * a + b.dag() * b + 0.25 * (a + a.dag()) * (b + b.dag())

    c_ops = [sqrt(0.1) * a, sqrt(0.075) * a.dag(), sqrt(0.1) * b]
    
    return H, c_ops



In [11]:

    
bm_mat = benchmark_steadystate_solvers(10, solvers, bm_problem1)
show_bm_mat(bm_mat, solvers)









    Out[11]:




Solver Elapsed time Ratio
power use_umfpack=True 0.02930117 1.86
power use_umfpack=False 0.02769184 1.76
direct sparse use_umfpack=True 0.01657772 1.05
direct sparse use_umfpack=False 0.01577711 1.00
iterative use_precond=True 0.01887131 1.20
iterative use_precond=False 0.42715192 27.07
iterative-bicg use_precond=True 0.01824522 1.16
iterative-bicg use_precond=False 0.05967665 3.78
direct dense use_umfpack=True 0.03680706 2.33
direct dense use_umfpack=False 0.03597355 2.28
svd_dense 0.43549585 27.60
lu 0.01581359 1.00



In [12]:

    
bm_mat = benchmark_steadystate_solvers(50, large_solvers, bm_problem1)
show_bm_mat(bm_mat, large_solvers)









    Out[12]:




Solver Elapsed time Ratio
power use_umfpack=True 0.30860066 2.16
power use_umfpack=False 0.30906749 2.16
direct sparse use_umfpack=True 0.26223445 1.83
direct sparse use_umfpack=False 0.30435157 2.13
iterative use_precond=True 0.14319158 1.00
iterative-bicg use_precond=True 0.14433670 1.01

Test problem 2: high temperature harmonic oscillator



In [13]:

    
def bm_problem2(N):
    
    a = destroy(N)
    H = a.dag() * a
    nth = N / 4
    gamma = 0.05
    c_ops = [sqrt(gamma * (nth + 1)) * a, sqrt(gamma * nth) * a.dag()]

    return H, c_ops



In [14]:

    
bm_mat = benchmark_steadystate_solvers(50, solvers, bm_problem2)
show_bm_mat(bm_mat, solvers)









    Out[14]:




Solver Elapsed time Ratio
power use_umfpack=True 0.03091383 2.33
power use_umfpack=False 0.02793431 2.11
direct sparse use_umfpack=True 0.01574659 1.19
direct sparse use_umfpack=False 0.01326680 1.00
iterative use_precond=True 0.01814747 1.37
iterative use_precond=False 0.08728552 6.58
iterative-bicg use_precond=True 0.01731229 1.30
iterative-bicg use_precond=False 0.03140640 2.37
direct dense use_umfpack=True 3.13456941 236.27
direct dense use_umfpack=False 3.13574004 236.36
svd_dense 71.83850408 5414.91
lu 0.01391077 1.05

Test problem 3: Coupled oscillators



In [15]:

    
def bm_problem3(N):
    
    a = tensor(destroy(N), identity(N))
    b = tensor(identity(N), destroy(N))
    
    H = a.dag() * a + 0.25 * b.dag() * b + 0.05 * a.dag() * a * (b + b.dag()) + 0.1 * (a + a.dag()) 

    c_ops = [sqrt(0.05) * a, sqrt(0.015) * a.dag(), sqrt(0.1) * b, sqrt(0.075) * b.dag()]

    return H, c_ops



In [16]:

    
bm_mat = benchmark_steadystate_solvers(10, large_solvers, bm_problem3)
show_bm_mat(bm_mat, large_solvers)









    Out[16]:




Solver Elapsed time Ratio
power use_umfpack=True 5.09422636 7.99
power use_umfpack=False 35.79435349 56.16
direct sparse use_umfpack=True 5.80584502 9.11
direct sparse use_umfpack=False 30.20628524 47.40
iterative use_precond=True 0.71662021 1.12
iterative-bicg use_precond=True 0.63732290 1.00

Test problem 4: a two level system



In [17]:

    
def bm_problem4(args=None):

    sz = sigmaz()    
    sx = sigmax()
    
    H = sz
    c_ops = [sqrt(0.05) * sx]

    return H, c_ops



In [18]:

    
bm_mat = benchmark_steadystate_solvers(None, solvers, bm_problem4)
show_bm_mat(bm_mat, solvers)









    Out[18]:




Solver Elapsed time Ratio
power use_umfpack=True 0.01852179 3.05
power use_umfpack=False 0.01754975 2.89
direct sparse use_umfpack=True 0.00707531 1.17
direct sparse use_umfpack=False 0.00687790 1.13
iterative use_precond=True 0.00709677 1.17
iterative use_precond=False 0.00711107 1.17
iterative-bicg use_precond=True 0.01673436 2.76
iterative-bicg use_precond=False 0.00678062 1.12
direct dense use_umfpack=True 0.00611544 1.01
direct dense use_umfpack=False 0.00606394 1.00
svd_dense 0.00769496 1.27
lu 0.00651836 1.07

Test problem 5: spin chain



In [19]:

    
def bm_problem5(N=1):

    H = 0
    for m in range(N):
        H += tensor([sigmaz() if n == m else identity(2) for n in range(N)])

    for m in range(N-1):
        H += tensor([sigmax() if n in [m,m+1] else identity(2) for n in range(N)])      
        
    c_ops = [sqrt(0.05) * tensor([sigmam() if n == m else identity(2) for n in range(N)])
             for m in range(N)]
   
    return H, c_ops



In [20]:

    
bm_mat = benchmark_steadystate_solvers(2, solvers, bm_problem5)
show_bm_mat(bm_mat, solvers)









    Out[20]:




Solver Elapsed time Ratio
power use_umfpack=True 0.02150130 2.35
power use_umfpack=False 0.02094078 2.28
direct sparse use_umfpack=True 0.01023006 1.12
direct sparse use_umfpack=False 0.00997591 1.09
iterative use_precond=True 0.01020479 1.11
iterative use_precond=False 0.01033449 1.13
iterative-bicg use_precond=True 0.00971484 1.06
iterative-bicg use_precond=False 0.01010346 1.10
direct dense use_umfpack=True 0.00918984 1.00
direct dense use_umfpack=False 0.00916553 1.00
svd_dense 0.01097131 1.20
lu 0.00957203 1.04



In [21]:

    
bm_mat = benchmark_steadystate_solvers(4, solvers, bm_problem5)
show_bm_mat(bm_mat, solvers)









    Out[21]:




Solver Elapsed time Ratio
power use_umfpack=True 0.03119278 1.59
power use_umfpack=False 0.03097486 1.58
direct sparse use_umfpack=True 0.01957011 1.00
direct sparse use_umfpack=False 0.01961446 1.00
iterative use_precond=True 0.02565002 1.31
iterative use_precond=False 0.27903152 14.26
iterative-bicg use_precond=True 0.02416730 1.23
iterative-bicg use_precond=False 0.06678843 3.41
direct dense use_umfpack=True 0.02458715 1.26
direct dense use_umfpack=False 0.02418828 1.24
svd_dense 0.16517234 8.44
lu 0.01960063 1.00



In [22]:

    
bm_mat = benchmark_steadystate_solvers(6, large_solvers, bm_problem5)
show_bm_mat(bm_mat, large_solvers)









    



/usr/local/lib/python3.3/dist-packages/qutip/steadystate.py:314: UserWarning: Preconditioning failed. Continuing without.
  UserWarning)






    Out[22]:




Solver Elapsed time Ratio
power use_umfpack=True 1.09917188 1.08
power use_umfpack=False 9.79997540 9.60
direct sparse use_umfpack=True 1.02135420 1.00
direct sparse use_umfpack=False 9.18890977 9.00
iterative use_precond=True 13.48026681 13.20
iterative-bicg use_precond=True 2.67135262 2.62

Software versions



In [23]:

    
from qutip.ipynbtools import version_table

version_table()









    Out[23]:




Software Version
OS posix [linux]
IPython 1.0.dev
Numpy 1.8.0.dev-c5694c5
QuTiP 2.3.0.dev-169f358
Cython 0.19.1
SciPy 0.13.0.dev-38faf7c
Python 3.3.1 (default, Apr 17 2013, 22:30:32) 
[GCC 4.7.3]
matplotlib 1.4.x
Fri Jul 19 14:08:27 2013 JST

Solver	Elapsed time	Ratio
power use_umfpack=True	0.02930117	1.86
power use_umfpack=False	0.02769184	1.76
direct sparse use_umfpack=True	0.01657772	1.05
direct sparse use_umfpack=False	0.01577711	1.00
iterative use_precond=True	0.01887131	1.20
iterative use_precond=False	0.42715192	27.07
iterative-bicg use_precond=True	0.01824522	1.16
iterative-bicg use_precond=False	0.05967665	3.78
direct dense use_umfpack=True	0.03680706	2.33
direct dense use_umfpack=False	0.03597355	2.28
svd_dense	0.43549585	27.60
lu	0.01581359	1.00

Software	Version
OS	posix [linux]
IPython	1.0.dev
Numpy	1.8.0.dev-c5694c5
QuTiP	2.3.0.dev-169f358
Cython	0.19.1
SciPy	0.13.0.dev-38faf7c
Python	3.3.1 (default, Apr 17 2013, 22:30:32) [GCC 4.7.3]
matplotlib	1.4.x
Fri Jul 19 14:08:27 2013 JST