[WIP] override XERBLA to raise exception

[bnavigator]

Fixes #126 (eventually)

Trying to wrap the Fortran XERBLA call. Still having some issues with the parameter list.

Python 3.8.2 (default, Apr  8 2020, 14:31:25) 
[GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import slycot
>>> print(slycot._wrapper.xerbla.__doc__)
xerbla(srname,info,raise_xerbla,[raise_xerbla_extra_args])

Wrapper for ``xerbla``.

Parameters
----------
srname : input string(len=-1)
info : input int
raise_xerbla : call-back function

Other Parameters
----------------
raise_xerbla_extra_args : input tuple, optional
    Default: ()

Notes
-----
Call-back functions::

  def raise_xerbla(srname,info): return 
  Required arguments:
    srname : input string(len=-1)
    info : input int

>>> def raise_xerbla_2_arguments(srname, info):
...     raise Exception("If I was called with ({}, {}) you would see this".format(srname, info))
... 
>>> slycot._wrapper.raise_xerbla = raise_xerbla_2_arguments
>>> slycot.sb10fd(1,1,1,0,1,1,[1],[1],[1],[1],1e-3,1)
capi_return is NULL
Call-back cb_raise_xerbla_in___user__routines failed.
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/ben/src/Slycot/_skbuild/linux-x86_64-3.8/cmake-install/slycot/synthesis.py", line 2714, in sb10fd
    out = _wrapper.sb10fd(n, m, np, ncon, nmeas, gamma,
TypeError: raise_xerbla_2_arguments() missing 2 required positional arguments: 'srname' and 'info'

But:

Python 3.8.2 (default, Apr  8 2020, 14:31:25) 
[GCC 9.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import slycot
>>> def raise_no_arguments():
...     raise Exception("no arguments provided")
... 
>>> slycot._wrapper.raise_xerbla = raise_no_arguments
>>> slycot.sb10fd(1,1,1,0,1,1,[1],[1],[1],[1],1e-3,1)
capi_return is NULL
Call-back cb_raise_xerbla_in___user__routines failed.
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/ben/src/Slycot/_skbuild/linux-x86_64-3.8/cmake-install/slycot/synthesis.py", line 2714, in sb10fd
    out = _wrapper.sb10fd(n, m, np, ncon, nmeas, gamma,
  File "<stdin>", line 2, in raise_no_arguments
Exception: no arguments provided
>>> 

[bnavigator]

numpy/numpy#16277

[repagh]

OK, is this anyhow possible? How does the exception propagate back through the Fortran calls? AFAIK Fortran has no exception handling.
Then how are Python exceptions handled, and would this not mess with the call stack?

After some searching, I see that the original LAPACK routine calls STOP on purpose. MKL does not do this. XERBLA seems to be called from just about every SLICOT function, and in MKL it only prints a message then returns, and in the LAPACK implementation it stops.

It seems safer to me to just provide a neutralize Fortran XERBLA implementation for Slycot that behaves in an MKL way.

[bnavigator]

The idea is to use F2PY callback arguments. The Fortran routine XERBLA calls RAISE_XERBLA which is supposed to be the Python function slycot._wrapper.raise_xerbla which can then raise the Exception.

The simpler solution would be to just make XERBLA noop to avoid the annoying print to stdout and exit. We can always let the python code handle the INFO parameter as we already do.

[repagh]

I don't have the insight in Python's exception mechanism to determine if it is even possible and safe to pass an exception with Fortran routines on the F2PY call stack.

If you have the experience for that, by all means, but I would take the safe route.

[bnavigator]

If you have the experience for that,

I don't ;)

You have a point here:

I tested by executing

>>> import slycot
>>> import numpy as np
>>> n=1000; m=3; p=2
>>> A = np.random.randn(n,n); B = np.random.randn(n,m)
>>> C = np.random.randn(p,n); D=np.random.randn(p,m)
>>> for i in range(500):
...     try:
...         slycot.sb10fd(n,m,p,0,1,1,A,B,C,D,1e-3,1)
...     except TypeError:
...         pass
... 

For this test, it doesn't matter if the exception currently raised is the desired SlycotError or the TypeError from my initial post above. This allocates a lot of memory during the loop and never frees it until the python process is ended.

So you are probably right: Raising exceptions in a callback and thus never returning to the Fortran call stack is possible, but maybe not a good idea.

I wonder if there will be any comments by F2PY people, reacting to my question in numpy/numpy#16277