The three levels of interface in the NAG Library for Python

Posted on
3 May 2019

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 semi-automatically 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 in-development C++ interface.

When we used this technology to re-engineer 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 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.

  • _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$

In [1]:
from numpy.random import rand
import numpy as np
import timeit

We first create a function that creates random example problems for us

In [2]:
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$

In [3]:
(A,b) = generate_linear_problem(matrix_size=5,seed=2)
A
Out[3]:
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]])
In [4]:
b
Out[4]:
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

In [5]:
#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

In [6]:
A
Out[6]:
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$

In [7]:
b
Out[7]:
array([[0.59674531],
       [0.226012  ],
       [0.10694568],
       [0.22030621],
       [0.34982629]])

The solution is in its own vector, as compared to a Fortran-like interface that might have overwritten one of the input matrices with the solution.

In [8]:
x
Out[8]:
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.

In [9]:
A @ x
Out[9]:
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 Fortran-like 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 set-up 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.

In [10]:
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.

In [11]:
A
Out[11]:
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$

In [12]:
b
Out[12]:
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 /site-packages/naginterfaces/_primitive/ within your Python install. For example, The _primitive version of dgesv is in /site-packages/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

In [13]:
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:

In [14]:
b_c
Out[14]:
<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.

In [15]:
np.ctypeslib.as_array(b_c,(N,1))
Out[15]:
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

In [16]:
#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.
In [17]:
#base interface - Fortran-like 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.
In [18]:
#_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

Add new comment