Python/NAG Part 4 – Structures
This is part 4 of a series of articles devoted to demonstrating how to call the Numerical Algorithms Group (NAG) C library from Python. Click here for the index to this series.
The NAG C library makes extensive use of C structures throughout its code base and so we will need to know how to represent them in Python. For example, the NAG routine d01ajc is a general purpose one dimensional numerical integrator and one of it’s arguments is a structure called Nag_QuadProgress with the following prototype (found in the nag_types.h file that comes with the NAG C library)
typedef struct { Integer num_subint; Integer fun_count; double *sub_int_beg_pts; double *sub_int_end_pts; double *sub_int_result; double *sub_int_error; } Nag_QuadProgress;
following the ctypes documentation on Structures and Unions, I came up with the following Python representation of this structure and it seems to work just fine.
class Nag_QuadProgress(Structure): _fields_ = [("num_subint", c_int), ("funcount", c_int), ("sub_int_beg_pts",POINTER(c_double) ), ("sub_int_end_pts",POINTER(c_double)), ("sub_int_result",POINTER(c_double) ), ("sub_int_error",POINTER(c_double) ) ]
This is the only extra information we need in order to write a Python program that calculates the numerical approximation to an integral using the NAG routine d01ajc. Let’s calculate the integral of the function 4/(1.0+x*x) over the interval [0,1]. The code is as follows:
#!/usr/bin/env python #Example 3 using NAG C library and ctypes #d01ajc - one dimensional quadrature #Mike Croucher - April 2008 #Concept introduced : Creating and using Python versions of structures such as Nag_QuadProgress from ctypes import * libnag = cdll.LoadLibrary("/opt/NAG/cllux08dgl/lib/libnagc_nag.so.8") d01ajc = libnag.d01ajc d01ajc.restype=None class Nag_QuadProgress(Structure): _fields_ = [("num_subint", c_int), ("funcount", c_int), ("sub_int_beg_pts",POINTER(c_double) ), ("sub_int_end_pts",POINTER(c_double)), ("sub_int_result",POINTER(c_double) ), ("sub_int_error",POINTER(c_double) ) ] #define python function def py_func(x): return 4./(1.0+x*x) #create the prototype for the c callback function C_FUNC = CFUNCTYPE(c_double,c_double ) #now create the c callback function c_func = C_FUNC(py_func) #set up the problem a = c_double(0.0) b = c_double(1.0) nlimit = c_int(100) epsr = c_double(0.00001) ifail = c_int(0) max_num_subint=c_int(200) nlimit = c_int(0) result=c_double(0.0) epsabs= c_double(0.0) epsrel=c_double(0.0001); abserr=c_double(0.0); qp=Nag_QuadProgress() d01ajc(c_func,a,b,epsabs,epsrel,max_num_subint,byref(result) ,byref(abserr),byref(qp),ifail ) print "Approximation of integral= %s" %(result.value) print "Number of function evaluations=%s" % qp.funcount print "Number of subintervals used= %s" % qp.num_subint
Copy and paste (or download) this to a file called d01ajc.py, make it executable and run it by typing the following at a bash prompt.
chmod +x ./d01ajc.py
./d01ajc.py
If all has gone well then you should get following result
Approximation of integral= 3.14159265359 Number of function evaluations=21 Number of subintervals used= 1
Note that in order to create an instance of the Nag_QuadProgress structure I just do
qp=Nag_QuadProgress()
It is also useful to compare the C prototype of d01ajc
void nag_1d_quad_gen (double (*f)(double x),
double a, double b, double epsabs, double epsrel,
Integer max_num_subint, double *result, double *abserr,
Nag_QuadProgress *qp, NagError *fail)
with the way I called it in Python
d01ajc(c_func,a,b,epsabs,epsrel,max_num_subint,byref(result) ,byref(abserr),byref(qp),ifail )
Note how pointer arguments such as
double *abserr
are called as follows
byref(abserr)
Yet again I have cheated with the NagError parameter and just used c_int(0) instead of a proper NagError structure. NagError will be treated properly in a future article.