I a writing a Python module to access fortran library that has elemental functions. I came up with the following code (an example):

The library function is fthe(t, p)

I created a fortran wrapper:

! calculate equivalent potential temperature given temperature and

! pressure

subroutine c_fthe(t, nt, p, np, the, nthe) bind(c)

integer(c_int), intent(in):: nt, np, nthe

real(c_float), intent(in):: t(nt), p(np)

real(c_float), intent(out):: the(nthe)

real(krealfp):: pk(np)

pk = fpkap(p) ! another elemental function, not relevant to the discussion

if (nt.eq.np) then

the = fthe(t, pk)

else if (nt.eq.1) then

the = fthe(t(1), pk)

else if (np.eq.1) then

the = fthe(t, pk(1))

end if

end subroutine

In *.pyx file I have

cdef extern void c_fthe(float *t, int *nt, float *pk, int *npk, float *the, int *nthe)

def pt2the(object p, object t):

'''Compute equivalent potential temperature at the LCL

from pressure and temperature.

Temperature range is 183.16 to 303.16 Kelvin, pressure

4000 to 110000 Pa.

Input argument list:

t: LCL temperature in Kelvin

p: LCL pressure

Returns equivalent potential temperature in Kelvin.

Output values outside valid range are masked.

'''

t = np.clip(t, 183.16, 303.16)

p = np.clip(p, 4.0e3, 1.1e5)

return generic2arg(t, p, &c_fthe)

and the following monstrosity (generic3arg is much worse):

cdef generic2arg(object o1, object o2, f3arg f):

cdef np.ndarray[np.float32_t, ndim=1, mode='c'] o1c, o2c, oc

cdef int n1, n2, n

n1 = o1.size

n2 = o2.size

mask1 = o1.mask if hasattr(o1, 'mask') else None

mask2 = o2.mask if hasattr(o2, 'mask') else None

mask = None

if o1.shape == o2.shape:

n = n1

shape = o1.shape

if mask1 is not None and mask2 is not None:

# if objects are arrays, both must be masked or not masked

mask = mask1 | mask2

elif n1 == 1:

n = n2

shape = o2.shape

mask = mask2

elif n2 == 1:

n = n1

shape = o1.shape

mask = mask1

else:

raise ValueError('Invalid dimensions: %d, %d' % (n1, n2))

if shape == (0,):

return np.array([], dtype=o1.dtype)

o1c = np.ascontiguousarray(o1.ravel(), dtype=np.float32)

o2c = np.ascontiguousarray(o2.ravel(), dtype=np.float32)

oc = np.ascontiguousarray(np.zeros((n,)), dtype=np.float32)

f(&o1c[0], &n1, &o2c[0], &n2, &oc[0], &n)

if mask is not None:

return np.ma.array(oc.reshape(shape), mask=mask)

else:

return oc.reshape(shape)

Is there a better way of handling this? I know C does not know about elemental functions. Googling for "cython fortran elemental functions" produced no hits.

George

--

---

You received this message because you are subscribed to the Google Groups "cython-users" group.

To unsubscribe from this group and stop receiving emails from it, send an email to

cython-users+unsubscribe <at> googlegroups.com.

For more options, visit

https://groups.google.com/d/optout.