As an IMSL user, you now have a clear path to access your favorite IMSL C functions from Python with the release of the IMSL Python Numerical Library (PyNL) 1.0. PyNL contains a set of Python wrappers for a subset of IMSL C functions and you can create additional wrappers to access all the other functions as needed. In this post, we’ll discuss the key steps in wrapping an IMSL C function using PyNL.

PyNL’s underlying architecture and integration with the IMSL C library makes adding new wrappers straightforward. PyNL also integrates with NumPy and SciPy, and is available with Python 2.7.8 and Python 3.4.2. For dependencies and available platforms, see our website.

**Leveraging an existing PyNL wrapper**

Since many of IMSL C’s 500+ APIs have common design features, the best way to create a new wrapper for an IMSL C function is to leverage an existing wrapper for a similar function, if possible. The current set of wrappers include a representative set of algorithms that serve as reference implementations for most API types. For example, to create a new wrapper for the function `lin_sol_posdef()`

, we can reuse much of the existing wrapper for a similar function, `lin_sol_gen()`

. The existing wrapper is a Python class named LU, so let’s walk through the steps.

**Developing the wrapper **

The recent white paper, **How to access IMSL C functions from Python**, describes wrapping an IMSL C function with a Python function. While Python functions work in most cases, sometimes it’s better to use a Python class so you can maintain state in the Python object, which we’ll show here. The complete wrapper code is provided at the end of this blog for closer analysis.

After you install PyNL, all wrapper code is available as Python source within the “imsl” package in the site-packages area of your Python installation. This source code is available to allow additional wrappers to be added to PyNL.

**Setup**

We recommend creating a new package within the PyNL area in your Python site-packages area to allow multiple user-written wrappers to coexist with the existing PyNL wrappers. Use the following steps to create a new package, where `PYNL_AREA`

refers to the `site-packages/imsl`

folder in your Python installation.

1. Create a new folder named `PYNL_AREA/user_lib`

.

2. Add a file named `__init__.py`

to the `PYNL_AREA/user_lib`

folder, with the following contents:

```
"""Initialize imsl.user_lib package."""
from imsl.user_lib.cholesky import Cholesky
__all__ = ('Cholesky',)
```

3. Add the package name ‘user_lib’ to the file `PYNL_AREA/__init__.py`

.

4. Expose the name of the IMSL C function to be wrapped in `PYNL_AREA/_imsllib.py`

. For example, to expose `imsl_d_lin_sol_posdef()`

, add the following line at the end of the set of other functions that are exposed:

` self.expose("imsl_d_lin_sol_posdef ", None, _err_check, "math") `

This sets up the exception handling of PyNL to throw Python exceptions when errors are detected within IMSL C.

Note: If you are wrapping a C/Stat function, the final argument to `self.expose`

should be “`stat`

” rather than “`math`

“.

**The wrapper**

Again, the IMSL C API for `lin_sol_posdef()`

is similar to the interface for `lin_sol_gen()`

, so we’ll leverage the design and source code of the PyNL LU class to create a PyNL class for `lin_sol_posdef()`

. Note that the existing PyNL LU class supports both real and complex versions of `lin_sol_gen()`

but for the purposes of this posting, we’ll wrap just the real-valued IMSL C function `imsl_d_lin_sol_posdef()`

.

The PyNL LU class has the following characteristics:

• **class LU(a)**

Parameters

a: An NxN array containing the input matrix

Methods

cond() | Compute the L1 norm condition number of the matrix |

factor() | Compute the pivoted LU factorization of the matrix |

factor_full() | Compute the pivoted LU factorization of the matrix |

inv() | Compute the inverse of the matrix |

solve(b[, transpose]) | Solve a system of linear equations using right-hand side b |

These characteristics are similar to those required by the IMSL C function `imsl_d_lin_sol_posdef()`

. Using the PyNL LU class as a basis, we can create a Python class with the following characteristics for `imsl_d_lin_sol_posdef()`

:

• **class Cholesky(a)**

Parameters

a: An NxN array containing the input matrix

Methods

cond() | Compute the L1 norm condition number of the matrix |

factor() | Compute the Cholesky factorization of the matrix |

inv() | Compute the inverse of the matrix |

solve(b) | Solve a system of linear equations using right-hand side b |

This new class shares many of the methods of the LU class as well as the shape of the input array, a. It was created using the LU class source code (located in the `site-packages/imsl/linalg`

folder in your Python installation) as a starting point. Note that the method `LU.factor_full()`

is unique to the LU algorithm so is not needed for the Cholesky class.

Once we have defined the design of the new class, the next steps include:

• Validating the input data types

• Performing data conversion in preparation for calling the IMSL C function

• Setting up arguments and calling the IMSL C function

Let’s look at some code snippets to see how to accomplish these steps.

**Validating input data types**

Since IMSL C is a C library that expects data to be of a specific type, the Python wrapper determines whether the input data can be converted to the target precision. This excerpt is from the wrapper method `Cholesky.__init__`

:

```
# attempt to promote matrix a to a compatible type.
common_type = _numpy.promote_types(_numpy.float64, self._a.dtype)
self._a = _numpy.asarray(self._a, dtype=common_type)
if (not _numpy.issubdtype(self._a.dtype, _numpy.float64)):
raise ValueError("array type {} not supported".format(
self._a.dtype.name))
```

**Performing data conversion in preparation for calling the IMSL C function**

Again, since IMSL C expects specific datatypes, all data passed to the IMSL C function is converted to the desired type using the NumPy method `asarray`

. For example, in the wrapper method `Cholesky.solve`

:

` b = _numpy.asarray(b, dtype=self._a.dtype, order='C')`

**Setting up arguments to pass to the IMSL C function**

The Python wrapper creates an argument list and passes it to the IMSL C function. This argument list corresponds to the set of required and optional arguments necessary to perform the desired action. Note that IMSL C requires a terminating zero at the end of the optional argument list. Consider the method `Cholesky.inv`

:

```
args = []
row_dim = self._a.shape[0]
args.append(row_dim)
args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
args.append(0) # b is ignored
self._inverse = _numpy.empty([row_dim, row_dim],
dtype=self._a.dtype)
args.append(_constants.IMSL_INVERSE_ONLY)
args.append(_constants.IMSL_INVERSE_USER)
args.append(self._inverse.ctypes.data_as(_ctypes.c_void_p))
args.append(0)
func = _imsllib.imsllib.imsl_d_lin_sol_posdef
func(*args)
```

**Sample usage**

Let’s see how you can use the Python interface to the IMSL C function `imsl_d_lin_sol_posdef()`

to solve a system of three linear equations with a symmetric positive definite coefficient matrix. We’ll also demonstrate how to compute the inverse of the coefficient matrix, the *L _{1}* norm condition number, and Cholesky factors. The equations are listed below:

_{x1}− 3

_{x2}+ 2

_{x3}= 27

−3_{x1} + 10_{x2} − 5_{x3} = −78

2_{x1} − 5_{x2} + 6_{x3} = 64

```
""" Cholesky class example. """
import imsl.user_lib.cholesky as chol
import numpy as np
a = [[1.0, -3.0, 2.0], [-3.0, 10.0, -5.0], [2.0, -5.0, 6.0]]
b = [27.0, -78.0, 64.0]
chol = chol.Cholesky(a)
x = chol.solve(b)
print("Solution:")
print(x)
print("")
inv = chol.inv()
print("Inverse:")
print(inv)
print("")
# Compute a*inv(a) (=I)
prod = np.dot(a, inv)
print("a*inv(a):")
print(prod)
print("")
print("L1 condition number:")
print(chol.cond())
print("")
L = chol.factor()
print("Cholesky factors:")
print(L)
print("")
# Compute a-L*trans(L) (=0)
for i in range(x.size):
L[i, i+1:] = 0
print("a - L*trans(L):")
print(a - np.dot(L, np.transpose(L)))
```

**Output**

```
Solution:
[ 1. -4. 7.]
Inverse:
[[ 35. 8. -5.]
[ 8. 2. -1.]
[ -5. -1. 1.]]
a*inv(a):
[[ 1. 0. 0.]
[ 0. 1. 0.]
[ 0. 0. 1.]]
L1 condition number:
673.838709677
Cholesky factors:
[[ 1. -3. 2.]
[-3. 1. 1.]
[ 2. 1. 1.]]
a - L*trans(L):
[[ 0. 0. 0.]
[ 0. 0. 0.]
[ 0. 0. 0.]]
```

**Complete wrapper code for imsl_d_lin_sol_posdef()**

The complete wrapper for the IMSL C function `imsl_d_lin_sol_posdef()`

, using a Python class named Cholesky is provided below.

```
"""Cholesky factorization related class, methods."""
import ctypes as _ctypes
import numpy as _numpy
import imsl._constants as _constants
import imsl._imsllib as _imsllib
class Cholesky():
"""Solve a real symmetric positive definite system of linear equations."""
def __init__(self, a):
"""Instantiate Cholesky class."""
self._a = a
self._factor_factor = None
self._cond = None
self._inverse = None
if self._a is None:
raise TypeError("None not supported")
self._a = _numpy.array(a, order='C')
# attempt to promote matrix a to a compatible type.
common_type = _numpy.promote_types(_numpy.float64, self._a.dtype)
self._a = _numpy.asarray(self._a, dtype=common_type)
if (not _numpy.issubdtype(self._a.dtype, _numpy.float64)):
raise ValueError("array type {} not supported".format(
self._a.dtype.name))
if self._a.ndim != 2:
raise ValueError("array of dimension {} not"
" supported".format(self._a.ndim))
if self._a.size == 0:
raise ValueError("empty array not supported")
if (self._a.shape[0] != self._a.shape[1]):
raise ValueError("input matrix must be square")
def solve(self, b):
r"""Solve a real symmetric positive definite system of linear equations
using right-hand side `b`."""
if b is None:
raise TypeError("None not supported")
b = _numpy.asarray(b, dtype=self._a.dtype, order='C')
if b.ndim != 1:
raise ValueError("array of dimension {} not"
" supported".format(b.ndim))
if b.size != self._a.shape[0]:
raise ValueError("dependent variable values length ({}) does not"
" match coefficient matrix row count"
" ({})".format(b.size, self._a.shape[0]))
args = []
row_dim = self._a.shape[0]
args.append(row_dim)
args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
args.append(b.ctypes.data_as(_ctypes.c_void_p))
solution = _numpy.empty_like(b)
args.append(_constants.IMSL_RETURN_USER)
args.append(solution.ctypes.data_as(_ctypes.c_void_p))
if self._factor_factor is None:
self._factor_factor = _numpy.empty([row_dim, row_dim],
dtype=self._a.dtype)
else:
args.append(_constants.IMSL_SOLVE_ONLY)
args.append(_constants.IMSL_FACTOR_USER)
args.append(self._factor_factor.ctypes.data_as(_ctypes.c_void_p))
if self._cond is None:
args.append(_constants.IMSL_CONDITION)
cond = _ctypes.c_double()
args.append(_ctypes.byref(cond))
args.append(0)
func = _imsllib.imsllib.imsl_d_lin_sol_posdef
func(*args)
if self._cond is None:
self._cond = cond.value
return solution
def factor(self):
"""Compute the Cholesky factorization of the matrix."""
if (self._factor_factor is None or self._cond is None):
args = []
row_dim = self._a.shape[0]
args.append(row_dim)
args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
args.append(0) # b is ignored
self._factor_factor = _numpy.empty([row_dim, row_dim],
dtype=self._a.dtype)
args.append(_constants.IMSL_FACTOR_ONLY)
args.append(_constants.IMSL_FACTOR_USER)
args.append(self._factor_factor.ctypes.data_as(_ctypes.c_void_p))
args.append(_constants.IMSL_CONDITION)
cond = _ctypes.c_double()
args.append(_ctypes.byref(cond))
args.append(0)
func = _imsllib.imsllib.imsl_d_lin_sol_posdef
func(*args)
self._cond = cond.value
return (_numpy.copy(self._factor_factor))
def inv(self):
"""Compute the inverse of the matrix."""
if self._inverse is None:
args = []
row_dim = self._a.shape[0]
args.append(row_dim)
args.append(self._a.ctypes.data_as(_ctypes.c_void_p))
args.append(0) # b is ignored
self._inverse = _numpy.empty([row_dim, row_dim],
dtype=self._a.dtype)
args.append(_constants.IMSL_INVERSE_ONLY)
args.append(_constants.IMSL_INVERSE_USER)
args.append(self._inverse.ctypes.data_as(_ctypes.c_void_p))
args.append(0)
func = _imsllib.imsllib.imsl_d_lin_sol_posdef
func(*args)
return _numpy.copy(self._inverse)
def cond(self):
"""Compute the L_1 norm condition number of the matrix."""
if self._cond is None:
self.factor()
return self._cond
```