SLICOT routines alter input parameters

[bnavigator]

Opening a new issue because the wrapper function specific issue and the PR are closed already.

Slycot must be more careful with F2PY mappings of in/out parameters. See #200 for discussion.

[KybernetikJo]

I was investigating this quite a bit. The behavior of f2py can be surprising.

I still think

intent(in,out,copy)

is the best option. It makes the API save and clean. However, a downside is that all arrays are copied. These could lead to a drop in runtime performance, especially for very big arrays.

Another option could be to use

intent(in,overwrite)

Now the behavior depends on the memory order of the arrays. In numpy there are two memory types of arrays:

• C-contiguous, data is stored row-wise. This is the default
• FORTRAN-contiguous, data is stored column-wise.
  https://en.wikipedia.org/wiki/Row-_and_column-major_order

The f2py behavior depends on the memory order. The underlying SLICOT routines need FORTRAN-contiguous arrays. Without intent(copy) FORTRAN-contiguous numpy arrays go through without a copy. These arrays are altered by Slycot routines.
C-contiguous arrays will trigger a copy, in order to get FORTRAN-contiguous arrays. These arrays are NOT altered by Slycot routines.
All numpy vectors (pure, row, column) are both FORTRAN-contiguous and C-contiguous at the same time. These vectors are altered by Slycot routines.

There are more things, that must be considered, like data type mismatch (e.g. np.float32 and FORTRAN double precision), which can also trigger a copy.

Maybe I will reach out to the f2py developers at some point, in order to full grasp the details.

I will do a detail analysis at some point in the future and sum it up in a report.

[ilayn]

The Fortran behavior is always pass by reference. Hence you have to give the right order and right dtype to have the bindings to work correctly. You can't escape the in-between copies when you are using Python. Unless you specifically start with Fortran arrays with exact dtypes according to <s, d, c, z> flavors.

If the target routine and the call argument does not match and/or the layout is not Fortran contiguous (note that A[1:3, 1:3] is not contiguous even if A is), f2py up/down casts (if possible) and creates a copy of the array. Even if you write overwrite it will overwrite this interim array and you get the eventual copy back.

This is the price we pay for using Fortran in Python via f2py. SLICOT routines have no inherent performance considerations and typically a copy is a negligible cost compared to other, very often O(n**3) operations.

Like we discussed elsewhere previously, I am seriously considering rewriting SLICOT in a different language possibly in pure Python even which I did some already in harold without sacrificing too much performance. The initial response is always that's a lot of work or that's a bad idea and so on. I am aware. But once you start doing it it gets going, see scipy/scipy#18566 for 30K LOC of F77 already converted.