f90wrap
f90wrap copied to clipboard
jameskermode
Allocate Fortran arrays from Python
@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
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).
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.
@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.
As far as I remember f90wrap doesn't treat allocatable arrays specially, maybe that's the problem.
@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 = []
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.
@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
Ah, ok, yes, I see that. Maybe try removing the setter routine altogether.
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.
Thanks, I am not in a rush. When trying to access self._arrays, I got ValueError: array is NULL.
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.
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.
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
