Copyright (C) 2013, Paul D. Nation & Robert J. Johansson
%pylab inline
Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline]. For more information, type 'help(pylab)'.
from qutip import *
import time
reps = 1
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')
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
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')],
]
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)]
]
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
N_vec = array([2, 5, 10, 15, 20, 25, 30, 35])
bm_mat = benchmark_steadystate_solvers(N_vec, solvers, bm_problem1)
plot_bm_mat(bm_mat, 2 * N_vec, solvers)
N_vec = arange(25, 400, 25)
bm_mat = benchmark_steadystate_solvers(N_vec, large_solvers, bm_problem1)
plot_bm_mat(bm_mat, 2 * N_vec, large_solvers)
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
N_vec = array([2, 5, 10, 15, 20, 25])
bm_mat = benchmark_steadystate_solvers(N_vec, solvers, bm_problem2)
plot_bm_mat(bm_mat, N_vec, solvers)
N_vec = arange(25, 300, 25)
bm_mat = benchmark_steadystate_solvers(N_vec, large_solvers, bm_problem2)
plot_bm_mat(bm_mat, N_vec, large_solvers)
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
N_vec = array([2, 4, 6, 8, 10])
bm_mat = benchmark_steadystate_solvers(N_vec, large_solvers, bm_problem3)
plot_bm_mat(bm_mat, N_vec ** 2, large_solvers)
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
N_vec = array([2, 3, 4, 5, 6])
bm_mat = benchmark_steadystate_solvers(N_vec, large_solvers, bm_problem3)
plot_bm_mat(bm_mat, N_vec, large_solvers)
H, c_ops = bm_problem1(300)
%%prun -q -T steadystate.txt
steadystate(H, c_ops)
*** Profile printout saved to text file 'steadystate.txt'.
!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}
%%prun -q -T steadystate_direct.txt
steadystate(H, c_ops, method="direct")
*** Profile printout saved to text file 'steadystate_direct.txt'.
!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}
%%prun -q -T steadystate_iterative.txt
steadystate(H, c_ops, method="iterative")
*** Profile printout saved to text file 'steadystate_iterative.txt'.
!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}
%%prun -q -T steadystate_iterative.txt
steadystate(H, c_ops, method="iterative-bicg")
*** Profile printout saved to text file 'steadystate_iterative.txt'.
!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)
from qutip.ipynbtools import version_table
version_table()
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 |