{ "metadata": { "name": "development-steadystate-solver-benchmarks-2" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# QuTiP development notebook for testing steadystate solvers\n", "\n", "Copyright (C) 2013, Paul D. Nation & Robert J. Johansson" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%pylab inline" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "\n", "Welcome to pylab, a matplotlib-based Python environment [backend: module://IPython.kernel.zmq.pylab.backend_inline].\n", "For more information, type 'help(pylab)'.\n" ] } ], "prompt_number": 1 }, { "cell_type": "code", "collapsed": false, "input": [ "from qutip import *" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "code", "collapsed": false, "input": [ "import time" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setup" ] }, { "cell_type": "code", "collapsed": false, "input": [ "reps = 1" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 4 }, { "cell_type": "code", "collapsed": false, "input": [ "from IPython.display import HTML" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "code", "collapsed": false, "input": [ "def show_bm_mat(bm_mat, solvers):\n", " \n", " m = bm_mat[bm_mat > 0].min()\n", " \n", " html = \"\"\n", "\n", " html += \"\"\n", " \n", " for idx, (desc, func) in enumerate(solvers):\n", " \n", " if bm_mat[idx] == m:\n", " html += \"\" % \\\n", " (desc, bm_mat[idx], bm_mat[idx]/m)\n", " else:\n", " html += \"\" % \\\n", " (desc, bm_mat[idx], bm_mat[idx]/m)\n", " \n", " \n", " html += \"
SolverElapsed timeRatio
%s%.8f%.2f
%s%.8f%.2f
\"\n", "\n", " return HTML(html)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "def benchmark_steadystate_solvers(args, solvers, problem_func):\n", "\n", " bm_mat = zeros(len(solvers))\n", "\n", " H, c_ops = problem_func(args)\n", " \n", " for sol_idx, solver in enumerate(solvers):\n", " solver_name, solver_func = solver\n", " try:\n", " t1 = time.time()\n", " for r in range(reps):\n", " rhoss = solver_func(H, c_ops)\n", " t2 = time.time()\n", " bm_mat[sol_idx] = (t2 - t1)/reps\n", " \n", " except:\n", " bm_mat[sol_idx] = nan\n", " \n", " return bm_mat" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 }, { "cell_type": "code", "collapsed": false, "input": [ "solvers = [[\"power use_umfpack=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='power', use_umfpack=True)],\n", " [\"power use_umfpack=False\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='power', use_umfpack=False)],\n", " [\"direct sparse use_umfpack=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, use_umfpack=True, sparse=True)],\n", " [\"direct sparse use_umfpack=False\", \n", " lambda H, c_ops: steadystate(H, c_ops, use_umfpack=False, sparse=True)],\n", " [\"iterative use_precond=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True)],\n", " [\"iterative use_precond=False\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False)],\n", " [\"iterative-bicg use_precond=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True)],\n", " [\"iterative-bicg use_precond=False\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False)],\n", " [\"direct dense use_umfpack=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, use_umfpack=True, sparse=False)],\n", " [\"direct dense use_umfpack=False\", \n", " lambda H, c_ops: steadystate(H, c_ops, use_umfpack=False, sparse=False)],\n", " [\"svd_dense\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='svd')],\n", " [\"lu\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='lu')],\n", " ]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "code", "collapsed": false, "input": [ "large_solvers = [\n", " [\"power use_umfpack=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='power', use_umfpack=True)],\n", " [\"power use_umfpack=False\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='power', use_umfpack=False)],\n", " [\"direct sparse use_umfpack=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, use_umfpack=True, sparse=True)],\n", " [\"direct sparse use_umfpack=False\", \n", " lambda H, c_ops: steadystate(H, c_ops, use_umfpack=False, sparse=True)],\n", " [\"iterative use_precond=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True)],\n", " [\"iterative-bicg use_precond=True\", \n", " lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True)],\n", " ]" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Test problem 1" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def bm_problem1(N):\n", "\n", " a = tensor(destroy(N), identity(2))\n", " b = tensor(identity(N), destroy(2))\n", "\n", " H = a.dag() * a + b.dag() * b + 0.25 * (a + a.dag()) * (b + b.dag())\n", "\n", " c_ops = [sqrt(0.1) * a, sqrt(0.075) * a.dag(), sqrt(0.1) * b]\n", " \n", " return H, c_ops" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "code", "collapsed": false, "input": [ "bm_mat = benchmark_steadystate_solvers(10, solvers, bm_problem1)\n", "show_bm_mat(bm_mat, solvers)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
SolverElapsed timeRatio
power use_umfpack=True0.029301171.86
power use_umfpack=False0.027691841.76
direct sparse use_umfpack=True0.016577721.05
direct sparse use_umfpack=False0.015777111.00
iterative use_precond=True0.018871311.20
iterative use_precond=False0.4271519227.07
iterative-bicg use_precond=True0.018245221.16
iterative-bicg use_precond=False0.059676653.78
direct dense use_umfpack=True0.036807062.33
direct dense use_umfpack=False0.035973552.28
svd_dense0.4354958527.60
lu0.015813591.00
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 11, "text": [ "" ] } ], "prompt_number": 11 }, { "cell_type": "code", "collapsed": false, "input": [ "bm_mat = benchmark_steadystate_solvers(50, large_solvers, bm_problem1)\n", "show_bm_mat(bm_mat, large_solvers)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
SolverElapsed timeRatio
power use_umfpack=True0.308600662.16
power use_umfpack=False0.309067492.16
direct sparse use_umfpack=True0.262234451.83
direct sparse use_umfpack=False0.304351572.13
iterative use_precond=True0.143191581.00
iterative-bicg use_precond=True0.144336701.01
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 12, "text": [ "" ] } ], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Test problem 2: high temperature harmonic oscillator" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def bm_problem2(N):\n", " \n", " a = destroy(N)\n", " H = a.dag() * a\n", " nth = N / 4\n", " gamma = 0.05\n", " c_ops = [sqrt(gamma * (nth + 1)) * a, sqrt(gamma * nth) * a.dag()]\n", "\n", " return H, c_ops" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 13 }, { "cell_type": "code", "collapsed": false, "input": [ "bm_mat = benchmark_steadystate_solvers(50, solvers, bm_problem2)\n", "show_bm_mat(bm_mat, solvers)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
SolverElapsed timeRatio
power use_umfpack=True0.030913832.33
power use_umfpack=False0.027934312.11
direct sparse use_umfpack=True0.015746591.19
direct sparse use_umfpack=False0.013266801.00
iterative use_precond=True0.018147471.37
iterative use_precond=False0.087285526.58
iterative-bicg use_precond=True0.017312291.30
iterative-bicg use_precond=False0.031406402.37
direct dense use_umfpack=True3.13456941236.27
direct dense use_umfpack=False3.13574004236.36
svd_dense71.838504085414.91
lu0.013910771.05
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 14, "text": [ "" ] } ], "prompt_number": 14 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Test problem 3: Coupled oscillators" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def bm_problem3(N):\n", " \n", " a = tensor(destroy(N), identity(N))\n", " b = tensor(identity(N), destroy(N))\n", " \n", " H = a.dag() * a + 0.25 * b.dag() * b + 0.05 * a.dag() * a * (b + b.dag()) + 0.1 * (a + a.dag()) \n", "\n", " c_ops = [sqrt(0.05) * a, sqrt(0.015) * a.dag(), sqrt(0.1) * b, sqrt(0.075) * b.dag()]\n", "\n", " return H, c_ops" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, "input": [ "bm_mat = benchmark_steadystate_solvers(10, large_solvers, bm_problem3)\n", "show_bm_mat(bm_mat, large_solvers)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
SolverElapsed timeRatio
power use_umfpack=True5.094226367.99
power use_umfpack=False35.7943534956.16
direct sparse use_umfpack=True5.805845029.11
direct sparse use_umfpack=False30.2062852447.40
iterative use_precond=True0.716620211.12
iterative-bicg use_precond=True0.637322901.00
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 16, "text": [ "" ] } ], "prompt_number": 16 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Test problem 4: a two level system " ] }, { "cell_type": "code", "collapsed": false, "input": [ "def bm_problem4(args=None):\n", "\n", " sz = sigmaz() \n", " sx = sigmax()\n", " \n", " H = sz\n", " c_ops = [sqrt(0.05) * sx]\n", "\n", " return H, c_ops" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 17 }, { "cell_type": "code", "collapsed": false, "input": [ "bm_mat = benchmark_steadystate_solvers(None, solvers, bm_problem4)\n", "show_bm_mat(bm_mat, solvers)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
SolverElapsed timeRatio
power use_umfpack=True0.018521793.05
power use_umfpack=False0.017549752.89
direct sparse use_umfpack=True0.007075311.17
direct sparse use_umfpack=False0.006877901.13
iterative use_precond=True0.007096771.17
iterative use_precond=False0.007111071.17
iterative-bicg use_precond=True0.016734362.76
iterative-bicg use_precond=False0.006780621.12
direct dense use_umfpack=True0.006115441.01
direct dense use_umfpack=False0.006063941.00
svd_dense0.007694961.27
lu0.006518361.07
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 18, "text": [ "" ] } ], "prompt_number": 18 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Test problem 5: spin chain" ] }, { "cell_type": "code", "collapsed": false, "input": [ "def bm_problem5(N=1):\n", "\n", " H = 0\n", " for m in range(N):\n", " H += tensor([sigmaz() if n == m else identity(2) for n in range(N)])\n", "\n", " for m in range(N-1):\n", " H += tensor([sigmax() if n in [m,m+1] else identity(2) for n in range(N)]) \n", " \n", " c_ops = [sqrt(0.05) * tensor([sigmam() if n == m else identity(2) for n in range(N)])\n", " for m in range(N)]\n", " \n", " return H, c_ops" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 19 }, { "cell_type": "code", "collapsed": false, "input": [ "bm_mat = benchmark_steadystate_solvers(2, solvers, bm_problem5)\n", "show_bm_mat(bm_mat, solvers)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
SolverElapsed timeRatio
power use_umfpack=True0.021501302.35
power use_umfpack=False0.020940782.28
direct sparse use_umfpack=True0.010230061.12
direct sparse use_umfpack=False0.009975911.09
iterative use_precond=True0.010204791.11
iterative use_precond=False0.010334491.13
iterative-bicg use_precond=True0.009714841.06
iterative-bicg use_precond=False0.010103461.10
direct dense use_umfpack=True0.009189841.00
direct dense use_umfpack=False0.009165531.00
svd_dense0.010971311.20
lu0.009572031.04
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 20, "text": [ "" ] } ], "prompt_number": 20 }, { "cell_type": "code", "collapsed": false, "input": [ "bm_mat = benchmark_steadystate_solvers(4, solvers, bm_problem5)\n", "show_bm_mat(bm_mat, solvers)" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
SolverElapsed timeRatio
power use_umfpack=True0.031192781.59
power use_umfpack=False0.030974861.58
direct sparse use_umfpack=True0.019570111.00
direct sparse use_umfpack=False0.019614461.00
iterative use_precond=True0.025650021.31
iterative use_precond=False0.2790315214.26
iterative-bicg use_precond=True0.024167301.23
iterative-bicg use_precond=False0.066788433.41
direct dense use_umfpack=True0.024587151.26
direct dense use_umfpack=False0.024188281.24
svd_dense0.165172348.44
lu0.019600631.00
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 21, "text": [ "" ] } ], "prompt_number": 21 }, { "cell_type": "code", "collapsed": false, "input": [ "bm_mat = benchmark_steadystate_solvers(6, large_solvers, bm_problem5)\n", "show_bm_mat(bm_mat, large_solvers)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stderr", "text": [ "/usr/local/lib/python3.3/dist-packages/qutip/steadystate.py:314: UserWarning: Preconditioning failed. Continuing without.\n", " UserWarning)\n" ] }, { "html": [ "
SolverElapsed timeRatio
power use_umfpack=True1.099171881.08
power use_umfpack=False9.799975409.60
direct sparse use_umfpack=True1.021354201.00
direct sparse use_umfpack=False9.188909779.00
iterative use_precond=True13.4802668113.20
iterative-bicg use_precond=True2.671352622.62
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 22, "text": [ "" ] } ], "prompt_number": 22 }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Software versions" ] }, { "cell_type": "code", "collapsed": false, "input": [ "from qutip.ipynbtools import version_table\n", "\n", "version_table()" ], "language": "python", "metadata": {}, "outputs": [ { "html": [ "
SoftwareVersion
OSposix [linux]
IPython1.0.dev
Numpy1.8.0.dev-c5694c5
QuTiP2.3.0.dev-169f358
Cython0.19.1
SciPy0.13.0.dev-38faf7c
Python3.3.1 (default, Apr 17 2013, 22:30:32) \n", "[GCC 4.7.3]
matplotlib1.4.x
Fri Jul 19 14:08:27 2013 JST
" ], "metadata": {}, "output_type": "pyout", "prompt_number": 23, "text": [ "" ] } ], "prompt_number": 23 } ], "metadata": {} } ] }