{
 "metadata": {
  "name": "",
  "signature": "sha256:d11bd246f68abb950adb5879d727cb8eb58128a5c7c950a747324e287d3ac190"
 },
 "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": [
        "Populating the interactive namespace from numpy and matplotlib\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 = 10"
     ],
     "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 = \"<table>\"\n",
      "\n",
      "    html += \"<tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr>\"\n",
      "    \n",
      "    for idx, (desc, func) in enumerate(solvers):\n",
      "    \n",
      "        if bm_mat[idx] == m:\n",
      "            html += \"<tr><td><b>%s</b></td><td><b>%.8f</b></td><td><b>%.2f</b></td></tr>\" % \\\n",
      "                (desc, bm_mat[idx], bm_mat[idx]/m)\n",
      "        else:\n",
      "            html += \"<tr><td>%s</td><td>%.8f</td><td>%.2f</td></tr>\" % \\\n",
      "                (desc, bm_mat[idx], bm_mat[idx]/m)\n",
      "            \n",
      "    \n",
      "    html += \"</table>\"\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 Exception as e:\n",
      "            bm_mat[sol_idx] = nan\n",
      "            print(\"Failure in %s: %s\" % (solver_name, str(e)))\n",
      "\n",
      "    return bm_mat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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",
      "\n",
      "    [\"direct sparse use_umfpack=True\",  \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='direct', use_umfpack=True, sparse=True)],\n",
      "    [\"direct sparse use_umfpack=False\", \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='direct', use_umfpack=False, sparse=True)],\n",
      "    [\"direct dense use_umfpack=True\",   \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='direct', use_umfpack=True, sparse=False)],\n",
      "    [\"direct dense use_umfpack=False\",  \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='direct', use_umfpack=False, sparse=False)],\n",
      "\n",
      "    [\"iterative use_precond=True, use_rcm=True, sym=True\",      \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=True, sym=True)],\n",
      "    [\"iterative use_precond=True, use_rcm=False, sym=True\",   \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=False, sym=True)],\n",
      "    [\"iterative use_precond=False, use_rcm=True, sym=True\",     \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False, use_rcm=True, sym=True)],\n",
      "    [\"iterative use_precond=False, use_rcm=False, sym=True\",     \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False, use_rcm=False, sym=True)],\n",
      "    [\"iterative use_precond=True, use_rcm=True, sym=False\",      \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=True, sym=False)],\n",
      "    [\"iterative use_precond=True, use_rcm=False, sym=False\",   \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=False, sym=False)],\n",
      "    [\"iterative use_precond=False, use_rcm=True, sym=False\",     \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False, use_rcm=True, sym=False)],\n",
      "    [\"iterative use_precond=False, use_rcm=False, sym=False\",     \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=False, use_rcm=False, sym=False)],\n",
      "\n",
      "    [\"iterative-bicg use_precond=True, use_rcm=True, sym=True\",  \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=True, sym=True)],\n",
      "    [\"iterative-bicg use_precond=False, use_rcm=False, sym=True\", \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False, use_rcm=False, sym=True)],\n",
      "    [\"iterative-bicg use_precond=True, use_rcm=True, sym=True\",  \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=True, sym=True)],\n",
      "    [\"iterative-bicg use_precond=False, use_rcm=False, sym=True\", \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False, use_rcm=False, sym=True)],\n",
      "    [\"iterative-bicg use_precond=True, use_rcm=True, sym=False\",      \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=True, sym=False)],\n",
      "    [\"iterative-bicg use_precond=True, use_rcm=False, sym=False\",   \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=False, sym=False)],\n",
      "    [\"iterative-bicg use_precond=False, use_rcm=True, sym=False\",     \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False, use_rcm=True, sym=False)],\n",
      "    [\"iterative-bicg use_precond=False, use_rcm=False, sym=False\",     \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=False, use_rcm=False, sym=False)],\n",
      "\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",
      "\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",
      "\n",
      "    [\"iterative use_precond=True, use_rcm=True, sym=False\",      \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=True, sym=False)],\n",
      "    [\"iterative use_precond=True, use_rcm=False, sym=False\",      \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative', use_precond=True, use_rcm=False, sym=False)],\n",
      "\n",
      "    [\"iterative-bicg use_precond=True, use_rcm=True, sym=False\",      \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=True, sym=False)],\n",
      "    [\"iterative-bicg use_precond=True, use_rcm=False, sym=False\",      \n",
      "     lambda H, c_ops: steadystate(H, c_ops, method='iterative-bicg', use_precond=True, use_rcm=False, sym=False)],\n",
      "\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": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Failure in iterative use_precond=False, use_rcm=True, sym=True: Steadystate solver did not reach tolerance after 1000 steps.\n",
        "Failure in iterative use_precond=False, use_rcm=False, sym=True: Steadystate solver did not reach tolerance after 1000 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Failure in iterative use_precond=False, use_rcm=True, sym=False: Steadystate solver did not reach tolerance after 1000 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Failure in iterative use_precond=False, use_rcm=False, sym=False: Steadystate solver did not reach tolerance after 1000 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "/home/rob/py-envs/py3-stable/lib/python3.3/site-packages/qutip/steadystate.py:153: UserWarning: The use of use_umfpack is deprecated.\n",
        "  warnings.warn(\"The use of use_umfpack is deprecated.\")\n",
        "-c:3: RuntimeWarning: invalid value encountered in greater\n"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "html": [
        "<table><tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr><tr><td>power use_umfpack=True</td><td>0.02337635</td><td>1.47</td></tr><tr><td>power use_umfpack=False</td><td>0.02320838</td><td>1.46</td></tr><tr><td>direct sparse use_umfpack=True</td><td>0.01801844</td><td>1.14</td></tr><tr><td>direct sparse use_umfpack=False</td><td>0.01740201</td><td>1.10</td></tr><tr><td>direct dense use_umfpack=True</td><td>0.04866381</td><td>3.07</td></tr><tr><td>direct dense use_umfpack=False</td><td>0.04866440</td><td>3.07</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=True</td><td>0.02051804</td><td>1.29</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=True</td><td>0.01893902</td><td>1.19</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=True</td><td>nan</td><td>nan</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=True</td><td>nan</td><td>nan</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=False</td><td>0.01817827</td><td>1.15</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=False</td><td>0.01661460</td><td>1.05</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=False</td><td>nan</td><td>nan</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=False</td><td>nan</td><td>nan</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.01868496</td><td>1.18</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.06821716</td><td>4.30</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.01880364</td><td>1.18</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.06844871</td><td>4.31</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=False</td><td>0.01848016</td><td>1.16</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=False, sym=False</td><td>0.01699817</td><td>1.07</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=True, sym=False</td><td>0.06557376</td><td>4.13</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=False</td><td>0.06788254</td><td>4.28</td></tr><tr><td><b>lu</b></td><td><b>0.01587019</b></td><td><b>1.00</b></td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "<IPython.core.display.HTML at 0x7f6597842d10>"
       ]
      }
     ],
     "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": [
        "<table><tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr><tr><td>power use_umfpack=True</td><td>0.27807622</td><td>1.19</td></tr><tr><td>power use_umfpack=False</td><td>0.28087587</td><td>1.20</td></tr><tr><td><b>direct sparse use_umfpack=True</b></td><td><b>0.23442681</b></td><td><b>1.00</b></td></tr><tr><td>direct sparse use_umfpack=False</td><td>0.28316166</td><td>1.21</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=False</td><td>0.27145686</td><td>1.16</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=False</td><td>0.27361963</td><td>1.17</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=False</td><td>0.27396398</td><td>1.17</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=False, sym=False</td><td>0.29887297</td><td>1.27</td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 12,
       "text": [
        "<IPython.core.display.HTML at 0x7f6597816dd0>"
       ]
      }
     ],
     "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": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Failure in iterative-bicg use_precond=True, use_rcm=False, sym=False: Steadystate solver failed with fatal error: -10.\n"
       ]
      },
      {
       "html": [
        "<table><tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr><tr><td>power use_umfpack=True</td><td>0.02198510</td><td>1.60</td></tr><tr><td>power use_umfpack=False</td><td>0.02195816</td><td>1.60</td></tr><tr><td>direct sparse use_umfpack=True</td><td>0.01739521</td><td>1.27</td></tr><tr><td>direct sparse use_umfpack=False</td><td>0.01502705</td><td>1.09</td></tr><tr><td>direct dense use_umfpack=True</td><td>5.38591220</td><td>391.95</td></tr><tr><td>direct dense use_umfpack=False</td><td>5.56329439</td><td>404.86</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=True</td><td>0.01940596</td><td>1.41</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=True</td><td>0.01877887</td><td>1.37</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=True</td><td>0.85532901</td><td>62.25</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=True</td><td>0.90501690</td><td>65.86</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=False</td><td>0.01919260</td><td>1.40</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=False</td><td>0.01824219</td><td>1.33</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=False</td><td>0.82320507</td><td>59.91</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=False</td><td>0.89947019</td><td>65.46</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.01872413</td><td>1.36</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.03652155</td><td>2.66</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.01884017</td><td>1.37</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.03594759</td><td>2.62</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=False</td><td>0.01916583</td><td>1.39</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=False, sym=False</td><td>nan</td><td>nan</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=True, sym=False</td><td>0.03907797</td><td>2.84</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=False</td><td>0.03581498</td><td>2.61</td></tr><tr><td><b>lu</b></td><td><b>0.01374128</b></td><td><b>1.00</b></td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 14,
       "text": [
        "<IPython.core.display.HTML at 0x7f6597854c50>"
       ]
      }
     ],
     "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": [
        "<table><tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr><tr><td>power use_umfpack=True</td><td>33.75028317</td><td>10.41</td></tr><tr><td>power use_umfpack=False</td><td>34.65209894</td><td>10.68</td></tr><tr><td>direct sparse use_umfpack=True</td><td>8.38044708</td><td>2.58</td></tr><tr><td>direct sparse use_umfpack=False</td><td>27.64443343</td><td>8.52</td></tr><tr><td><b>iterative use_precond=True, use_rcm=True, sym=False</b></td><td><b>3.24328792</b></td><td><b>1.00</b></td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=False</td><td>3.47413125</td><td>1.07</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=False</td><td>4.31313200</td><td>1.33</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=False, sym=False</td><td>3.44422872</td><td>1.06</td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 16,
       "text": [
        "<IPython.core.display.HTML at 0x7f659780f390>"
       ]
      }
     ],
     "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": [
        "<table><tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr><tr><td>power use_umfpack=True</td><td>0.01118724</td><td>1.95</td></tr><tr><td>power use_umfpack=False</td><td>0.01090584</td><td>1.90</td></tr><tr><td>direct sparse use_umfpack=True</td><td>0.00697067</td><td>1.21</td></tr><tr><td>direct sparse use_umfpack=False</td><td>0.00657489</td><td>1.14</td></tr><tr><td>direct dense use_umfpack=True</td><td>0.00577335</td><td>1.00</td></tr><tr><td><b>direct dense use_umfpack=False</b></td><td><b>0.00575178</b></td><td><b>1.00</b></td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=True</td><td>0.00739422</td><td>1.29</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=True</td><td>0.00644894</td><td>1.12</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=True</td><td>0.00735872</td><td>1.28</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=True</td><td>0.00638795</td><td>1.11</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=False</td><td>0.00745142</td><td>1.30</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=False</td><td>0.00646114</td><td>1.12</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=False</td><td>0.00748370</td><td>1.30</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=False</td><td>0.00631652</td><td>1.10</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.00740299</td><td>1.29</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.00633924</td><td>1.10</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.00731163</td><td>1.27</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.00634141</td><td>1.10</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=False</td><td>0.00732563</td><td>1.27</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=False, sym=False</td><td>0.00647802</td><td>1.13</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=True, sym=False</td><td>0.00735517</td><td>1.28</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=False</td><td>0.00637131</td><td>1.11</td></tr><tr><td>lu</td><td>0.00636835</td><td>1.11</td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "<IPython.core.display.HTML at 0x7f659780fad0>"
       ]
      }
     ],
     "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": [
        "<table><tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr><tr><td>power use_umfpack=True</td><td>0.01427205</td><td>1.59</td></tr><tr><td>power use_umfpack=False</td><td>0.01437883</td><td>1.60</td></tr><tr><td>direct sparse use_umfpack=True</td><td>0.01016231</td><td>1.13</td></tr><tr><td>direct sparse use_umfpack=False</td><td>0.00989394</td><td>1.10</td></tr><tr><td>direct dense use_umfpack=True</td><td>0.00903437</td><td>1.00</td></tr><tr><td><b>direct dense use_umfpack=False</b></td><td><b>0.00900269</b></td><td><b>1.00</b></td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=True</td><td>0.01079981</td><td>1.20</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=True</td><td>0.00970068</td><td>1.08</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=True</td><td>0.01065261</td><td>1.18</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=True</td><td>0.00960047</td><td>1.07</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=False</td><td>0.01077218</td><td>1.20</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=False</td><td>0.00982673</td><td>1.09</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=False</td><td>0.01066728</td><td>1.18</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=False</td><td>0.00975940</td><td>1.08</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.01067920</td><td>1.19</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.00998452</td><td>1.11</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.01066823</td><td>1.19</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.01001103</td><td>1.11</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=False</td><td>0.01083741</td><td>1.20</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=False, sym=False</td><td>0.00956621</td><td>1.06</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=True, sym=False</td><td>0.01108155</td><td>1.23</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=False</td><td>0.00998340</td><td>1.11</td></tr><tr><td>lu</td><td>0.00949795</td><td>1.06</td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 20,
       "text": [
        "<IPython.core.display.HTML at 0x7f659780fdd0>"
       ]
      }
     ],
     "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": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Failure in iterative use_precond=False, use_rcm=True, sym=True: Steadystate solver did not reach tolerance after 1000 steps.\n",
        "Failure in iterative use_precond=False, use_rcm=False, sym=True: Steadystate solver did not reach tolerance after 1000 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Failure in iterative use_precond=False, use_rcm=True, sym=False: Steadystate solver did not reach tolerance after 1000 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Failure in iterative use_precond=False, use_rcm=False, sym=False: Steadystate solver did not reach tolerance after 1000 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "html": [
        "<table><tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr><tr><td>power use_umfpack=True</td><td>0.02471032</td><td>1.25</td></tr><tr><td>power use_umfpack=False</td><td>0.02463510</td><td>1.25</td></tr><tr><td>direct sparse use_umfpack=True</td><td>0.01991639</td><td>1.01</td></tr><tr><td>direct sparse use_umfpack=False</td><td>0.01971462</td><td>1.00</td></tr><tr><td>direct dense use_umfpack=True</td><td>0.02616692</td><td>1.33</td></tr><tr><td>direct dense use_umfpack=False</td><td>0.02603197</td><td>1.32</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=True</td><td>0.02279491</td><td>1.16</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=True</td><td>0.02091999</td><td>1.06</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=True</td><td>nan</td><td>nan</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=True</td><td>nan</td><td>nan</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=False</td><td>0.02344491</td><td>1.19</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=False</td><td>0.02175648</td><td>1.10</td></tr><tr><td>iterative use_precond=False, use_rcm=True, sym=False</td><td>nan</td><td>nan</td></tr><tr><td>iterative use_precond=False, use_rcm=False, sym=False</td><td>nan</td><td>nan</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.02228684</td><td>1.13</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.07718029</td><td>3.92</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=True</td><td>0.02243965</td><td>1.14</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=True</td><td>0.07755005</td><td>3.94</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=False</td><td>0.02243128</td><td>1.14</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=False, sym=False</td><td>0.02093019</td><td>1.06</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=True, sym=False</td><td>0.07641232</td><td>3.88</td></tr><tr><td>iterative-bicg use_precond=False, use_rcm=False, sym=False</td><td>0.07933478</td><td>4.03</td></tr><tr><td><b>lu</b></td><td><b>0.01969354</b></td><td><b>1.00</b></td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 21,
       "text": [
        "<IPython.core.display.HTML at 0x7f6597808490>"
       ]
      }
     ],
     "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": "stdout",
       "text": [
        "Failure in iterative use_precond=True, use_rcm=True, sym=False: Steadystate solver did not reach tolerance after 1000 steps.\n",
        "Failure in iterative use_precond=True, use_rcm=False, sym=False: Steadystate solver did not reach tolerance after 1000 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Failure in iterative-bicg use_precond=True, use_rcm=True, sym=False: Steadystate solver did not reach tolerance after 40960 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "Failure in iterative-bicg use_precond=True, use_rcm=False, sym=False: Steadystate solver did not reach tolerance after 40960 steps."
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      },
      {
       "html": [
        "<table><tr><td><b>Solver</b></td><td><b>Elapsed time</b></td><td><b>Ratio</b></td></tr><tr><td>power use_umfpack=True</td><td>10.34370623</td><td>8.90</td></tr><tr><td>power use_umfpack=False</td><td>10.37572424</td><td>8.93</td></tr><tr><td><b>direct sparse use_umfpack=True</b></td><td><b>1.16162479</b></td><td><b>1.00</b></td></tr><tr><td>direct sparse use_umfpack=False</td><td>9.75404005</td><td>8.40</td></tr><tr><td>iterative use_precond=True, use_rcm=True, sym=False</td><td>nan</td><td>nan</td></tr><tr><td>iterative use_precond=True, use_rcm=False, sym=False</td><td>nan</td><td>nan</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=True, sym=False</td><td>nan</td><td>nan</td></tr><tr><td>iterative-bicg use_precond=True, use_rcm=False, sym=False</td><td>nan</td><td>nan</td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "<IPython.core.display.HTML at 0x7f659779f150>"
       ]
      }
     ],
     "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": [
        "<table><tr><th>Software</th><th>Version</th></tr><tr><td>Python</td><td>3.3.2+ (default, Oct  9 2013, 14:50:09) \n",
        "[GCC 4.8.1]</td></tr><tr><td>matplotlib</td><td>1.3.1</td></tr><tr><td>QuTiP</td><td>3.0.0.dev-0f20b2a</td></tr><tr><td>Cython</td><td>0.20.1</td></tr><tr><td>SciPy</td><td>0.13.3</td></tr><tr><td>Numpy</td><td>1.8.0</td></tr><tr><td>IPython</td><td>2.0.0</td></tr><tr><td>OS</td><td>posix [linux]</td></tr><tr><td colspan='2'>Thu Apr 03 17:04:46 2014 JST</td></tr></table>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 23,
       "text": [
        "<IPython.core.display.HTML at 0x7f659780f190>"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    }
   ],
   "metadata": {}
  }
 ]
}