Array ordering when wrapping Fortran for Python using SWIG and numpy.i

author

2018-05-01 16:51

Unresolved SO Question: I have a Fortran subroutine similar to the following:

subroutine fsub(array, dim1, dim2) bind(c)
  use iso_c_binding, only: c_int, c_double
  integer(c_int), intent(in), value:: dim1, dim2
  real(c_double), intent(inout):: array(dim1, dim2)

  array(1, 1) = 1
  array(2, 1) = 2
  array(1, 2) = 100
end subroutine

If I wrap it using SWIG and numpy.i and the following typemaps:

%apply (double* INPLACE_FARRAY2, int DIM1, int DIM2) {(double* array, int dim1, int dim2)}

%inline %{
void fsub(double* array, int dim1, int dim2);
%}

Then I would have to allocate a ‘C’ order array to pass in:

In [1]: import numpy as np; \
        import fmod; \
        arrayF = np.empty((100, 100), dtype=np.float_, order='F'); \
        arrayC = np.empty((100, 100), dtype=np.float_)
In [2]: fmod.fsub(arrayF)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-fa80220ec8a8> in <module>()
----> 1 fmod.fsub(arrayF)

TypeError: Array must be contiguous.  A non-contiguous array was given

In [3]: fmod.fsub(arrayC)
In [4]: arrayC[0, 0]
Out[4]: 1.0
In [5]: arrayC[1, 0]
Out[5]: 2.0
In [6]: arrayC[0, 1]
Out[6]: 100.0

My questions are:

Shouldn’t the Fortran statement array(2, 1) = 2 sets arrayC[0, 1] instead?

If I had INPLACE_ARRAY2 instead of INPLACE_FARRAY2 in the %apply directive, then indeed I would have arrayC[0,1] = 2.0 after the call.

What happens exactly?

Why wasn’t arrayF allowed? If I use f2py, then arrayF must be used and arrayC isn’t allowed, which is intuitive.

Comments