nag_examples/python_3_interface_levels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The three levels of interface in the NAG Library for Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In an earlier blog post, we discussed the technology we use at NAG that [makes use of our XML documentation to create interfaces to the NAG Library Engine](https://www.nag.co.uk/content/making-new-nag-library-python-usage-tutorials). In this way, we can semi-automatically create idiomatic interfaces to our core algorithms, which are mainly written in Fortran and and C, in languages such as [MATLAB](https://www.nag.co.uk/nag-toolbox-matlab) and [Python](https://www.nag.co.uk/nag-library-python).  It is also the technology behind our [in-development C++ interface](https://www.nag.co.uk/content/creating-c-interfaces-nag-library-work-progress).\n",
    "\n",
    "When we used this technology to re-engineer the [NAG Library for Python](https://www.nag.co.uk/nag-library-python), we provided three levels of interface for each of NAG's 1900+ routines.\n",
    "\n",
    "* **library** is a streamlined Python wrapper that calls the 'base' (described below). This is the version that most people will want to use most of the time.  We are attempting to make it as Pythonic as possible.\n",
    "\n",
    "* **base** is a pure Python Fortran-like API calling the NAG Engine. Typically, input data is type checked, some input constraints are verified, and the Engine is called. Arrays are updated in place as they would be in the Fortran Library.\n",
    "\n",
    "* **_primitive** is a set of undocumented, direct ctypes wrappers for each routine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  A demonstration: Solving a linear system of equations\n",
    "\n",
    "To demonstrate the three different interface types available, we will consider the solution of the matrix equation \n",
    "\n",
    "$Ax = b$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy.random import rand\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import scipy.linalg\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_linear_problem(matrix_size=5,seed=2):\n",
    "    \"\"\"Creates example problems for Linear solvers that solve the matrix equation Ax=b\"\"\"\n",
    "    np.random.seed(seed)\n",
    "    A = rand(matrix_size,matrix_size)  \n",
    "    b = rand(matrix_size,1)\n",
    "    \n",
    "    return(A,b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.4359949 , 0.02592623, 0.54966248, 0.43532239, 0.4203678 ],\n",
       "       [0.33033482, 0.20464863, 0.61927097, 0.29965467, 0.26682728],\n",
       "       [0.62113383, 0.52914209, 0.13457995, 0.51357812, 0.18443987],\n",
       "       [0.78533515, 0.85397529, 0.49423684, 0.84656149, 0.07964548],\n",
       "       [0.50524609, 0.0652865 , 0.42812233, 0.09653092, 0.12715997]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(A,b) = generate_linear_problem(matrix_size=5,seed=2)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.59674531],\n",
       "       [0.226012  ],\n",
       "       [0.10694568],\n",
       "       [0.22030621],\n",
       "       [0.34982629]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have $A$ and $b$ so our task is to find $x$.  For this we will use the NAG routine **dgesv** which is taken from LAPACK"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# naginterfaces.library - easy to use and Pythonic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Generate problem\n",
    "(A,b) = generate_linear_problem(matrix_size = 5,seed=2)\n",
    "\n",
    "# Solve with the easiest to use version of NAG's general linear solver\n",
    "from naginterfaces.library.lapacklin import dgesv\n",
    "[a,ipiv,x] = dgesv(A, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.4359949 , 0.02592623, 0.54966248, 0.43532239, 0.4203678 ],\n",
       "       [0.33033482, 0.20464863, 0.61927097, 0.29965467, 0.26682728],\n",
       "       [0.62113383, 0.52914209, 0.13457995, 0.51357812, 0.18443987],\n",
       "       [0.78533515, 0.85397529, 0.49423684, 0.84656149, 0.07964548],\n",
       "       [0.50524609, 0.0652865 , 0.42812233, 0.09653092, 0.12715997]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# As you might expect, the input matrices are unchanged.  Here is A\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.59674531],\n",
       "       [0.226012  ],\n",
       "       [0.10694568],\n",
       "       [0.22030621],\n",
       "       [0.34982629]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# And here is b\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.62236818],\n",
       "       [-1.41921342],\n",
       "       [ 0.2098716 ],\n",
       "       [ 1.03787095],\n",
       "       [-0.48761155]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The solution is in it's own vector\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.59674531],\n",
       "       [0.226012  ],\n",
       "       [0.10694568],\n",
       "       [0.22030621],\n",
       "       [0.34982629]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Demonstrate that the matrix product Ax = b\n",
    "A @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "#There is documentation for every routine in naginterfaces.library\n",
    "?dgesv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#You can see the source code for the wrapper\n",
    "??dgesv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# naginterfaces.base - More like Fortran"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from naginterfaces.base.lapacklin import dgesv as dgesv_base\n",
    "#Looks more like a Fortran routine in that it overwrites the input matrices \n",
    "#A litttle harder to use and requires more set up by the user but can be faster\n",
    "\n",
    "#Generate problem\n",
    "(A,b) = generate_linear_problem(matrix_size = 5,seed = 2)\n",
    "\n",
    "#Have to define more input variables than the library version\n",
    "N=5  # In the the Fortran-like interface, matrices don't know how big they are\n",
    "ipiv = np.zeros(N,dtype='int64');\n",
    "\n",
    "#Note that there is no explicit output. The matrices in the input arguments are overwritten \n",
    "dgesv_base(N,1,A,ipiv,b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, A is no longer the original matrix A, it now contains the the factors L and U from the factorization A = PLU; the unit diagonal elements of L are not stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.78533515,  0.85397529,  0.49423684,  0.84656149,  0.07964548],\n",
       "       [ 0.64335092, -0.48411929,  0.1101546 , -0.4481052 ,  0.07591998],\n",
       "       [ 0.42062911,  0.3192565 ,  0.37621299,  0.08662677,  0.20908812],\n",
       "       [ 0.55517049,  0.92575459,  0.46064501,  0.34026769,  0.20955231],\n",
       "       [ 0.79091562,  0.30215756, -0.76978663,  0.13548722,  0.2310688 ]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, the original input vector `b` has been overwritten with the solution `x`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.62236818],\n",
       "       [-1.41921342],\n",
       "       [ 0.2098716 ],\n",
       "       [ 1.03787095],\n",
       "       [-0.48761155]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can already imagine that this might be more efficient with respect to memory since we are reusing memory wherever possible rather than allocating new memory for solutions.\n",
    "\n",
    "We will see what this might do for us later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the `naginterfaces.base` versions of all functions are also fully documented"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "?dgesv_base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "??dgesv_base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# naginterfaces._primitive - close to the metal\n",
    "Primitive, and difficult to use but possibly useful.  These are the ctypes headers to the compiled library functions.\n",
    "\n",
    "There is no documentation.\n",
    "\n",
    "All you have is the function definitions themselves at `/site-packages/naginterfaces/_primitive/`.  For example, The `_primitive` version of dgesv is in `/site-packages/naginterfaces/_primitive/lapacklin.py` and contains the following\n",
    "\n",
    "```\n",
    "dgesv = utils._EngineState._get_symbol('dgesv')\n",
    "dgesv.restype = None\n",
    "dgesv.argtypes = (\n",
    "    utils._BYREF_INT, # Input, n.\n",
    "    utils._BYREF_INT, # Input, nrhs.\n",
    "    utils._EngineFloat64ArrayType.c_type, # Input/Output, a[lda, *].\n",
    "    utils._BYREF_INT, # Input, lda.\n",
    "    utils._EngineIntArrayType.c_type, # Output, ipiv[n].\n",
    "    utils._EngineFloat64ArrayType.c_type, # Input/Output, b[ldb, *].\n",
    "    utils._BYREF_INT, # Input, ldb.\n",
    "    utils._BYREF_INT, # Output, info.\n",
    ")\n",
    "f07aaf = utils._EngineState._get_symbol('f07aaf')\n",
    "f07aaf.restype = dgesv.restype\n",
    "f07aaf.argtypes = dgesv.argtypes\n",
    "```\n",
    "\n",
    "Here's how you use it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from naginterfaces._primitive.lapacklin import dgesv as dgesv_primitive\n",
    "from naginterfaces.base import utils\n",
    "import ctypes\n",
    "#Close to the metal function that requires knowledge of ctypes, extra setup and care that you don't shoot off your own foot\n",
    "\n",
    "#Generate problem\n",
    "(A,b) = generate_linear_problem(matrix_size = 5,seed=2)\n",
    "\n",
    "#Create ctypes input and output arguments\n",
    "N_c = utils.EngineIntCType(N)\n",
    "one_c = utils.EngineIntCType(1)\n",
    "A_c = A.ctypes.data_as(ctypes.POINTER(ctypes.c_double))\n",
    "ipiv_c = ipiv.ctypes.data_as(ctypes.POINTER(ctypes.c_int))\n",
    "b_c = b.ctypes.data_as(ctypes.POINTER(ctypes.c_double))\n",
    "info = ctypes.c_int(0);\n",
    "\n",
    "dgesv_primitive(N_c, one_c ,A_c ,N_c ,ipiv_c, b_c, N_c,info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is not directly viewable since it's just a pointer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<naginterfaces.base.utils.LP_c_double at 0x1b50142d6c8>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b_c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to convert this pointer to a numpy array if we are to view the results.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.07647678],\n",
       "       [ 0.37618264],\n",
       "       [ 1.70266994],\n",
       "       [-0.90186141],\n",
       "       [ 0.3097503 ]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.ctypeslib.as_array(b_c,(N,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There may be situations, however, where we might prefer to omit the final step. For example, if we want to pipe the result into additional _primitive functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Speed differences between interface levels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The natural qestion that arises when considering these different interfaces regards speed. One may expect that the closer you are to the pure, compiled Fortran code, the faster you'll get the result.  In the case of **dgesv** this turns out to be true  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 10.998698 seconds.\n"
     ]
    }
   ],
   "source": [
    "#library interface - Easy to use\n",
    "from pytictoc import TicToc\n",
    "timer = TicToc()\n",
    "matrix_size=12000\n",
    "(A,b) = generate_linear_problem(matrix_size,seed=2)\n",
    "timer.tic()\n",
    "[a,ipiv,x] = dgesv(A, b)\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 10.422995 seconds.\n"
     ]
    }
   ],
   "source": [
    "#base interface - Fortran-like API\n",
    "(A,b) = generate_linear_problem(matrix_size,seed=2)\n",
    "#Have to define more input variables for this version of the API\n",
    "N=matrix_size\n",
    "ipiv = np.zeros(N,dtype='int64');\n",
    "\n",
    "timer.tic()\n",
    "dgesv_base(N,1,A,ipiv,b)\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 9.022218 seconds.\n"
     ]
    }
   ],
   "source": [
    "#_primitive interface - ugly but fast\n",
    "(A,b) = generate_linear_problem(matrix_size,seed=2)\n",
    "\n",
    "N_c = utils.EngineIntCType(N)\n",
    "one_c = utils.EngineIntCType(1)\n",
    "A_c = A.ctypes.data_as(ctypes.POINTER(ctypes.c_double))\n",
    "ipiv_c = ipiv.ctypes.data_as(ctypes.POINTER(ctypes.c_int))\n",
    "b_c = b.ctypes.data_as(ctypes.POINTER(ctypes.c_double))\n",
    "info = ctypes.c_int(0);\n",
    "\n",
    "timer.tic()\n",
    "dgesv_primitive(N_c, one_c ,A_c ,N_c ,ipiv_c, b_c, N_c,info);\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above run, there is almost 20% difference between the fastest and slowest versions of the interface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Another trick for speed - Fortran arrays"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Numpy arrays have C-type storage by default.  That is, the **rows** are contiguous in memory.  We can see this with the **flags** attribute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  C_CONTIGUOUS : True\n",
       "  F_CONTIGUOUS : False\n",
       "  OWNDATA : True\n",
       "  WRITEABLE : True\n",
       "  ALIGNED : True\n",
       "  WRITEBACKIFCOPY : False\n",
       "  UPDATEIFCOPY : False"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mymat = np.zeros((2,2))\n",
    "mymat.flags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "LAPACK routines such as **dgesv** require that input matrices are stored with Fortran-type storage where the **columns** are contiguous in memory.  If they get the wrong type of matrix, they need to perform a potentially costly transpose step.  \n",
    "\n",
    "We can convert numpy arrays to follow the Fortran convetion as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  C_CONTIGUOUS : False\n",
       "  F_CONTIGUOUS : True\n",
       "  OWNDATA : True\n",
       "  WRITEABLE : True\n",
       "  ALIGNED : True\n",
       "  WRITEBACKIFCOPY : False\n",
       "  UPDATEIFCOPY : False"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mymat = np.asfortranarray(mymat)\n",
    "mymat.flags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This conversion takes time of course but it may be worth it if your algorithm contains many matrix operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "matrix_size=12000\n",
    "#Create C-ordered problem\n",
    "(A_c,b_c) = generate_linear_problem(matrix_size,seed=2)\n",
    "# Create Fortran ordered versions\n",
    "A_f = np.asfortranarray(A_c);\n",
    "b_f = np.asfortranarray(b_c);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 11.171874 seconds.\n"
     ]
    }
   ],
   "source": [
    "from pytictoc import TicToc\n",
    "timer = TicToc()\n",
    "\n",
    "timer.tic()\n",
    "[a,ipiv,x] = dgesv(A_c, b_c) # C-ordered input matrices\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 8.956717 seconds.\n"
     ]
    }
   ],
   "source": [
    "timer.tic()\n",
    "[a,ipiv,x] = dgesv(A_f, b_f) # Fortran-ordered matrices\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NAG has all of LAPACK - Mixed precision solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Along with its own algorithmic content, The NAG library for Python gives us access to all of LAPACK including **dsgesv** which uses single precision internally whenever it can while aiming for full double precision in the output.\n",
    "This can be significantly faster in some cases.\n",
    "\n",
    "NAG is currently actively researching new mixed precision routines in collaboration with the University of Manchester "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 5.870411 seconds.\n"
     ]
    }
   ],
   "source": [
    "from naginterfaces.library.lapacklin import dsgesv\n",
    "(A,b) = generate_linear_problem(matrix_size,seed=2)\n",
    "# Chnage to Fortran ordering\n",
    "A = np.asfortranarray(A);\n",
    "b = np.asfortranarray(b);\n",
    "\n",
    "timer.tic()\n",
    "[a,ipiv,x,itera] = dsgesv(A, b)\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also have access to the **base** and **_primitive** interfaces to **dsgesv** so we can try those out too."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 5.148290 seconds.\n"
     ]
    }
   ],
   "source": [
    "from naginterfaces.base.lapacklin import dsgesv as dsgesv_base\n",
    "(A,b) = generate_linear_problem(matrix_size,seed=2)\n",
    "# Chnage to Fortran ordering\n",
    "A = np.asfortranarray(A);\n",
    "b = np.asfortranarray(b);\n",
    "\n",
    "N=matrix_size\n",
    "ipiv = np.zeros(N,dtype='int64');\n",
    "x = np.zeros((N,1))\n",
    "\n",
    "timer.tic()\n",
    "dsgesv_base(N, 1, A, ipiv, b, x)\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 5.247763 seconds.\n"
     ]
    }
   ],
   "source": [
    "from naginterfaces._primitive.lapacklin import dsgesv as dsgesv_primitive\n",
    "from naginterfaces.base import utils\n",
    "import ctypes\n",
    "\n",
    "(A,b) = generate_linear_problem(matrix_size,seed=2)\n",
    "# Chnage to Fortran ordering\n",
    "A = np.asfortranarray(A);\n",
    "b = np.asfortranarray(b);\n",
    "\n",
    "N = matrix_size\n",
    "x = np.zeros((N,1))\n",
    "ipiv = np.zeros(N,dtype='int64');\n",
    "work = np.zeros((matrix_size,1))\n",
    "swork = np.zeros((matrix_size*(matrix_size+1),1),dtype=np.float32)\n",
    "work = np.asfortranarray(work);\n",
    "swork = np.asfortranarray(swork);\n",
    "\n",
    "N_c = utils.EngineIntCType(N)\n",
    "one_c = utils.EngineIntCType(1)\n",
    "A_c = A.ctypes.data_as(ctypes.POINTER(ctypes.c_double))\n",
    "ipiv_c = ipiv.ctypes.data_as(ctypes.POINTER(ctypes.c_int))\n",
    "b_c = b.ctypes.data_as(ctypes.POINTER(ctypes.c_double))\n",
    "x_c = x.ctypes.data_as(ctypes.POINTER(ctypes.c_double))\n",
    "work_c = work.ctypes.data_as(ctypes.POINTER(ctypes.c_double))\n",
    "swork_c = swork.ctypes.data_as(ctypes.POINTER(ctypes.c_float))\n",
    "info = ctypes.c_int(0);\n",
    "itera = ctypes.c_int(0);\n",
    "\n",
    "timer.tic()\n",
    "dsgesv_primitive(N_c, one_c ,A_c ,N_c ,ipiv_c, b_c, N_c,x_c,N_c,work_c,swork_c,itera,info);\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above, you'll note that the **_primitive** interface is slightly slower than the **base**.  This is because the overhead in **base** is so small for this routine that what we seeing here is the natural small variation between individual runs.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# More robust timings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The timings that we have demonstrated above are not very robust since many factors can affect a single computation. For example, a laptop CPU may only be able to maintain its fastest clock speed for a short time before it has to slow down for thermal reasons.  There may be other programs running or operating system tasks happening in the background that affect individual timings.\n",
    "\n",
    "In an attempt to control for this kind of thing, there is a `compare_solvers` routine in the same directory as this notebook that uses Python's timieit module for more robust timings.\n",
    "\n",
    "Below, we run each computation 10 times and report the average result in seconds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>MatSize</th>\n",
       "      <th>Scipy</th>\n",
       "      <th>NAG:dgesv</th>\n",
       "      <th>NAG:dgesv_base</th>\n",
       "      <th>NAG:dgesv_primitive</th>\n",
       "      <th>NAG:dsgesv</th>\n",
       "      <th>NAG:dsgesv_base</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10000.0</td>\n",
       "      <td>8.672804</td>\n",
       "      <td>8.565575</td>\n",
       "      <td>8.277488</td>\n",
       "      <td>7.011921</td>\n",
       "      <td>5.623696</td>\n",
       "      <td>4.727431</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   MatSize     Scipy  NAG:dgesv  NAG:dgesv_base  NAG:dgesv_primitive  \\\n",
       "0  10000.0  8.672804   8.565575        8.277488             7.011921   \n",
       "\n",
       "   NAG:dsgesv  NAG:dsgesv_base  \n",
       "0    5.623696         4.727431  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = [compare_solvers(mat_size,number=10,fortran_order=False,res_form='raw') for mat_size in [10000]]\n",
    "result = pd.DataFrame(result)\n",
    "result.columns=[\"MatSize\",\"Scipy\",\"NAG:dgesv\",\"NAG:dgesv_base\",\"NAG:dgesv_primitive\",\"NAG:dsgesv\",\"NAG:dsgesv_base\"]\n",
    "result.set_index('MatSize')\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we try the same again with Fortran-ordered matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>MatSize</th>\n",
       "      <th>Scipy</th>\n",
       "      <th>NAG:dgesv</th>\n",
       "      <th>NAG:dgesv_base</th>\n",
       "      <th>NAG:dgesv_primitive</th>\n",
       "      <th>NAG:dsgesv</th>\n",
       "      <th>NAG:dsgesv_base</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>10000.0</td>\n",
       "      <td>8.188727</td>\n",
       "      <td>7.889164</td>\n",
       "      <td>7.17841</td>\n",
       "      <td>7.992964</td>\n",
       "      <td>4.273695</td>\n",
       "      <td>4.641136</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   MatSize     Scipy  NAG:dgesv  NAG:dgesv_base  NAG:dgesv_primitive  \\\n",
       "0  10000.0  8.188727   7.889164         7.17841             7.992964   \n",
       "\n",
       "   NAG:dsgesv  NAG:dsgesv_base  \n",
       "0    4.273695         4.641136  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = [compare_solvers(10000,number=10,fortran_order=True,res_form='raw')]\n",
    "result = pd.DataFrame(result)\n",
    "result.columns=[\"MatSize\",\"Scipy\",\"NAG:dgesv\",\"NAG:dgesv_base\",\"NAG:dgesv_primitive\",\"NAG:dsgesv\",\"NAG:dsgesv_base\"]\n",
    "result.set_index('MatSize')\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Making use of Matrix structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NAG library often contains specialised varations for any given mathematical problem.  In the case of linear solvers, one set of variations caters for coefficient matrices that have a certain structure.  When your problem exhibits such structure, it often pays divideds to use the specialised solver.\n",
    "\n",
    "## Solving linear systems of Toeplitz matrices - An example of exploiting stucture\n",
    "\n",
    "Toeplitz matrices occur in various applications. A couple of examples are given at http://jack.valmadre.net/notes/2015/03/28/symmetric-positive-toeplitz/\n",
    "\n",
    "Here we solve a linear system where the coefficient matrix is a Toeplitz matrix.  First, we make no use of the underlying structure and just use NAG's **dgesv**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 12.231292 seconds.\n"
     ]
    }
   ],
   "source": [
    "from naginterfaces.library.lapacklin import dgesv   # A general solver\n",
    "from naginterfaces.library.linsys import real_toeplitz_solve # A toeplitz solver\n",
    "from pytictoc import TicToc\n",
    "import numpy as np\n",
    "import scipy.linalg\n",
    "\n",
    "timer = TicToc()\n",
    "\n",
    "# Construct a real, symmetric, positive definite toeplitz matrix \n",
    "matrix_size = 5000\n",
    "t = np.arange(0,matrix_size);\n",
    "a = np.exp(-np.abs(t)/10);\n",
    "# The toeplitz matrix is defined by it's diagonals.  We can construct the full matrix from the diagonals using scipy\n",
    "A = scipy.linalg.toeplitz(a, a)\n",
    "# Construct and a random Righ hand side\n",
    "np.random.seed(2)\n",
    "b = np.random.rand(matrix_size,1)\n",
    "\n",
    "timer.tic()\n",
    "[asol,ipiv,x_nag_gen] = dgesv(A, b)\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now solve the same problem but make use of the function **real_toeplitz_solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 0.068450 seconds.\n"
     ]
    }
   ],
   "source": [
    "# NAG's toeplitz solver requires that b be of dimension N instead of N,1\n",
    "# So we reconstruct it with the correct dimension\n",
    "np.random.seed(2)\n",
    "b = np.random.rand(matrix_size)\n",
    "\n",
    "timer.tic()\n",
    "[x_nag_toe,p] = real_toeplitz_solve(matrix_size,a,b,wantp=False)\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "127.17948717948718"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's see how many times faster it is to use the toeplitz solver\n",
    "4.96/0.039"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is **well over 100x faster** to use the Toeplitz solver.  We should check that we get the same answers?  \n",
    "Annoyingly, the solution matrices are are different shapes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5000, 1)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Solution from the general solver\n",
    "x_nag_gen.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5000,)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Solution from the toeplitz solver\n",
    "x_nag_toe.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but the values agree to high accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1191048088221578e-13"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(abs(x_nag_gen - np.reshape(x_nag_toe,(matrix_size,1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other linear solvers in the NAG library for Python that make use of various matrix structures \n",
    "\n",
    "Your application may not be able to make use of Toeplitz solvers but it may well be able to take advantage of one of the other specialised linear solvers in the NAG library. In addition to the mathematical structure of the matrix, other specialisations include the use of **packed storage** and **mixed precision arithmetic**. If you have multiple right-hand sides to compute with the same coefficient matrix **A**, there are some specialised solvers for that situation too.  Here are **some** of the available solvers with links to their Python documentation:\n",
    "\n",
    "\n",
    "### Real matrices\n",
    "\n",
    "* [library.linsys.real_toeplitz_solve](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.linsys.html#naginterfaces.library.linsys.real_toeplitz_solve) Solution of **real symmetric positive definite Toeplitz** system of linear equations.\n",
    "* [library.lapacklin.dgbsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dgbsv) Computes the solution to a **real banded** system of linear equations.\n",
    "* [library.lapacklin.dsgesv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dsgesv) Computes the solution to a **real** system of linear equations using **mixed precision arithmetic**.\n",
    "* [library.lapacklin.dsysv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dsysv) Computes the solution to a **real symmetric** system of linear equations.\n",
    "* [library.lapacklin.dspsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dspsv) Computes the solution to a **real symmetric** system of linear equations, **packed storage**.\n",
    "* [library.lapacklin.dpbsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dpbsv) Computes the solution to a **real symmetric positive definite banded** system of linear equations.\n",
    "* [library.lapacklin.dposv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dposv) Computes the solution to a **real symmetric positive definite** system of linear equations.\n",
    "* [library.lapacklin.dppsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dppsv) Computes the solution to a **real symmetric** positive definite system of linear equations, **packed storage**.\n",
    "* [library.lapacklin.dsposv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dsposv) Computes the solution to a **real symmetric positive definite** system of linear equations using **mixed precision arithmetic**.\n",
    "* [library.lapacklin.dptsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dptsv) Computes the solution to a **real symmetric positive definite tridiagonal** system of linear equations.\n",
    "* [library.lapacklin.dtbtrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dtbtrs) Solution of **real band triangular** system of linear equations, **multiple right-hand sides**.\n",
    "* [library.lapacklin.dtrtrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dtrtrs) Solution of **real triangular** system of **linear equations, multiple right-hand sides**.\n",
    "* [library.lapacklin.dtptrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dtptrs) Solution of **real triangular** system of linear equations, **multiple right-hand sides, packed storage**.\n",
    "* [library.lapacklin.dgtsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dgtsv) Computes the solution to a **real tridiagonal** system of linear equations.\n",
    "\n",
    "### Complex matrices\n",
    "\n",
    "* [library.lapacklin.zgbsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zgbsv) Computes the solution to a **complex banded** system of linear equations.\n",
    "* [library.lapacklin.zhesv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zhesv) Computes the solution to a **complex Hermitian** system of linear equations.\n",
    "* [library.lapacklin.zhpsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zhpsv) Computes the solution to  **complex Hermitian** of linear equations, **packed storage.**\n",
    "* [library.lapacklin.zpbsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zpbsv) Computes the solution to a **complex Hermitian positive definite banded** system of linear equations.\n",
    "* [library.lapacklin.zposv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zposv) Computes the solution to a **complex Hermitian positive definite** system of linear equations.\n",
    "* [library.lapacklin.zppsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zppsv) Computes the solution to a **complex Hermitian positive definite** system of linear equations, **packed storage**.\n",
    "* [library.lapacklin.zcposv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zcposv) Computes the solution to a **complex Hermitian positive definite** system of linear equations using **mixed precision** arithmetic.\n",
    "* [library.lapacklin.zptsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zptsv) Computes the solution to a **complex Hermitian positive definite tridiagonal** system of linear equations.\n",
    "* [library.lapacklin.zgesv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zgesv) Computes the solution to a **complex** system of linear equations.\n",
    "* [library.lapacklin.zcgesv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zcgesv) Computes the solution to a **complex** system of linear equations using **mixed precision** arithmetic.\n",
    "* [library.lapacklin.zsysv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zsysv) Computes the solution to a **complex symmetric** system of linear equations. \n",
    "* [library.lapacklin.zspsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zspsv) Computes the solution to a **complex symmetric** system of linear equations, **packed storage**.\n",
    "* [ibrary.lapacklin.ztbtrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.ztbtrs) Solution of **complex band triangular** system of linear equations, **multiple right-hand sides**.\n",
    "* [library.lapacklin.ztrtrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.ztrtrs) Solution of **complex triangular** system of linear equations, **multiple right-hand sides**.\n",
    "* [library.lapacklin.ztptrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.ztptrs) Solution of **complex triangular** system of linear equations, **multiple right-hand sides**, **packed storage**.\n",
    "* [library.lapacklin.zgtsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zgtsv) Computes the solution to a **complex tridiagonal** system of linear equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}