1. qutip-notebooks
  2. development

QuTiP development notebook for testing steadystate solvers

Copyright (C) 2013, Paul D. Nation & Robert J. Johansson

In [1]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib
In [2]:
from qutip import *
In [3]:
import time

Setup

In [4]:
reps = 10
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 Exception as e:
            bm_mat[sol_idx] = nan
            print("Failure in %s: %s" % (solver_name, str(e)))

    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, method='direct', use_umfpack=True, sparse=True)],
    ["direct sparse use_umfpack=False", 
     lambda H, c_ops: steadystate(H, c_ops, method='direct', use_umfpack=False, sparse=True)],
    ["direct dense use_umfpack=True",   
     lambda H, c_ops: steadystate(H, c_ops, method='direct', use_umfpack=True, sparse=False)],
    ["direct dense use_umfpack=False",  
     lambda H, c_ops: steadystate(H, c_ops, method='direct', use_umfpack=False, sparse=False)],

    ["iterative use_precond=True, use_rcm=True, sym=True",      
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=True, sym=True)],
    ["iterative use_precond=True, use_rcm=False, sym=True",   
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=False, sym=True)],
    ["iterative use_precond=False, use_rcm=True, sym=True",     
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False, use_rcm=True, sym=True)],
    ["iterative use_precond=False, use_rcm=False, sym=True",     
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False, use_rcm=False, sym=True)],
    ["iterative use_precond=True, use_rcm=True, sym=False",      
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=True, sym=False)],
    ["iterative use_precond=True, use_rcm=False, sym=False",   
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=False, sym=False)],
    ["iterative use_precond=False, use_rcm=True, sym=False",     
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False, use_rcm=True, sym=False)],
    ["iterative use_precond=False, use_rcm=False, sym=False",     
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False, use_rcm=False, sym=False)],

    ["iterative-bicg use_precond=True, use_rcm=True, sym=True",  
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=True, sym=True)],
    ["iterative-bicg use_precond=False, use_rcm=False, sym=True", 
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False, use_rcm=False, sym=True)],
    ["iterative-bicg use_precond=True, use_rcm=True, sym=True",  
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=True, sym=True)],
    ["iterative-bicg use_precond=False, use_rcm=False, sym=True", 
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False, use_rcm=False, sym=True)],
    ["iterative-bicg use_precond=True, use_rcm=True, sym=False",      
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=True, sym=False)],
    ["iterative-bicg use_precond=True, use_rcm=False, sym=False",   
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=False, sym=False)],
    ["iterative-bicg use_precond=False, use_rcm=True, sym=False",     
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False, use_rcm=True, sym=False)],
    ["iterative-bicg use_precond=False, use_rcm=False, sym=False",     
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False, use_rcm=False, sym=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, use_rcm=True, sym=False",      
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=True, sym=False)],
    ["iterative use_precond=True, use_rcm=False, sym=False",      
     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=False, sym=False)],

    ["iterative-bicg use_precond=True, use_rcm=True, sym=False",      
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=True, sym=False)],
    ["iterative-bicg use_precond=True, use_rcm=False, sym=False",      
     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=False, sym=False)],

   ]

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)

Failure in iterative use_precond=False, use_rcm=True, sym=True: Steadystate solver did not reach tolerance after 1000 steps.
Failure in iterative use_precond=False, use_rcm=False, sym=True: Steadystate solver did not reach tolerance after 1000 steps.
Failure in iterative use_precond=False, use_rcm=True, sym=False: Steadystate solver did not reach tolerance after 1000 steps.
Failure in iterative use_precond=False, use_rcm=False, sym=False: Steadystate solver did not reach tolerance after 1000 steps.
/home/rob/py-envs/py3-stable/lib/python3.3/site-packages/qutip/steadystate.py:153: UserWarning: The use of use_umfpack is deprecated.
  warnings.warn("The use of use_umfpack is deprecated.")
-c:3: RuntimeWarning: invalid value encountered in greater









Out[11]:




SolverElapsed timeRatio
power use_umfpack=True0.023376351.47
power use_umfpack=False0.023208381.46
direct sparse use_umfpack=True0.018018441.14
direct sparse use_umfpack=False0.017402011.10
direct dense use_umfpack=True0.048663813.07
direct dense use_umfpack=False0.048664403.07
iterative use_precond=True, use_rcm=True, sym=True0.020518041.29
iterative use_precond=True, use_rcm=False, sym=True0.018939021.19
iterative use_precond=False, use_rcm=True, sym=Truenannan
iterative use_precond=False, use_rcm=False, sym=Truenannan
iterative use_precond=True, use_rcm=True, sym=False0.018178271.15
iterative use_precond=True, use_rcm=False, sym=False0.016614601.05
iterative use_precond=False, use_rcm=True, sym=Falsenannan
iterative use_precond=False, use_rcm=False, sym=Falsenannan
iterative-bicg use_precond=True, use_rcm=True, sym=True0.018684961.18
iterative-bicg use_precond=False, use_rcm=False, sym=True0.068217164.30
iterative-bicg use_precond=True, use_rcm=True, sym=True0.018803641.18
iterative-bicg use_precond=False, use_rcm=False, sym=True0.068448714.31
iterative-bicg use_precond=True, use_rcm=True, sym=False0.018480161.16
iterative-bicg use_precond=True, use_rcm=False, sym=False0.016998171.07
iterative-bicg use_precond=False, use_rcm=True, sym=False0.065573764.13
iterative-bicg use_precond=False, use_rcm=False, sym=False0.067882544.28
lu0.015870191.00
















In [12]:






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















Out[12]:




SolverElapsed timeRatio
power use_umfpack=True0.278076221.19
power use_umfpack=False0.280875871.20
direct sparse use_umfpack=True0.234426811.00
direct sparse use_umfpack=False0.283161661.21
iterative use_precond=True, use_rcm=True, sym=False0.271456861.16
iterative use_precond=True, use_rcm=False, sym=False0.273619631.17
iterative-bicg use_precond=True, use_rcm=True, sym=False0.273963981.17
iterative-bicg use_precond=True, use_rcm=False, sym=False0.298872971.27


















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)



















Failure in iterative-bicg use_precond=True, use_rcm=False, sym=False: Steadystate solver failed with fatal error: -10.









Out[14]:




SolverElapsed timeRatio
power use_umfpack=True0.021985101.60
power use_umfpack=False0.021958161.60
direct sparse use_umfpack=True0.017395211.27
direct sparse use_umfpack=False0.015027051.09
direct dense use_umfpack=True5.38591220391.95
direct dense use_umfpack=False5.56329439404.86
iterative use_precond=True, use_rcm=True, sym=True0.019405961.41
iterative use_precond=True, use_rcm=False, sym=True0.018778871.37
iterative use_precond=False, use_rcm=True, sym=True0.8553290162.25
iterative use_precond=False, use_rcm=False, sym=True0.9050169065.86
iterative use_precond=True, use_rcm=True, sym=False0.019192601.40
iterative use_precond=True, use_rcm=False, sym=False0.018242191.33
iterative use_precond=False, use_rcm=True, sym=False0.8232050759.91
iterative use_precond=False, use_rcm=False, sym=False0.8994701965.46
iterative-bicg use_precond=True, use_rcm=True, sym=True0.018724131.36
iterative-bicg use_precond=False, use_rcm=False, sym=True0.036521552.66
iterative-bicg use_precond=True, use_rcm=True, sym=True0.018840171.37
iterative-bicg use_precond=False, use_rcm=False, sym=True0.035947592.62
iterative-bicg use_precond=True, use_rcm=True, sym=False0.019165831.39
iterative-bicg use_precond=True, use_rcm=False, sym=Falsenannan
iterative-bicg use_precond=False, use_rcm=True, sym=False0.039077972.84
iterative-bicg use_precond=False, use_rcm=False, sym=False0.035814982.61
lu0.013741281.00


















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]:




SolverElapsed timeRatio
power use_umfpack=True33.7502831710.41
power use_umfpack=False34.6520989410.68
direct sparse use_umfpack=True8.380447082.58
direct sparse use_umfpack=False27.644433438.52
iterative use_precond=True, use_rcm=True, sym=False3.243287921.00
iterative use_precond=True, use_rcm=False, sym=False3.474131251.07
iterative-bicg use_precond=True, use_rcm=True, sym=False4.313132001.33
iterative-bicg use_precond=True, use_rcm=False, sym=False3.444228721.06


















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]:




SolverElapsed timeRatio
power use_umfpack=True0.011187241.95
power use_umfpack=False0.010905841.90
direct sparse use_umfpack=True0.006970671.21
direct sparse use_umfpack=False0.006574891.14
direct dense use_umfpack=True0.005773351.00
direct dense use_umfpack=False0.005751781.00
iterative use_precond=True, use_rcm=True, sym=True0.007394221.29
iterative use_precond=True, use_rcm=False, sym=True0.006448941.12
iterative use_precond=False, use_rcm=True, sym=True0.007358721.28
iterative use_precond=False, use_rcm=False, sym=True0.006387951.11
iterative use_precond=True, use_rcm=True, sym=False0.007451421.30
iterative use_precond=True, use_rcm=False, sym=False0.006461141.12
iterative use_precond=False, use_rcm=True, sym=False0.007483701.30
iterative use_precond=False, use_rcm=False, sym=False0.006316521.10
iterative-bicg use_precond=True, use_rcm=True, sym=True0.007402991.29
iterative-bicg use_precond=False, use_rcm=False, sym=True0.006339241.10
iterative-bicg use_precond=True, use_rcm=True, sym=True0.007311631.27
iterative-bicg use_precond=False, use_rcm=False, sym=True0.006341411.10
iterative-bicg use_precond=True, use_rcm=True, sym=False0.007325631.27
iterative-bicg use_precond=True, use_rcm=False, sym=False0.006478021.13
iterative-bicg use_precond=False, use_rcm=True, sym=False0.007355171.28
iterative-bicg use_precond=False, use_rcm=False, sym=False0.006371311.11
lu0.006368351.11


















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]:




SolverElapsed timeRatio
power use_umfpack=True0.014272051.59
power use_umfpack=False0.014378831.60
direct sparse use_umfpack=True0.010162311.13
direct sparse use_umfpack=False0.009893941.10
direct dense use_umfpack=True0.009034371.00
direct dense use_umfpack=False0.009002691.00
iterative use_precond=True, use_rcm=True, sym=True0.010799811.20
iterative use_precond=True, use_rcm=False, sym=True0.009700681.08
iterative use_precond=False, use_rcm=True, sym=True0.010652611.18
iterative use_precond=False, use_rcm=False, sym=True0.009600471.07
iterative use_precond=True, use_rcm=True, sym=False0.010772181.20
iterative use_precond=True, use_rcm=False, sym=False0.009826731.09
iterative use_precond=False, use_rcm=True, sym=False0.010667281.18
iterative use_precond=False, use_rcm=False, sym=False0.009759401.08
iterative-bicg use_precond=True, use_rcm=True, sym=True0.010679201.19
iterative-bicg use_precond=False, use_rcm=False, sym=True0.009984521.11
iterative-bicg use_precond=True, use_rcm=True, sym=True0.010668231.19
iterative-bicg use_precond=False, use_rcm=False, sym=True0.010011031.11
iterative-bicg use_precond=True, use_rcm=True, sym=False0.010837411.20
iterative-bicg use_precond=True, use_rcm=False, sym=False0.009566211.06
iterative-bicg use_precond=False, use_rcm=True, sym=False0.011081551.23
iterative-bicg use_precond=False, use_rcm=False, sym=False0.009983401.11
lu0.009497951.06
















In [21]:






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



















Failure in iterative use_precond=False, use_rcm=True, sym=True: Steadystate solver did not reach tolerance after 1000 steps.
Failure in iterative use_precond=False, use_rcm=False, sym=True: Steadystate solver did not reach tolerance after 1000 steps.
Failure in iterative use_precond=False, use_rcm=True, sym=False: Steadystate solver did not reach tolerance after 1000 steps.
Failure in iterative use_precond=False, use_rcm=False, sym=False: Steadystate solver did not reach tolerance after 1000 steps.









Out[21]:




SolverElapsed timeRatio
power use_umfpack=True0.024710321.25
power use_umfpack=False0.024635101.25
direct sparse use_umfpack=True0.019916391.01
direct sparse use_umfpack=False0.019714621.00
direct dense use_umfpack=True0.026166921.33
direct dense use_umfpack=False0.026031971.32
iterative use_precond=True, use_rcm=True, sym=True0.022794911.16
iterative use_precond=True, use_rcm=False, sym=True0.020919991.06
iterative use_precond=False, use_rcm=True, sym=Truenannan
iterative use_precond=False, use_rcm=False, sym=Truenannan
iterative use_precond=True, use_rcm=True, sym=False0.023444911.19
iterative use_precond=True, use_rcm=False, sym=False0.021756481.10
iterative use_precond=False, use_rcm=True, sym=Falsenannan
iterative use_precond=False, use_rcm=False, sym=Falsenannan
iterative-bicg use_precond=True, use_rcm=True, sym=True0.022286841.13
iterative-bicg use_precond=False, use_rcm=False, sym=True0.077180293.92
iterative-bicg use_precond=True, use_rcm=True, sym=True0.022439651.14
iterative-bicg use_precond=False, use_rcm=False, sym=True0.077550053.94
iterative-bicg use_precond=True, use_rcm=True, sym=False0.022431281.14
iterative-bicg use_precond=True, use_rcm=False, sym=False0.020930191.06
iterative-bicg use_precond=False, use_rcm=True, sym=False0.076412323.88
iterative-bicg use_precond=False, use_rcm=False, sym=False0.079334784.03
lu0.019693541.00
















In [22]:






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



















Failure in iterative use_precond=True, use_rcm=True, sym=False: Steadystate solver did not reach tolerance after 1000 steps.
Failure in iterative use_precond=True, use_rcm=False, sym=False: Steadystate solver did not reach tolerance after 1000 steps.
Failure in iterative-bicg use_precond=True, use_rcm=True, sym=False: Steadystate solver did not reach tolerance after 40960 steps.
Failure in iterative-bicg use_precond=True, use_rcm=False, sym=False: Steadystate solver did not reach tolerance after 40960 steps.









Out[22]:




SolverElapsed timeRatio
power use_umfpack=True10.343706238.90
power use_umfpack=False10.375724248.93
direct sparse use_umfpack=True1.161624791.00
direct sparse use_umfpack=False9.754040058.40
iterative use_precond=True, use_rcm=True, sym=Falsenannan
iterative use_precond=True, use_rcm=False, sym=Falsenannan
iterative-bicg use_precond=True, use_rcm=True, sym=Falsenannan
iterative-bicg use_precond=True, use_rcm=False, sym=Falsenannan


















Software versions












In [23]:






from qutip.ipynbtools import version_table

version_table()















Out[23]:




SoftwareVersion
Python3.3.2+ (default, Oct 9 2013, 14:50:09)
[GCC 4.8.1]
matplotlib1.3.1
QuTiP3.0.0.dev-0f20b2a
Cython0.20.1
SciPy0.13.3
Numpy1.8.0
IPython2.0.0
OSposix [linux]
Thu Apr 03 17:04:46 2014 JST
















In [23]:






 















In [23]: