The three levels of interface in the NAG Library for Python
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. In this way, we can semiautomatically create idiomatic interfaces to our core algorithms, which are mainly written in Fortran and C, in languages such as MATLAB and Python. It is also the technology behind our indevelopment C++ interface.
When we used this technology to reengineer the NAG Library for Python, we provided three levels of interface for each of NAG's 1900+ routines.

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.

base is a pure Python Fortranlike 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.

_primitive is a set of undocumented, direct ctypes wrappers for each routine.
A demonstration: Solving a linear system of equations
To demonstrate the three different interface types available, we will consider the solution of the matrix equation
$Ax = b$
from numpy.random import rand import numpy as np import timeit
We first create a function that creates random example problems for us
def generate_linear_problem(matrix_size=5,seed=2): """Creates example problems for Linear solvers that solve the matrix equation Ax=b""" np.random.seed(seed) A = rand(matrix_size,matrix_size) b = rand(matrix_size,1) return(A,b)
Let's use this function to create an example matrix $A$ and vector $b$
(A,b) = generate_linear_problem(matrix_size=5,seed=2) A
array([[0.4359949 , 0.02592623, 0.54966248, 0.43532239, 0.4203678 ], [0.33033482, 0.20464863, 0.61927097, 0.29965467, 0.26682728], [0.62113383, 0.52914209, 0.13457995, 0.51357812, 0.18443987], [0.78533515, 0.85397529, 0.49423684, 0.84656149, 0.07964548], [0.50524609, 0.0652865 , 0.42812233, 0.09653092, 0.12715997]])
b
array([[0.59674531], [0.226012 ], [0.10694568], [0.22030621], [0.34982629]])
We have $A$ and $b$ so our task is to find $x$. For this we will use NAG's interface to dgesv which is taken from LAPACK. NAG has been a collaborator in the LAPACK project since 1987 so naturally we provide access to all of it from within our library. Furthermore, we offer the NAG Library for Python linked against the Intel MKL as well as our own implementation of LAPACK so you can be sure that you are making use of highly optimized routines.
naginterfaces.library  easy to use and Pythonic
#Generate problem (A,b) = generate_linear_problem(matrix_size = 5,seed=2) # Solve with the easiest to use version of NAG's general linear solver from naginterfaces.library.lapacklin import dgesv [a,ipiv,x] = dgesv(A, b)
As you would expect from a Python function, the input matrices are unchanged. For evidence, here is $A$ after the call to dgesv
A
array([[0.4359949 , 0.02592623, 0.54966248, 0.43532239, 0.4203678 ], [0.33033482, 0.20464863, 0.61927097, 0.29965467, 0.26682728], [0.62113383, 0.52914209, 0.13457995, 0.51357812, 0.18443987], [0.78533515, 0.85397529, 0.49423684, 0.84656149, 0.07964548], [0.50524609, 0.0652865 , 0.42812233, 0.09653092, 0.12715997]])
and here is $b$
b
array([[0.59674531], [0.226012 ], [0.10694568], [0.22030621], [0.34982629]])
The solution is in its own vector, as compared to a Fortranlike interface that might have overwritten one of the input matrices with the solution.
x
array([[ 0.62236818], [1.41921342], [ 0.2098716 ], [ 1.03787095], [0.48761155]])
Finally, let's ensure that if we compute the matrix product $Ax$ then we recover $b$. That is, let's show that we really have solved the problem.
A @ x
array([[0.59674531], [0.226012 ], [0.10694568], [0.22030621], [0.34982629]])
In short, everything works as you would expect a Python function to work. The fact that you are really calling a Fortran routine is completely hidden to you as the user and we have worked to do this with all 1900+ of the routines in the NAG Library.
Let's now go a little deeper into the interfaces available.
naginterfaces.base  more like Fortran
The routines in naginterfaces.base
look a lot more Fortranlike which in this case means that the input matrices and vectors are modifed in place with no explicit return values. They are usually more difficult to use and typically require the user to perform more setup but they can be faster than naginterfaces.library
in certain circumstances.
Recall that to import dgesv
from naginterfaces.library
we did
from naginterfaces.library.lapacklin import dgesv
To import the version of dgesv
from base
, all we need to do is change library
to base
in the import
from naginterfaces.base.lapacklin import dgesv
However, we'll import using the name dgesv_base
to allow us to differentiate from the dgesv
we previously used.
from naginterfaces.base.lapacklin import dgesv as dgesv_base #Generate problem (A,b) = generate_linear_problem(matrix_size = 5,seed = 2) #Have to define more input variables compared to the Library version N = 5 # Size of matrix A ipiv = np.zeros(N,dtype='int64') # space for the pivot indices. Refer to NAG documentation for what these are nrhs = 1 # Number of RHS vectors in b #Note that there is no explicit output. The matrices in the input arguments are overwritten dgesv_base(N,nrhs,A,ipiv,b)
The behaviour of this version of dgesv
is very different from the previous one. Input variables may have been modified in place and we suggest that you refer to the documentation to see which ones might have been.
In this case, $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.
A
array([[ 0.78533515, 0.85397529, 0.49423684, 0.84656149, 0.07964548], [ 0.64335092, 0.48411929, 0.1101546 , 0.4481052 , 0.07591998], [ 0.42062911, 0.3192565 , 0.37621299, 0.08662677, 0.20908812], [ 0.55517049, 0.92575459, 0.46064501, 0.34026769, 0.20955231], [ 0.79091562, 0.30215756, 0.76978663, 0.13548722, 0.2310688 ]])
Similarly, the original input vector $b$ has been overwritten with the solution $x$
b
array([[ 0.62236818], [1.41921342], [ 0.2098716 ], [ 1.03787095], [0.48761155]])
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.
It turns out that this is the case and, for large matrices, it can be demonstrably faster to use the base
interfaces compared to the Library
interfaces.
naginterfaces._primitive  close to the metal
Finally we have the _primitive
interface. As the name implies, these are primitive and difficult to use but possibly useful. These are the ctypes headers to the compiled library functions that are used by the higher level interfaces.
There is no documentation for this interface level!
All you have is the function definitions themselves which you'll find /sitepackages/naginterfaces/_primitive/
within your Python install. For example, The _primitive
version of dgesv
is in /sitepackages/naginterfaces/_primitive/lapacklin.py
and contains the following
dgesv = utils._EngineState._get_symbol('dgesv')
dgesv.restype = None
dgesv.argtypes = (
utils._BYREF_INT, # Input, n.
utils._BYREF_INT, # Input, nrhs.
utils._EngineFloat64ArrayType.c_type, # Input/Output, a[lda, *].
utils._BYREF_INT, # Input, lda.
utils._EngineIntArrayType.c_type, # Output, ipiv[n].
utils._EngineFloat64ArrayType.c_type, # Input/Output, b[ldb, *].
utils._BYREF_INT, # Input, ldb.
utils._BYREF_INT, # Output, info.
)
f07aaf = utils._EngineState._get_symbol('f07aaf')
f07aaf.restype = dgesv.restype
f07aaf.argtypes = dgesv.argtypes
Here's how you use it
from naginterfaces._primitive.lapacklin import dgesv as dgesv_primitive from naginterfaces.base import utils import ctypes #Close to the metal function that requires knowledge of ctypes, extra setup and care that you don't shoot off your own foot #Generate problem (A,b) = generate_linear_problem(matrix_size = 5,seed=2) #Create ctypes input and output arguments N_c = utils.EngineIntCType(N) one_c = utils.EngineIntCType(1) A_c = A.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) ipiv_c = ipiv.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) b_c = b.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) info = ctypes.c_int(0); dgesv_primitive(N_c, one_c ,A_c ,N_c ,ipiv_c, b_c, N_c,info)
The solution is not directly viewable since it's just a pointer:
b_c
<naginterfaces.base.utils.LP_c_double at 0x286aab71948>
We need to convert this pointer to a numpy array if we are to view the results.
np.ctypeslib.as_array(b_c,(N,1))
array([[0.07647678], [ 0.37618264], [ 1.70266994], [0.90186141], [ 0.3097503 ]])
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.
Speed differences between interface levels
The natural question 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
#library interface  Easy to use from pytictoc import TicToc timer = TicToc() matrix_size=12000 (A,b) = generate_linear_problem(matrix_size,seed=2) timer.tic() [a,ipiv,x] = dgesv(A, b) timer.toc()
Elapsed time is 12.050469 seconds.
#base interface  Fortranlike API (A,b) = generate_linear_problem(matrix_size,seed=2) #Have to define more input variables for this version of the API N=matrix_size ipiv = np.zeros(N,dtype='int64'); timer.tic() dgesv_base(N,1,A,ipiv,b) timer.toc()
Elapsed time is 11.606593 seconds.
#_primitive interface  ugly but fast (A,b) = generate_linear_problem(matrix_size,seed=2) N_c = utils.EngineIntCType(N) one_c = utils.EngineIntCType(1) A_c = A.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) ipiv_c = ipiv.ctypes.data_as(ctypes.POINTER(ctypes.c_int)) b_c = b.ctypes.data_as(ctypes.POINTER(ctypes.c_double)) info = ctypes.c_int(0); timer.tic() dgesv_primitive(N_c, one_c ,A_c ,N_c ,ipiv_c, b_c, N_c,info); timer.toc()
Elapsed time is 10.575022 seconds.
In the above run, there is an almost 20% difference between calls to the fastest and slowest versions of the interface. This is just as you might expect since lower interface levels do a lot less for you.
For more robust timings, and additional advice concerning speed and algorithm choice for this particular problem, we refer you to our Jupyter notebook on this subject.
Materials relating to this blog post
 NAG Library for Python examples on GitHub Includes a list of the various areas of mathematical funcitonality provided by the Library with links to the documentation.
 Exploiting Matrix Structure for the solution of linear systems A demonstration of how using specialized solvers for this problem can offer significant speedup compared to the general case.
Add new comment