Python/NAG Part 4 – Structures

April 16th, 2009 | Categories: NAG Library, programming, python | Tags:

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.

No comments yet.