{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### General information\n",
    "\n",
    "http://slepc.upv.es/slepc4py-current/docs/usrman/install.html#using-pip-or-easy-install\n",
    "\n",
    "\n",
    "http://slepc.upv.es/handson/handson1.html\n",
    "\n",
    "\n",
    "\n",
    "**How to calculate the smallest eigenvalues (instead of the biggest)**:\n",
    "\n",
    "http://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetWhichEigenpairs.html#EPSSetWhichEigenpairs\n",
    "\n",
    "`opts.setValue(\"eps_smallest_magnitude\", 1)`\n",
    "\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, slepc4py\n",
    "slepc4py.init(sys.argv)\n",
    "\n",
    "# Next, we have to import the relevant modules. Normally, both PETSc and SLEPc modules have to be imported in all slepc4py programs. It may be useful to import NumPy as well:\n",
    "\n",
    "from petsc4py import PETSc\n",
    "from slepc4py import SLEPc\n",
    "import numpy\n",
    "np.set_printoptions(linewidth=270, precision=4)\n",
    "\n",
    "# At this point, we can use any petsc4py and slepc4py operations. For instance, the following lines allow the user to specify an integer command-line argument n with a default value of 30 (see the next section for example usage of command-line options):\n",
    "\n",
    "opts = PETSc.Options()\n",
    "n = opts.getInt('n', 3000)\n",
    "opts.setValue(\"eps_max_it\", 100000)\n",
    "opts.setValue(\"eps_smallest_magnitude\", 1)\n",
    "\n",
    "\n",
    "opts.setValue(\"eps_nev\", 10)\n",
    "\n",
    "\n",
    "# It is necessary to build a matrix to define an eigenproblem (or two in the case of generalized eigenproblems). The following fragment of code creates the matrix object and then fills the non-zero elements one by one. The matrix of this particular example is tridiagonal, with value 2 in the diagonal, and -1 in off-diagonal positions. See petsc4py documentation for details about matrix objects:\n",
    "\n",
    "A = PETSc.Mat().create()\n",
    "A.setSizes([n, n])\n",
    "A.setFromOptions()\n",
    "A.setUp()\n",
    "\n",
    "rstart, rend = A.getOwnershipRange()\n",
    "\n",
    "# first row\n",
    "if rstart == 0:\n",
    "    A[0, :2] = [2, -1]\n",
    "    rstart += 1\n",
    "# last row\n",
    "if rend == n:\n",
    "    A[n-1, -2:] = [-1, 2]\n",
    "    rend -= 1\n",
    "# other rows\n",
    "for i in range(rstart, rend):\n",
    "    A[i, i-1:i+2] = [-1, 2, -1]\n",
    "\n",
    "A.assemble()\n",
    "\n",
    "# The solver object is created in a similar way as other objects in petsc4py:\n",
    "\n",
    "E = SLEPc.EPS(); E.create()\n",
    "\n",
    "# Once the object is created, the eigenvalue problem must be specified. At least one matrix must be provided. The problem type must be indicated as well, in this case it is HEP (Hermitian eigenvalue problem). Apart from these, other settings could be provided here (for instance, the tolerance for the computation). After all options have been set, the user should call the setFromOptions() operation, so that any options specified at run time in the command line are passed to the solver object:\n",
    "\n",
    "E.setOperators(A)\n",
    "E.setProblemType(SLEPc.EPS.ProblemType.HEP)\n",
    "E.setFromOptions()\n",
    "\n",
    "# After that, the solve() method will run the selected eigensolver, keeping the solution stored internally:\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.81 s, sys: 2.7 s, total: 6.51 s\n",
      "Wall time: 1.73 s\n"
     ]
    }
   ],
   "source": [
    "%time E.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "******************************\n",
      "*** SLEPc Solution Results ***\n",
      "******************************\n",
      "\n",
      "Number of iterations of the method: 2206\n",
      "Solution method: krylovschur\n",
      "Number of requested eigenvalues: 10\n",
      "Stopping condition: tol=1e-08, maxit=100000\n"
     ]
    }
   ],
   "source": [
    "# Print = PETSc.Sys.Print\n",
    "Print = print\n",
    "\n",
    "Print()\n",
    "Print(\"******************************\")\n",
    "Print(\"*** SLEPc Solution Results ***\")\n",
    "Print(\"******************************\")\n",
    "Print()\n",
    "\n",
    "its = E.getIterationNumber()\n",
    "Print(\"Number of iterations of the method: %d\" % its)\n",
    "\n",
    "eps_type = E.getType()\n",
    "Print(\"Solution method: %s\" % eps_type)\n",
    "\n",
    "nev, ncv, mpd = E.getDimensions()\n",
    "Print(\"Number of requested eigenvalues: %d\" % nev)\n",
    "\n",
    "tol, maxit = E.getTolerances()\n",
    "Print(\"Stopping condition: tol=%.4g, maxit=%d\" % (tol, maxit))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of converged eigenpairs 11\n",
      "\n",
      "        k          ||Ax-kx||/||kx|| \n",
      "----------------- ------------------\n",
      "     0.000001       2.44656e-08\n",
      "     0.000004       7.39994e-09\n",
      "     0.000010       2.55254e-09\n",
      "     0.000018       4.22676e-09\n",
      "     0.000027         1.597e-09\n",
      "     0.000039       3.96912e-09\n",
      "     0.000054       1.32783e-09\n",
      "     0.000070       1.59263e-09\n",
      "     0.000089       1.45552e-09\n",
      "     0.000110        4.0202e-09\n",
      "     0.000133       6.58014e-09\n",
      "\n"
     ]
    }
   ],
   "source": [
    "nconv = E.getConverged()\n",
    "Print(\"Number of converged eigenpairs %d\" % nconv)\n",
    "kk = []\n",
    "if nconv > 0:\n",
    "    # Create the results vectors\n",
    "    vr, wr = A.getVecs()\n",
    "    vi, wi = A.getVecs()\n",
    "    #\n",
    "    Print()\n",
    "    Print(\"        k          ||Ax-kx||/||kx|| \")\n",
    "    Print(\"----------------- ------------------\")\n",
    "    for i in range(nconv):\n",
    "        k = E.getEigenpair(i, vr, vi)\n",
    "        kk.append(k)\n",
    "        error = E.computeError(i)\n",
    "        if k.imag != 0.0:\n",
    "            Print(\" %9f%+9f j %12g\" % (k.real, k.imag, error))\n",
    "        else:\n",
    "            Print(\" %12f      %12g\" % (k.real, error))\n",
    "    \n",
    "    kk = np.array(kk)\n",
    "    Print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sc\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert to numpy array\n",
    "Anp = A[:,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 51.9 s, sys: 20 s, total: 1min 11s\n",
      "Wall time: 20 s\n"
     ]
    }
   ],
   "source": [
    "%time r0 = np.linalg.eigvals(Anp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f6d3e8bfac8>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl8VNXdx/HPLztLCAQiS1jC5oIgIJFVrLgiWreiYqtFq1KsWq36+KBVqz5u1dZatdVStVpbd62ioChbFRUwrIKA7BDWsCVsgSzn+WOGkGWSTMLs+b5fr7y499wzM7/LTH45c+6555hzDhERiS1x4Q5AREQCT8ldRCQGKbmLiMQgJXcRkRik5C4iEoOU3EVEYpCSu4hIDFJyFxGJQUruIiIxKCFcL9yqVSuXlZUVrpcXEYlKc+fO3e6cy6itXtiSe1ZWFjk5OeF6eRGRqGRm6/ypp24ZEZEYpOQuIhKDlNxFRGKQkruISAxSchcRiUF+J3czizez+Wb2sY9jyWb2lpmtNLPZZpYVyCBFRKRu6tJyvxVYWs2x64BdzrluwJ+A3x9tYCIiUn9+jXM3s/bA+cAjwO0+qlwEPODdfhd4zszMBWkNv0++28zSzQWM/3I1vTLTOP24Y5i1egcJccaALi0pOFBE2+aN6JjemL2FxTRNSSCrZWNaN0shzoxS50hOiGNLQSHpTZJITogPRpgiImHj701MTwN3AanVHM8ENgA454rNLB9oCWwvX8nMxgBjADp27FifeHHOcfvbCzlQVALAt2t38e3aXWXHpy/Pq/Nzrn50BHFxVq94REQiUa3dMmZ2AbDNOTe3pmo+yqq02p1z451z2c657IyMWu+e9enLFdvLEnug/N/E75n03WYKA/y8IiLh4k/LfQhwoZmNAFKAZmb2L+fcVeXq5AIdgFwzSwDSgJ0BjxYodY42zVJ44MITGfsvz9+brhlNWJW3jy4ZTbh6YCf2HyqhbVoK3Y5pyqHiUt7JyeXBi04kJbFi98urX6/ldxOW8I+v1vKPr9YCEGfQqWUTkuLjuP2cY3EOlm4uoFmjRC7pm0l6k6RgnJaISEBZXbrFzex04E7n3AWVym8CejnnxprZKOBS59zlNT1Xdna2q8/cMs45zALXheKcY/aanTwwYQnLtuypse41g7P43Y97sGHnATq2bBywGERE/GVmc51z2bXVq/fEYWb2EJDjnJsAvAS8ZmYr8bTYR9X3ef143YA/38AuLfn7z7MZ+sT0Gut2SG9MzrpdXPbCNwBktWzM2h37GXfe8VzcJ5OlmwvISE2mVdNk2qSlBDROEZG6qFPLPZDq23IPFuccL81cQ2pKAie2S6N1sxTSGiWSlBDHrn2H6Pt/n3Pv+SfQMb0xY16r6fKDx5MjT+Ky7A4hiFxEGpKgt9xjjZlx/dAuPo8lJ3quOxeXOopL/ftj+PtPlyu5i0jYKLn7Id47TPKVr9aSmuL7v+zVX/Rna34hd723CICO6Y1CFp+ISGVK7n5IjPO03LcUFNKyaTNuGtaVId1aMbBzyyrj4y/Lbk/nuyexp7CYq1+azYFDJbwzdlDArxWIiNREyd0PcXHGv64bQGpKAie1T6sxUR8+tmLbXtZs30dxqaPfw1O49czujB6cFaKIRaSh0wXVIPj8+60YcLC4lJten1fh2E8HdOTRS3qFJzARiXr+XlBVcg8i5xx7DhYzc8V27v1gMTv3HSo79uszu3P72ceGMToRiUZK7hGmtNTx/H9X8eTk5RXKh3RryW9H9KBHu2ZhikxEoomSewRbsimf85+ZWaV88m2nkRBvdM1oGoaoRCQaaJx7BDuxXZrP8nOf/gKAlk2SuHvECTROimdEr7ahDE1EYoSW2YtAO/Yd4s53FvKrf8/jjrcXlpX/Z34uWeMmMmdNUOZkE5EYopZ7hHtvXi7vzcutUHb5375h+cPDtciIiFRLLfcotSW/MNwhiEgEU3KPUuZzfRQREQ8l9ygw7LgM3v7loAplxaWlYYpGRKKB+tzDJDHeKCqpOgz1cPmFvdvxh8t6k5Rw5O/vmsdG8OKXa3hk0lLy9hzkmakruHpQJ/p1Sg9l6CISBWpN7maWAnwBJHvrv+uc+12lOtcATwIbvUXPOedeDGyoseWxS0/iiU+X8bMBneiS0YRuxzSlReOkGhf5MDM6eGebvGL8LAAOlZTSr1M6O/YeJDkxnqbJ+nstIv613A8CZzjn9ppZIjDTzD5xzs2qVO8t59zNgQ8xNo3s156R/drX+XH7DnoW8U5JjKOwqJT0Jknc+8F3/GvWegDWPn5+QOMUkehUa3J3nltY93p3E70/4bmtVTi3Zxu2FBTyiyGdOenByWVJXUSkPL++w5tZPDAX6Ab8xTk320e1n5jZacAPwG+ccxsCF6Yc1jQ5gZuGdQMgNSWRnfsOcenJmXy/ybN+69SlW7nu1SPTOqglL9Iw+TVaxjlX4pzrA7QH+ptZz0pVPgKynHMnAVOAV309j5mNMbMcM8vJy8s7mrgFeOOGgUy/83SeurwPqSkJlJS6CokdYE9hUZiiE5FwqtNQSOfcbmAGMLxS+Q7n3EHv7t+BftU8frxzLts5l52RkVGPcKW849qk0rlVEwCKShxfr9pRpc7hPnoRaVhqTe5mlmFmzb3bjYCzgGWV6pSf3epCYGkgg5TaLdiwG4BbzujGykfO43/OPQ6AgY9N5dWv14YxMhEJB39a7m2B6Wa2CPgW+Nw597GZPWRmF3rr/NrMlpjZQuDXwDXBCVeqc+2QLAZ0TueOc44jIT6urEUP8LsJS8gaN5FtBZqyQKSh0HzuMaqopJQvV+Txi1eq/h+//6vB5O8v4t15ufzp8j4VbpQSkcim+dwbuMT4OM44vrXPY5f+9euy7YmLNvP9Q+fSOEkfBZFYoiabkH9AI2pEYo2Su1BSqnvSRGKNknsDdGK7Zgzq0rJsHppNuwvJ36/Wu0gsUUdrjDu+TSrLtuwhOSGOP4/qy/CebcqOfbRwE7e8MZ/L//YNP+7djmev7AuAcw4zzRcvEs2U3GPcx7ecSkK87y9oTVOOvP3Tl20ja9zEsv3nftqXC05qF/T4RCQ41C0T46pL7ABDu7XinbGDGNA5nb0Hiyscm7p0W7BDE5EgUnJvwBLi4zglK93naJlWTZPCEJGIBIqSu9C9dSoAD154Irec4Zlx8u9frmHsa3NZlLs7nKGJSD3pDlVhT2ERh4pLadk0mc35Bxj02LQKx+89/wR+PihLd7KKRAB/71DVb6uQmpJIy6bJAMT5GCXz8MSlHHvvJzz12fJQhyYi9aTkLhXUNALymWkreeLTZWzcfQCAgsIiLn/hG37YuidE0YmIvzQUUirw1XIv768zVvHXGatonBTP/kOeueIf/2QZL19zSijCExE/qeUuFdSW3A87nNgBktUXLxJx9FspFcTV48bU9i0aBT4QETkqSu5SQX2mHSgpDUIgInJU/FlmL8XM5pjZQu9qSw/6qJNsZm+Z2Uozm21mWcEIVoLPn5Z7UkIc7904iDvPORaA/YeKa3mEiISaPxdUDwJnOOf2mlkiMNPMPnHOzSpX5zpgl3Oum5mNAn4PXBGEeCXIqrvr4aZhXbnu1C40ToonJTEegH6d0nlv3kbe/HYD89bv4tNbTyOuPv06IhJwtSZ357nLaa93N9H7UzkHXAQ84N1+F3jOzMyF6w4pqbfyqblLqyb84fLe9GnfvNqkvWb7PgB+2LqXDxdupH/nluzad4j1O/czoldbn48RkeDzayikmcUDc4FuwF+cc7MrVckENgA454rNLB9oCWwPYKwSAqkpibz482x6ZqbRJi2l1vo/Obk9783LBeA3by2scGzNYyM0dbBImPh1QdU5V+Kc6wO0B/qbWc9KVXz9BldptZvZGDPLMbOcvLy8ukcrIXFWj9Z+JXaAP17em4cvrvxx8JiwcFMgwxKROqjTaBnn3G5gBjC80qFcoAOAmSUAacBOH48f75zLds5lZ2Rk1CtgiTyJ8b5b509OXs7s1TtCHI2IgH+jZTLMrLl3uxFwFrCsUrUJwGjv9khgmvrbG47yNz7dcfaxZdMF5+46wBXjZ1FS6vhwwcayUTVFJaVc8485SvwiQeRPn3tb4FVvv3sc8LZz7mMzewjIcc5NAF4CXjOzlXha7KOCFrFEnBPaNuP4Nqk8emkvTu7YgnN7tuGcP31RdrzrPZMAuLhPO54e1ZetBYXMWJ7Hotx85t13drjCFolp/oyWWQT09VF+f7ntQuCywIYm0aJnZhqf3nZa2X51Uxh8sGATqSmJvDZrXY31ROTo6Q5VCbiahrofTuwA2/ceDEE0Ig2TkrsEXLxuZBIJOyV3CTh1t4iEn5K7BJxyu0j4abEOCThf3TJDu7eiU8vG9GyXxrj3vysrd87pLlaRIFByl4Czcjcs//rM7tx+9rEVkvio/h25+/1FvDFnA1f8bRav/OIUGifpoygSSOqWkYBLTfEk6t+OOIHbz/ZMC1y5dX7Au5LTnLU76XH/ZN6cs56vV27n6Sk/8P2mgtAGLBKDLFw3kmZnZ7ucnJywvLYEX23dLSu27uHscjc6Vbb28fODEZZI1DOzuc657NrqqeUuQVFbP3r31qkMO07zC4kEi5K7hE1CvD5+IsGi3y4Jm9V5e6s9tqewKISRiMQeJXcJm1V5nlWcXr9+AKsfHVGhn/2+DxazafcBxr42l6xxE5m4aHO4whSJShp/JmHz4U1D2FpQyOBurcrK+nVqwdx1u8jddYDBj08rK3/hv6s4/yQt2yfiL7XcJWx6d2jOOSe2qVD28uhTAMhZt6tCue5zEqkbJXeJKIfHyFfWp0PzEEciEt2U3CWixFWauuDK/h0BmLliO4Mfm6ppgkX85M8yex3MbLqZLTWzJWZ2q486p5tZvpkt8P7c7+u5RPyVnBDH9w+dy2OX9qJNsxRWb9/HpvxCsh+ewiV//Src4YlEPH8uqBYDdzjn5plZKjDXzD53zn1fqd6XzrkLAh+iNDRTbv8R7ZqnlM03U3kesvnrd5M1biKz7zmTZ6etoElSApMWb2bybadpjhoRL3+W2dsMbPZu7zGzpUAmUDm5iwREt2OaVtiv3FVz2AXPziRvz5FumhVb99JbffMiQB373M0sC896qrN9HB5kZgvN7BMzO7Gax48xsxwzy8nLy6tzsNIwVbf4R/nEDpCSGB+KcESigt/J3cyaAu8BtznnKk/bNw/o5JzrDTwLfODrOZxz451z2c657IwMzSsi/vF32b5b35xP7q79QY5GJDr4ldzNLBFPYv+3c+79ysedcwXOub3e7UlAopm1qlxPpD78HeO+bMseHp20NLjBiEQJf0bLGPASsNQ591Q1ddp462Fm/b3PuyOQgUrDFa87mETqzJ+hBUOAq4HvzGyBt+weoCOAc+4FYCRwo5kVAweAUS5cE8VLzDlYXFph/x/XnkKzlAT6dUrnwKESTrj/07JjCzfks3b7PrJaNQGguKSUP37+A2NP60pa48SQxi0STv6MlpkJ1Nh0cs49BzwXqKBEylu/09OPfvvZx5LdqUWFuWgaJcXzzJV9+fUb8wHYuPsAp/9hBl+NO4O2zVKYvGQrz89YxfY9B3nyst5hiV8kHDQoWCLe737cg0aJ8Yzy3q1ama/rrUMen0a7tBQ25RcCsKqG6YVFYpGSu0S8a4d0rvG4VfPF8nBiB5i3fndAYxKJdJpbRqKenyMlRRoUJXeJehpMI1KVkrtEvdoW4xZpiJTcJepVnp7g9RsGVKmjOWekodEFVYkpY3/UlcFdWzH33rNIiI+jSVI8V/59FgeLS3lu2goGd2vFyR1bhDtMkaBTcpeotyX/AABn92jNuPOOB6Bl0+Sy481SEpm6bBuLcvPpMn8j0+44PRxhioSUumUk6m3fewiAnu3SfB6ftnxb2fbqvH1s21Pos55ILFFyl6g3enAWl/TN5BenZvk8/n8X9WTYcUdmIe3/yFSembqC7XsP8vSUH5i4aDMbdx8IUbQioWHhmgImOzvb5eTkhOW1pWHKGjex2mPNGyey4P5zQhiNSP2Y2VznXHZt9dRyFwF27y8KdwgiAaXkLiISg5TcRURikJK7NFivXz+AL+8aVra/OV8XVSV2+LMSUwczm25mS81siZnd6qOOmdkzZrbSzBaZ2cnBCVfk6L01ZiBrHz+fwd1a0SG9cVn5xEWbKS3VGjMSG/xpuRcDdzjnTgAGAjeZWY9Kdc4Dunt/xgDPBzRKkQD4/U968dLobAZ0aVmh/IWrPG2Rhycu5Ud/mB6O0EQCrtbk7pzb7Jyb593eAywFMitVuwj4p/OYBTQ3s7YBj1bkKFxxSkfOPKF1lfLhPdvy5MiTANiw8wCvzVrHn6esIGftzlCHKBIwdepzN7MsoC8wu9KhTGBDuf1cqv4BEIlYl2V34B/XnALAfR8s5k9TfmDkC9+EOSqR+vM7uZtZU+A94DbnXEHlwz4eUqXz0szGmFmOmeXk5eXVLVKRIGvWqOoC2u/kbPBRUyTy+ZXczSwRT2L/t3PufR9VcoEO5fbbA5sqV3LOjXfOZTvnsjMyMiofFgmreB9LOv3Pu4tYu31fGKIROTr+jJYx4CVgqXPuqWqqTQB+7h01MxDId85tDmCcIkHXu30avTs059K+mXx8y6ll5af/YQYfLazSVhGJaP5M+TsEuBr4zswWeMvuAToCOOdeACYBI4CVwH7g2sCHKhJcZsaHNw0p2x92XAbTl3u6D295Yz7tmjeiaXICny3Zws1ndNMKUBLRNHGYSA1OemAyBYXFVcoXP3guTZO1HIKEniYOEwmA5MR4n+UFB3xPNJY1biJ//Gx5MEMS8YuSu0g9XPLXr1i6uYAt+Z6FPzbs3E9hUQkAz05bGc7QRAAtsydSL1sLDnLen78EYOUj5zH0iemcefwxYY5K5Ai13EWO0sMTlwIwddmR5fy27z0YrnBEACV3kaP2ytdrq5S9Nzc39IGIlKPkLhIEcRomKWGm5C5Sg/qm6L0Hqw6fFAklJXeRIOjRrlm4Q5AGTsld5Cide2LVaYSTEvSrJeGlT6BIDWrrOv/zqD787epsLj05k0v7ZtIuLQWAkhLHhIWbKNHKThImGucuUoPEeN/tny/vGkZm80bEeWeSfOryPgB8l5vPj5+byaOTlrJ6+z4KDhRx1cBOIYtX5DAld5Ea3DC0C7+bsIQRvdowsl97Tu2WQUKclSX1yuK8fwtWe6cJ3r3/UKhCFalAyV2kBqMHZ3H1wE7VJvPKKs8Jf3wbz4XVwqIS3pyznp8N7FTttwGRQFJyF6mFv4kdoHVqCm2apdC/czoTFm7iycnLueG1HA5PvvrAR9+z9vHzgxSpyBFqQogEUIsmScy650zGnNYFgOVb91B5Vm1dZJVQUHIXCYKa7lCdt35XCCORhsqfZfZeNrNtZra4muOnm1m+mS3w/twf+DBFoktcDb9Z6nOXUPDnU/YKMLyWOl865/p4fx46+rBEolt8uZb7u2MHcec5x5btX/3SbJ7Sgh4SZLUmd+fcF8DOEMQiEjMOX4RtkhRPdlY6N5/RnZev8ayMtqewmGemrWTDzv3hDFFiXKC+Hw4ys4Vm9omZnRig5xSJWo28y/MN79m2rOyY1JQKdYY+Mb1s9SaRQAtEcp8HdHLO9QaeBT6orqKZjTGzHDPLycvLC8BLi0Smds0b8fYvB/HIJT3LynpmprHogXO4pG9mWdnkJVsqPG7bnkJmrd4Rsjgldh11cnfOFTjn9nq3JwGJZtaqmrrjnXPZzrnsjIyMo31pkYjWv3M6KZUW2G6WksifrujDKVktALj1zQXMXXdk9Mz5z8xk1PhZIY1TYtNRJ3cza2PmuXpkZv29z6mmh0gNXrtuQNn2T57/mtmrd3DtP+aQt8ezPN+6HfvCFZrECH+GQr4BfAMcZ2a5ZnadmY01s7HeKiOBxWa2EHgGGOVc5ds2RKS8lMR4Hr74SJfNFeNnMX35ka7K2as1hkGOTq3TDzjnrqzl+HPAcwGLSEQoKCwKdwgS5XQ3hUiY1PT19qNFm0MWh8QmJXeRMEmKr36KgoUbdjN92Ta+WaXLV1I/mhVSJEyGdGvFOT1aM7hrSzq1asKw445h8GNT2ZRfCMC1r3wLoFkkpV6U3EXCpH2Lxoz/eXaFsngfrfmHP/6eey/o4fM5PvluMxt3H+D6oV2CEqNELyV3kQj34sw19GqfRvdjUmnXPIVGSfEkJ3jGz9/473kASu5ShZK7SBS49c0FZdv9s9J5e+ygMEYj0UDJXSTKzFm7k08XbyZv75H1Wact28oZx7cOY1QSaZTcRSKIv7f/jf3XvAr7L8xYreQuFWgopEgM6NGuWbhDkAij5C4SAzJSk3l55hqem7Yi3KFIhFC3jEgMKCl1PPTx9wDcfEb3MEcjkUAtd5EYUFKqufqkIrXcRaLIwC7pnN2jDQM6p1NYVELzxomc86cvNEWwVKHkLhIFnrmyLxf0alu2NmtlHyzYFOKIJNIpuYtEkB/3bsfzM1Yx9kdduWpgR9q3aFzrYyr3yKzZvo9GifFs2LWfrhlNSW+SFKRoJZJZuNbVyM7Odjk5OWF5bZFIVVrq2HuomGYpiX4/5sLnZrJj7yE27j7g8/iDF57I6MFZAYpQws3M5jrnsmur589KTC+b2TYzW1zNcTOzZ8xspZktMrOT6xOwiEBcnNUpsQN8eNMQXhxd/e/67yYsOdqwJAr5M1rmFWB4DcfPA7p7f8YAzx99WCLiLzMjvpq+eGm4ak3uzrkvgJoWdLwI+KfzmAU0N7O2gQpQRGpXPrmnJutSmgRmnHsmsKHcfq63rAozG2NmOWaWk5eX56uKiNRDvB1J7vf9uAeDu7akZ+aRKQme+mx5OMKSMApEcvf1fdDnVVrn3HjnXLZzLjsjIyMALy0iAKkpntZ6r8w0Ls/uwOs3DOTjW4aS3akFAM9MW8nd7y+iuKQ0nGFKCAXi+1su0KHcfntAg25FQqhl02QW3n8OaY0rXoy9amAnctbtAuCNORtIa5TEuPOOD0eIEmKBaLlPAH7uHTUzEMh3zmnpdpEQq5zYAS7um8m75Rb2eOG/q8gaN5FwDYGW0PFnKOQbwDfAcWaWa2bXmdlYMxvrrTIJWA2sBP4O/Cpo0YpInWVnpbP28fPp27F5WVnnuyexfe/BMEYlwaabmEQakKxxEyvsf3nXMDqk134XbE227SnkYFHpUT+P+CdgNzGJSOyY+OtT+emAjmX7Q5+YTv6BoqPqpun/yFSGPjE9EOFJACm5izQgJ7ZL49FLevGfXw0uK+v94Gfc+K95/Gd+Lm/OWR/G6CSQlNxFGqC+HVswelCnsv1Pl2zhN28tZNz731U7Rw3Aa7PWMWP5tlCEKEdJyV2kgTLzPWXBkMen8dqsdT4XALnvg8Vc849vgx2aBICSu4hUcd8Hi3l22gr2FBaRv7+oyvFSrfwU8ZTcRcSnp6esoNcDn9H7oc/4aOGmCiNtnp22MoyRiT+U3EUaqLqMkLnljfkV9j9dsiXQ4UiAKbmLSJ11bqUx7ZFOyV1E6iw1ufoFRVbl7WVzfvUjbiQ0NPGzSAN1NJdE4+KMklLnc5GQM//4XwDWPn7+UbyCHC213EXEL8e3SWVgl3QA3pizngufmxnmiKQmarmLNFA1XU/t3zmds044hn6d0umZ2YzkhPiyYwMfncqWgkKWbCoIQZSx58T7P2X04CzuGh7cqZeV3EUaqEPFRxbuuLRvJmmNE+nXqQVnndCalMT4ah+3paCwbPtgcUlQY4xFRSWOUNwmoOQu0kCN6t+Bd+flMvN/h9E2rVGdH9+pZWOOu/fTIEQW2xyOam4ODigld5EGqm/HFqx6dES9H79ux/4K+2u272NHuTninXPVTnHQkDnne23SQNMFVRGpkzduGEhSgid13DysG7ec0Q2AYX+YwcgXvimrV1ik9Vp9cRCSlrtfyd3MhpvZcjNbaWbjfBy/xszyzGyB9+f6wIcqIpFgUNeWvH/jYF6/YQB3nnsccdVkqv2HikMcWXRwzmEhaLvX2i1jZvHAX4Cz8SyG/a2ZTXDOfV+p6lvOuZuDEKOIRJiemWll2/kHqk4sBrA5v5AtBYUkJ8TTMb1xWWu/Onl7DpKRmhzQOCNVpPS59wdWOudWA5jZm8BFQOXkLiIN0OFE/8sfdeFXP+rGpMWbufv977jg2Yrj4Gu6qSln7U5GvvANf/3ZyYzo1Tao8YZbqObT9Ce5ZwIbyu3nAgN81PuJmZ0G/AD8xjm3oXIFMxsDjAHo2LFj5cMiEoVG9mvPRX3akRjvaZkfU03re/2O/WS2aOTzrtYV2/YCMH3ZtthP7hF0QdVXHJX/+HwEZDnnTgKmAK/6eiLn3HjnXLZzLjsjI6NukYpIxDqc2AEaJfkeI3/ak9Ppes8kznrqv9U+3tcCITEpBP0y/iT3XKBDuf32wKbyFZxzO5xzh8dA/R3oF5jwRCTaHB4z3zOzGV+NO4O3fzmowvGV2/Yy5fut5O7yDKXcvvdgWQuy9CgW6o4Gw5/+AghNy92fbplvge5m1hnYCIwCflq+gpm1dc5t9u5eCCwNaJQiEjU6t2rC3HvPIr1JEmbmc9Wm6/+Z4/OxsZ3aYdmWPUCEXFB1zhWb2c3AZCAeeNk5t8TMHgJynHMTgF+b2YVAMbATuCaIMYtIhGvZ9Ei/e0K8/5msTbOUYIQTcSJiKCSAc24SMKlS2f3ltu8G7g5saCISC3xdQK3O375YzV3Dj6/ymB17D5KakljrcEo5Qv9TIhJU8XXsg7j97QUV9otKSun/6FTufGdhhekNolnE3KEqIlJfCXG+08x7Nw7i2Sv7Vin/cMEmikuOTF2wp7CYklLHhIWb6PfwFLbkF7Jy2x4KCn3fPBUNIuWCqohIvTVOPjI08o6zj+XqQZ1okpxAYnwc/Tp5LsB2aNGY3g99Vlav228/Ye69ZzFl6VZen72+wvNtLSjkor98Rc/MZnx8y9CQnUcgRcQFVRGRo5EYH1fj3anlpzIor9/DU3yWHx5Rs3hj9C4WEorZMtUtIyJRpaRUs036Q8ldRCLSKVktfJZvK4i+i6qFRSUVriOEgpK7iESchy9HOt7gAAAJcUlEQVTuyVtjBnHbWd3p0bYZ15/auezYjf+eV+NjnXNMX7YNF0F3ux5/36dc+8q3ZfvqcxeRBmP2PWcSH2e0KncD1G1nHcttZx0LwIsz11R5TNa4iQDMuedMjvHeAPXhgk3c9tYCHr64J1cN7BSCyP3z5YrtZduh+LujlruIRITWzVIqJPbKXruuPyd3bM5HN59a5Vj/R6dSWuooKXXMXrMT8MxhE6mKS4Kf3dVyF5GoMLR7BkO7e2aTXfPYCP41ax33fbik7HiXeyrcRM+ewshdCao4BBeFldxFJOqYGVcPyuLqQVnsPVjMz16czcINuyvU2b3/UJiiq11RCFru6pYRkajWNDmBt8YMrFI+ddk2rvjbN2SNm0jWuIkhv8A6f/0unp+xyuexUIycUctdRKJeSmI8ax4bQVGJIzHeKDhQzBXjvynrfwc4/Q8zeOiinrRonMhJ7ZsHPaarXpzNvkMl/LR/1VXnDim5i4j4x8xISvCMMUxrnMgntw7l1a/XMn/Dbj5csIl1O/Yz+uU5AJzgHV75k37tj/p1txUU8vnSrfy0f8cKd54e/p7w+pz1VR6za3/w58VRt4yIxCQz45ohnfnzqL6sffx8/ufc48qOLd1cwB3vLOTsp/5L1riJPDN1Bde98i3LttR9SoOXv1rLb/+zmBPu/5TTnpheVt6/czpA2YpT5e07GPyLvX613M1sOPBnPIt1vOice7zS8WTgn3iW19sBXOGcWxvYUEVE6u+mYd24aVg31mzfx5Tvt/LIpKVlC3M/9fkPgKefvm/H5vzo2Ax+cnJ7OqQ3rvV5d+7z3DFbWFTK+p37mbliO6d2b0WSd13Yf8+u2nI/cKgkUKdVrVpb7mYWD/wFOA/oAVxpZj0qVbsO2OWc6wb8Cfh9oAMVEQmEzq2acMNpXWiUWHUh7+7HNGX++t08PWUFQ5+YTta4ifzilW95O2dDta3tymPzr3ppNh8u2EjT5OrbzvuLgp/c/Wm59wdWOudWA5jZm8BFwPfl6lwEPODdfhd4zszMRdL9vyIi5aQ1SuRApST71i8H4ZzjzW838OTk5QBMW7aNacu2cde7iwA4JjWZtmkpLMzNr/a5b31zQbXHAA4cCn63jD997pnAhnL7ud4yn3Wcc8VAPtAyEAGKiARD05SqbdvEeKNl02RuGtat2sdt23OwxsTuj+tP7XJUj/eHPy13X1PcVG6R+1MHMxsDjAHo2LHq8CARkVAZM7QL/12Rx/Rl2+iZmcaAzukVulKevbIvt7wxn6sHduLYNqlMW7qV6cvzGNA5nT4dm7PvYDG79hWRu/sA153amf/My+WPl/dhwoKNrN2xnx5tm9G9dVMu+evXvH7DANKbJHH1S3OY+b/DSE6o2iUUaFZbz4mZDQIecM6d692/G8A591i5OpO9db4xswRgC5BRU7dMdna2y8nJCcApiIg0HGY21zmXXVs9f7plvgW6m1lnM0sCRgETKtWZAIz2bo8Epqm/XUQkfGrtlnHOFZvZzcBkPEMhX3bOLTGzh4Ac59wE4CXgNTNbCezE8wdARETCxK9x7s65ScCkSmX3l9suBC4LbGgiIlJfukNVRCQGKbmLiMQgJXcRkRik5C4iEoOU3EVEYlCtNzEF7YXN8oB19Xx4K2B7rbWig84lMsXKucTKeYDO5bBOzrmM2iqFLbkfDTPL8ecOrWigc4lMsXIusXIeoHOpK3XLiIjEICV3EZEYFK3JfXy4AwggnUtkipVziZXzAJ1LnURln7uIiNQsWlvuIiJSg6hL7mY23MyWm9lKMxsX7nj8YWZrzew7M1tgZjnesnQz+9zMVnj/beEtNzN7xnt+i8zs5DDG/bKZbTOzxeXK6hy3mY321l9hZqN9vVaYzuUBM9vofV8WmNmIcsfu9p7LcjM7t1x52D9/ZtbBzKab2VIzW2Jmt3rLo+q9qeE8ou59MbMUM5tjZgu95/Kgt7yzmc32/v++5Z02HTNL9u6v9B7Pqu0c68w5FzU/eKYcXgV0AZKAhUCPcMflR9xrgVaVyp4Axnm3xwG/926PAD7Bs7rVQGB2GOM+DTgZWFzfuIF0YLX33xbe7RYRci4PAHf6qNvD+9lKBjp7P3PxkfL5A9oCJ3u3U4EfvDFH1XtTw3lE3fvi/b9t6t1OBGZ7/6/fBkZ5y18AbvRu/wp4wbs9CnirpnOsT0zR1nIvW6zbOXcIOLxYdzS6CHjVu/0qcHG58n86j1lAczNrG44AnXNf4Jmfv7y6xn0u8LlzbqdzbhfwOTA8+NFXVM25VOci4E3n3EHn3BpgJZ7PXkR8/pxzm51z87zbe4CleNYxjqr3pobzqE7Evi/e/9u93t1E748DzgDe9ZZXfk8Ov1fvAmeamVH9OdZZtCV3fxbrjkQO+MzM5ppnHVmA1s65zeD5kAPHeMsj/RzrGnekn8/N3q6Klw93YxBF5+L9Ot8XT0sxat+bSucBUfi+mFm8mS0AtuH5Q7kK2O2cK/YRV1nM3uP5QEsCeC7Rltz9Wog7Ag1xzp0MnAfcZGan1VA3Ws+xurgj+XyeB7oCfYDNwB+95VFxLmbWFHgPuM05V1BTVR9lEXM+Ps4jKt8X51yJc64P0B5Pa/sEX9W8/wb9XKItuecCHcrttwc2hSkWvznnNnn/3Qb8B88bv/Vwd4v3323e6pF+jnWNO2LPxzm31fsLWQr8nSNffyP+XMwsEU9C/Ldz7n1vcdS9N77OI5rfFwDn3G5gBp4+9+ZmdnjFu/JxlcXsPZ6Gp9swYOcSbcndn8W6I4qZNTGz1MPbwDnAYiouKj4a+NC7PQH4uXeEw0Ag//BX7QhR17gnA+eYWQvv1+tzvGVhV+laxiV43hfwnMso74iGzkB3YA4R8vnz9s2+BCx1zj1V7lBUvTfVnUc0vi9mlmFmzb3bjYCz8FxDmA6M9Far/J4cfq9GAtOc54pqdedYd6G8ohyIHzxX/n/A05/123DH40e8XfBc/V4ILDkcM57+tanACu+/6e7IVfe/eM/vOyA7jLG/gedrcRGeFsV19Ykb+AWeC0MrgWsj6Fxe88a6yPtL1bZc/d96z2U5cF4kff6AU/F8VV8ELPD+jIi296aG84i69wU4CZjvjXkxcL+3vAue5LwSeAdI9panePdXeo93qe0c6/qjO1RFRGJQtHXLiIiIH5TcRURikJK7iEgMUnIXEYlBSu4iIjFIyV1EJAYpuYuIxCAldxGRGPT/MxPt0ayqSEUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6d96d32828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.plot(r0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "A_csc = sc.sparse.csc_matrix(Anp)\n",
    "\n",
    "A_csr = sc.sparse.csr_matrix(Anp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 26.4 s, sys: 29.5 s, total: 56 s\n",
      "Wall time: 16.2 s\n",
      "CPU times: user 22.7 s, sys: 22.7 s, total: 45.4 s\n",
      "Wall time: 12.3 s\n",
      "(array([1.09589189e-06+0.j, 4.38356638e-06+0.j, 9.86301985e-06+0.j,\n",
      "       1.75342463e-05+0.j, 2.73972373e-05+0.j, 3.94519821e-05+0.j]), array([[ 2.70250059e-05+0.j,  5.40499823e-05+0.j, -8.10748994e-05+0.j,\n",
      "        -1.08099728e-04+0.j, -1.35124437e-04+0.j,  1.62148999e-04+0.j],\n",
      "       [ 5.40499823e-05+0.j,  1.08099728e-04+0.j, -1.62148999e-04+0.j,\n",
      "        -2.16197560e-04+0.j, -2.70245173e-04+0.j,  3.24291601e-04+0.j],\n",
      "       [ 8.10748994e-05+0.j,  1.62148999e-04+0.j, -2.43221499e-04+0.j,\n",
      "        -3.24291601e-04+0.j, -4.05358504e-04+0.j,  4.86421409e-04+0.j],\n",
      "       ...,\n",
      "       [ 8.10748994e-05+0.j, -1.62148999e-04+0.j, -2.43221499e-04+0.j,\n",
      "         3.24291601e-04+0.j, -4.05358504e-04+0.j, -4.86421409e-04+0.j],\n",
      "       [ 5.40499823e-05+0.j, -1.08099728e-04+0.j, -1.62148999e-04+0.j,\n",
      "         2.16197560e-04+0.j, -2.70245173e-04+0.j, -3.24291601e-04+0.j],\n",
      "       [ 2.70250059e-05+0.j, -5.40499823e-05+0.j, -8.10748994e-05+0.j,\n",
      "         1.08099728e-04+0.j, -1.35124437e-04+0.j, -1.62148999e-04+0.j]]))\n",
      "(array([1.09589189e-06+0.j, 4.38356638e-06+0.j, 9.86301985e-06+0.j,\n",
      "       1.75342463e-05+0.j, 2.73972373e-05+0.j, 3.94519821e-05+0.j]), array([[-2.70250059e-05+0.j, -5.40499823e-05+0.j, -8.10748994e-05+0.j,\n",
      "        -1.08099728e-04+0.j, -1.35124437e-04+0.j,  1.62148999e-04+0.j],\n",
      "       [-5.40499823e-05+0.j, -1.08099728e-04+0.j, -1.62148999e-04+0.j,\n",
      "        -2.16197560e-04+0.j, -2.70245173e-04+0.j,  3.24291601e-04+0.j],\n",
      "       [-8.10748993e-05+0.j, -1.62148999e-04+0.j, -2.43221499e-04+0.j,\n",
      "        -3.24291601e-04+0.j, -4.05358504e-04+0.j,  4.86421409e-04+0.j],\n",
      "       ...,\n",
      "       [-8.10748994e-05+0.j,  1.62148999e-04+0.j, -2.43221499e-04+0.j,\n",
      "         3.24291601e-04+0.j, -4.05358504e-04+0.j, -4.86421409e-04+0.j],\n",
      "       [-5.40499823e-05+0.j,  1.08099728e-04+0.j, -1.62148999e-04+0.j,\n",
      "         2.16197560e-04+0.j, -2.70245173e-04+0.j, -3.24291601e-04+0.j],\n",
      "       [-2.70250059e-05+0.j,  5.40499823e-05+0.j, -8.10748994e-05+0.j,\n",
      "         1.08099728e-04+0.j, -1.35124437e-04+0.j, -1.62148999e-04+0.j]]))\n"
     ]
    }
   ],
   "source": [
    "# https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.eigs.html\n",
    "# which=\"SM\" -> smalles magnitude\n",
    "\n",
    "# only the 6 smallest Eigenvalues are calculated\n",
    "\n",
    "%time r1 = sc.sparse.linalg.eigs(A_csc, which=\"SM\")\n",
    "%time r2 = sc.sparse.linalg.eigs(A_csr, which=\"SM\")\n",
    "\n",
    "print(r1)\n",
    "print(r2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-1.3004e-16+0.j -4.7375e-17+0.j -6.8432e-17+0.j -2.0634e-17+0.j -4.9589e-17+0.j  2.4645e-17+0.j]\n",
      "[-2.2420e-16+0.j -9.2711e-17+0.j -9.9823e-17+0.j -5.7771e-17+0.j -8.0817e-17+0.j -1.1676e-17+0.j]\n"
     ]
    }
   ],
   "source": [
    "# verify that the results are equivalent\n",
    "\n",
    "print(kk[:6] - r1[0])\n",
    "print(kk[:6] - r2[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Conclusion:** `selpc` is faster approx. by factor 30 than `sc.sparse.linalg` and even more compared to `np.linalg`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}