Allocate Fortran arrays from Python

[zhucaoxiang]

@jameskermode Hi James, is there a way to allocate Fortran arrays from the Python side? With the native f2py, it is realized by direct assignment in Python. Somehow, this ability was lost when interfacing with the *.fpp files. And I will get an error of ValueError: array is NULL.

Fortran 90 source code:

module mod
  real, allocatable, dimension(:,:) :: b 
contains
  subroutine foo
    integer k
    if (allocated(b)) then
       print*, "b=["
       do k = 1,size(b,1)
          print*, b(k,1:size(b,2))
       enddo
       print*, "]"
    else
       print*, "b is not allocated"
    endif
  end subroutine foo
end module mod

and wrap it using f2py -c -m allocarr allocarr.f90.

In Python:

>>> import allocarr
>>> print(allocarr.mod.__doc__)
b - 'f'-array(-1,-1), not allocated
foo - Function signature:
  foo()

>>> allocarr.mod.foo()  
 b is not allocated
>>> allocarr.mod.b = [[1, 2, 3], [4, 5, 6]]             # allocate/initialize b
>>> allocarr.mod.foo()
 b=[
   1.000000       2.000000       3.000000    
   4.000000       5.000000       6.000000    
 ]
>>> allocarr.mod.b                                      # b is Fortran-contiguous
array([[ 1.,  2.,  3.],
       [ 4.,  5.,  6.]], dtype=float32)
>>> allocarr.mod.b.flags.f_contiguous
True
>>> allocarr.mod.b = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]  # reallocate/initialize b
>>> allocarr.mod.foo()
 b=[
   1.000000       2.000000       3.000000    
   4.000000       5.000000       6.000000    
   7.000000       8.000000       9.000000    
 ]
>>> allocarr.mod.b = None                               # deallocate array
>>> allocarr.mod.foo()
 b is not allocated

[jameskermode]

If Fortran 'owns' the array then the allocate() call needs to be in Fortran - easiest way to achieve this is to wrap constructor/destructor routines.

However, it is possible to allocate a numpy array with Fortran memory ordering e.g. a = np.zeros((3, 2), order='F') and pass this in as an intent in/out argument which can be modified from Fortran. (For 1-D arrays order='C'and order='F' are equivalent so it doesn't matter).

[jameskermode]

Now that I read your code more carefully, it's interesting if this possible with native f2py. Would be interesting to figure out where this goes wrong with f90wrap.

[zhucaoxiang]

@jameskermode Yes, I tested it and it is possible for f2py. This comes from the documentation at the end of https://numpy.org/devdocs/f2py/python-usage.html.

[jameskermode]

As far as I remember f90wrap doesn't treat allocatable arrays specially, maybe that's the problem.

[zhucaoxiang]

@jameskermode It is also possible for f2py-f90wrap. The "problem" happens at generating the .fpp file will treat b array as a function.

    @property
    def b(self):
        """
        Element b ftype=real pytype=float

        Defined at library.fpp line 8

        """
        array_ndim, array_type, array_shape, array_handle = \
            _ExampleArray.f90wrap_library__array__b(f90wrap.runtime.empty_handle)
        if array_handle in self._arrays:
            b = self._arrays[array_handle]
        else:
            b = f90wrap.runtime.get_array(f90wrap.runtime.sizeof_fortran_t,
                                    f90wrap.runtime.empty_handle,
                                    _ExampleArray.f90wrap_library__array__b)
            self._arrays[array_handle] = b
        return b

    @b.setter
    def b(self, b):
        self.b[...] = b

    def __str__(self):
        ret = ['<library>{\n']
        ret.append('    b : ')
        ret.append(repr(self.b))
        ret.append('}')
        return ''.join(ret)

    _dt_array_initialisers = []

[jameskermode]

Aha. Would you be willing to remove the [...] on the line self.b[...] = b and see if this fixes the problem? You can edit the auto-generated code after running f90wrap.

[zhucaoxiang]

@jameskermode By doing so, I got the following error.

In [4]: lib.library.b = [[1, 2, 3], [4, 5, 6]]                                                                            
---------------------------------------------------------------------------
RecursionError                            Traceback (most recent call last)
<ipython-input-4-3b18ccdd94cd> in <module>
----> 1 lib.library.b = [[1, 2, 3], [4, 5, 6]]

~/Documents/Code/TMP/f90wrap_test/ExampleArray.py in b(self, b)
     98     @b.setter
     99     def b(self, b):
--> 100         self.b = b
    101 
    102     def __str__(self):

... last 1 frames repeated, from the frame below ...

~/Documents/Code/TMP/f90wrap_test/ExampleArray.py in b(self, b)
     98     @b.setter
     99     def b(self, b):
--> 100         self.b = b
    101 
    102     def __str__(self):

RecursionError: maximum recursion depth exceeded

[jameskermode]

Ah, ok, yes, I see that. Maybe try removing the setter routine altogether.

[zhucaoxiang]

Then it can not be set.

[jameskermode]

OK, this is a bit fiddly - I think I understand the problem at least, and I'll try to take a look myself, but won't manage it tonight. The raw array lives in self._arrays dictionary, you may be able to experiment with assigning to that.

[zhucaoxiang]

Thanks, I am not in a rush. When trying to access self._arrays, I got ValueError: array is NULL.

[jameskermode]

Yes, you need to get it via the property once first, then there will be a cached reference in self._arrays. Assigning to this may well do what you want.

[jameskermode]

If that does work we could change the setter to do self._arrays[array_handle] = b, but would need to look up the handle somehow.

[baojd42]

Yes, you need to get it via the property once first, then there will be a cached reference in self._arrays. Assigning to this may well do what you want.

I tried to access the property. But on the first access, it is still going to call the get_array method, which gives the error ValueError: array is NULL. @jameskermode