QuTiP development notebook for testing and benchmarking 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


    

In [4]:

    

reps = 1


    

In [5]:

    

def plot_bm_mat(bm_mat, N_vec, solvers):
    
    fig, ax = subplots(figsize=(12, 8))
    
    for idx, solver in enumerate(solvers):
        solver_name, solver_func = solver
        ax.plot(N_vec, bm_mat[:, idx], label=solver_name, linewidth=2)
        
    ax.set_yscale('log')
    ax.set_xlabel("Hilbert space size (N)", fontsize=18)
    ax.set_ylabel("Elapsed time", fontsize=18)
    ax.legend(loc=0)
    ax.axis('tight')


    

In [6]:

    

def benchmark_steadystate_solvers(N_vec, solvers, problem_func):

    bm_mat = zeros((len(N_vec), len(solvers)))

    for N_idx, N in enumerate(N_vec):
        H, c_ops = problem_func(N)
    
        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[N_idx, sol_idx] = (t2 - t1)/reps
        
            except:
                bm_mat[N_idx, sol_idx] = 0
            
    return bm_mat


    

In [7]:

    

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

    

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


    

In [9]:

    

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

    

N_vec = array([2, 5, 10, 15, 20, 25, 30, 35])


    

In [11]:

    

bm_mat = benchmark_steadystate_solvers(N_vec, solvers, bm_problem1)


    

In [12]:

    

plot_bm_mat(bm_mat, 2 * N_vec, solvers)


    

In [13]:

    

N_vec = arange(25, 400, 25)


    

In [14]:

    

bm_mat = benchmark_steadystate_solvers(N_vec, large_solvers, bm_problem1)


    

In [15]:

    

plot_bm_mat(bm_mat, 2 * N_vec, large_solvers)


    

In [16]:

    

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

    

N_vec = array([2, 5, 10, 15, 20, 25])


    

In [18]:

    

bm_mat = benchmark_steadystate_solvers(N_vec, solvers, bm_problem2)


    

In [19]:

    

plot_bm_mat(bm_mat, N_vec, solvers)


    

In [20]:

    

N_vec = arange(25, 300, 25)


    

In [21]:

    

bm_mat = benchmark_steadystate_solvers(N_vec, large_solvers, bm_problem2)


    

In [22]:

    

plot_bm_mat(bm_mat, N_vec, large_solvers)


    

In [23]:

    

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

    

N_vec = array([2, 4, 6, 8, 10])


    

In [25]:

    

bm_mat = benchmark_steadystate_solvers(N_vec, large_solvers, bm_problem3)


    

In [26]:

    

plot_bm_mat(bm_mat, N_vec ** 2, large_solvers)


    

In [27]:

    

def bm_problem4(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 [28]:

    

N_vec = array([2, 3, 4, 5, 6])


    

In [29]:

    

bm_mat = benchmark_steadystate_solvers(N_vec, large_solvers, bm_problem3)


    

In [30]:

    

plot_bm_mat(bm_mat, N_vec, large_solvers)


    

In [31]:

    

H, c_ops = bm_problem1(300)


    

In [32]:

    

%%prun -q -T steadystate.txt

steadystate(H, c_ops)


    

    




 
*** Profile printout saved to text file 'steadystate.txt'. 

In [33]:

    

!head -n 20 steadystate.txt


    

    




         9309 function calls (9275 primitive calls) in 59.146 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1   54.735   54.735   54.735   54.735 {built-in method umfpack_zi_numeric}
        1    1.624    1.624    1.624    1.624 {built-in method umfpack_zi_solve}
        1    1.518    1.518    1.518    1.518 {built-in method umfpack_zi_symbolic}
        7    0.277    0.040    0.277    0.040 {built-in method csr_minus_csr}
        1    0.170    0.170    0.170    0.170 {built-in method umfpack_zi_free_numeric}
        4    0.156    0.039    0.156    0.039 {built-in method csr_plus_csr}
       68    0.151    0.002    0.151    0.002 {method 'copy' of 'numpy.ndarray' objects}
       11    0.094    0.009    0.348    0.032 construct.py:272(kron)
       14    0.092    0.007    0.092    0.007 {built-in method coo_tocsr}
       41    0.061    0.001    0.061    0.001 {method 'repeat' of 'numpy.ndarray' objects}
       14    0.040    0.003    0.040    0.003 {built-in method csr_sum_duplicates}
        7    0.036    0.005    0.074    0.011 data.py:72(_mul_scalar)
        1    0.029    0.029    0.856    0.856 superoperator.py:59(liouvillian_fast)
      245    0.029    0.000    0.029    0.000 {method 'reduce' of 'numpy.ufunc' objects}
       17    0.026    0.002    0.026    0.002 {method 'astype' of 'numpy.ndarray' objects}

In [34]:

    

%%prun -q -T steadystate_direct.txt

steadystate(H, c_ops, method="direct")


    

    




 
*** Profile printout saved to text file 'steadystate_direct.txt'. 

In [35]:

    

!head -n 20 steadystate_direct.txt


    

    




         9309 function calls (9275 primitive calls) in 58.058 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1   53.669   53.669   53.669   53.669 {built-in method umfpack_zi_numeric}
        1    1.617    1.617    1.617    1.617 {built-in method umfpack_zi_solve}
        1    1.506    1.506    1.506    1.506 {built-in method umfpack_zi_symbolic}
        7    0.278    0.040    0.278    0.040 {built-in method csr_minus_csr}
        1    0.169    0.169    0.169    0.169 {built-in method umfpack_zi_free_numeric}
        4    0.154    0.038    0.154    0.038 {built-in method csr_plus_csr}
       68    0.135    0.002    0.135    0.002 {method 'copy' of 'numpy.ndarray' objects}
       11    0.096    0.009    0.362    0.033 construct.py:272(kron)
       14    0.095    0.007    0.095    0.007 {built-in method coo_tocsr}
       41    0.067    0.002    0.067    0.002 {method 'repeat' of 'numpy.ndarray' objects}
       14    0.039    0.003    0.039    0.003 {built-in method csr_sum_duplicates}
        7    0.037    0.005    0.075    0.011 data.py:72(_mul_scalar)
      245    0.032    0.000    0.032    0.000 {method 'reduce' of 'numpy.ufunc' objects}
        1    0.028    0.028    0.869    0.869 superoperator.py:59(liouvillian_fast)
       17    0.026    0.002    0.026    0.002 {method 'astype' of 'numpy.ndarray' objects}

In [36]:

    

%%prun -q -T steadystate_iterative.txt

steadystate(H, c_ops, method="iterative")


    

    




 
*** Profile printout saved to text file 'steadystate_iterative.txt'. 

In [37]:

    

!head -n 20 steadystate_iterative.txt


    

    




         10047 function calls (10015 primitive calls) in 7.414 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    3.391    3.391    3.391    3.391 {built-in method gstrf}
       17    1.577    0.093    1.577    0.093 {method 'solve' of 'factored_lu' objects}
        1    1.024    1.024    2.998    2.998 lgmres.py:20(lgmres)
       18    0.344    0.019    0.344    0.019 {built-in method csc_matvec}
        7    0.278    0.040    0.278    0.040 {built-in method csr_minus_csr}
        3    0.118    0.039    0.118    0.039 {built-in method csr_plus_csr}
       11    0.095    0.009    0.355    0.032 construct.py:272(kron)
       13    0.089    0.007    0.089    0.007 {built-in method coo_tocsr}
       41    0.067    0.002    0.067    0.002 {method 'repeat' of 'numpy.ndarray' objects}
        1    0.047    0.047    0.047    0.047 {built-in method csc_plus_csc}
        5    0.042    0.008    0.042    0.008 {built-in method csr_tocsc}
       13    0.039    0.003    0.039    0.003 {built-in method csr_sum_duplicates}
        7    0.037    0.005    0.075    0.011 data.py:72(_mul_scalar)
        1    0.030    0.030    6.538    6.538 steadystate.py:342(_steadystate_iterative)
      249    0.029    0.000    0.029    0.000 {method 'reduce' of 'numpy.ufunc' objects}

In [38]:

    

%%prun -q -T steadystate_iterative.txt

steadystate(H, c_ops, method="iterative-bicg")


    

    




 
*** Profile printout saved to text file 'steadystate_iterative.txt'. 

In [39]:

    

!head -n 20 steadystate_iterative.txt


    

    




         9150 function calls (9118 primitive calls) in 7.102 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    3.393    3.393    3.393    3.393 {built-in method gstrf}
       19    1.727    0.091    1.727    0.091 {method 'solve' of 'factored_lu' objects}
        1    0.567    0.567    2.702    2.702 iterative.py:156(bicgstab)
       19    0.357    0.019    0.357    0.019 {built-in method csc_matvec}
        7    0.274    0.039    0.274    0.039 {built-in method csr_minus_csr}
        3    0.118    0.039    0.118    0.039 {built-in method csr_plus_csr}
       11    0.096    0.009    0.355    0.032 construct.py:272(kron)
       13    0.090    0.007    0.090    0.007 {built-in method coo_tocsr}
       41    0.066    0.002    0.066    0.002 {method 'repeat' of 'numpy.ndarray' objects}
        1    0.047    0.047    0.047    0.047 {built-in method csc_plus_csc}
       26    0.045    0.002    0.045    0.002 {built-in method zeros}
        5    0.042    0.008    0.042    0.008 {built-in method csr_tocsc}
       13    0.039    0.003    0.039    0.003 {built-in method csr_sum_duplicates}
        7    0.037    0.005    0.075    0.011 data.py:72(_mul_scalar)
        1    0.028    0.028    0.859    0.859 superoperator.py:59(liouvillian_fast)

In [40]:

    

from qutip.ipynbtools import version_table

version_table()


    

    
Out[40]:




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