ENH: Add wrapper for LAPACK functions stgsyl and dtgsyl

[chjaesch]

Is your feature request related to a problem? Please describe.

I need the LAPACK functions stgsyl and dtgsyl for a dynamic system analyzing package.

Describe the solution you'd like.

A low-level wrapper in the sub-package scipy.linalg.lapack would be sufficient.

Describe alternatives you've considered.

Currently I use the python-julia interface and the sylvsys function in the MatrixEquations.jl package.
However, a pure Python solution would be advantageous.

Additional context (e.g. screenshots, GIFs)

No response

[mainak33]

Hi. I want to give this one a go. But since it's my first time going through the scipy source, I could use some help. I'm trying to base changes on 2aac3f3. But I have some questions about the *.pyf.src syntax and will likely need some help with testing too. Let me know if anyone can help.

[j-bowhay]

@Kai-Striega maybe you would be able to help here?

[ilayn]

I can also help. The initial steps are writing a <prefix2> wrapper since only s, and d flavors are relevant here. It can go into flapack_other.pyf.src. We will also need a wrapper for <prefix2>tgsyl_lwork to query the lwork parameter.

In the protoargument we <ctype2> and in the body of the wrapper wherever needed this is reflected as <ftype2> to generate wrapper both for s and d in one wrapper. These are the ones that need templating. The definitions of these are found in https://github.com/scipy/scipy/blob/main/scipy/linalg/flapack.pyf.src#L10-L25

A simpler example can be found here

https://github.com/scipy/scipy/blob/main/scipy/linalg/flapack_sym_herm.pyf.src#L239-L259

[Kai-Striega]

@mainak33 I've got two unmerged PRs for you to look at: #11890, #11985 these implement ?gbsvx and ?pbsvx respectively, including tests. These have not been merged because we didn't support the required LAPACK version. This could also be a problem for ?tgsyl (I don't know in what version it was introduced, off the top of my head). According to the Toolchain roadmap we currently use LAPACK 3.7.1. Does ?tgsyl fit within that requirement?

@ilayn LAPACK version 3.7.1 was released 25/06/17 - this seems a bit dated to me, do you think it would be worth bumping the version a bit?

[ilayn]

LAPACK version 3.7.1 was released 25/06/17 - this seems a bit dated to me, do you think it would be worth bumping the version a bit?

Historically, it depends on what Mac/Linux distros offer, by default, but because OpenBLAS has become more prominent, I think that situation got a bit better. I am following if we are missing very crucial ones on the new LAPACK versions but not much has changed substantially

[mainak33]

@ilayn Thanks for the information. I was unaware of the lapack version dependency and had started with the 3.11 doc, You can see the ongoing work here: main...mainak33:scipy:mainak33/ENH18404 . Lapack v3.11 also had c/ztgsyl in addition to s/dtgsyl. So, I was trying to wrap those as well as using instead of . However the IWORK dimension in the real vs complex types is M+N+6 vs M+N+2. Does that mean I have to write separate templates for real/complex versions or is there a way to check the template type within the function? Also, I thought into flapack_sym_herm.pyf.src because in the doc *sytrf function is in the same sections and that wrapper is in this file. I am also not clear about the usage of the _lwork wrappers routines? What is their purpose?

[mainak33]

@mainak33 I've got two unmerged PRs for you to look at: #11890, #11985 these implement ?gbsvx and ?pbsvx respectively, including tests. These have not been merged because we didn't support the required LAPACK version. This could also be a problem for ?tgsyl (I don't know in what version it was introduced, off the top of my head). According to the Toolchain roadmap we currently use LAPACK 3.7.1. Does ?tgsyl fit within that requirement?

@ilayn LAPACK version 3.7.1 was released 25/06/17 - this seems a bit dated to me, do you think it would be worth bumping the version a bit?

Thank you for the info. I'll take a look and maybe find some answers to my questions above. I cannot find the doc for Lapack v3.7.1 easily to see if the *tgsyl functions are there. If the functions are not there, is it worth pursuing a PR so that it can be merged in the future?

[ilayn]

Does that mean I have to write separate templates for real/complex versions

yes that's why we have <prefix2> and <prefix2c> distinctions. It's just painful to do it but unfortunately necessary. Also we hide all the interim variables instead of outputting them such as work and iwork and bunch of other things that it wants to spit out.

Also, I thought into flapack_sym_herm.pyf.src

Lapack classification is a complete chaos. So doxygen ordering is not to be trusted. These arrays are not symmetric so let's bunch them into ...other.pyf.src.

I am also not clear about the usage of the _lwork wrappers routines?

Because LAPACK is based on F77 there is this tendency to send in the variables and get a result in the same variables in place. So if you want to get the machine specific value of lwork you first invoke the function with lwork set to -1. That overrides the first element of work to the correct value you take that value and invoke the function again with lwork set to that value. Ideally it doesn't reference the remaining variables so essentially they can even be NULLs. It's just ... history. But there are also lower bounds to this value that allows the routines to do blocked calls and those are the numbers you are seeing as minimum numbers. So if you don't want to do lwork query then you use at least that much of amount.

But because we are using f2py which automatically generates those interim arrays, just to ask for lwork number it is very wasteful to create and destroy those big arrays. Instead we create tiny array just for LAPACK to fill in and it writes into that array. We read the number and ask again for the actual routine to work with all variables initialized.

You can read off from just above your work for the differences between sytrf wrapper and the sytrf_lwork wrapper.

If the functions are not there,

I just checked, they are quite old functions so no problem in that regard. In fact they are from ?xxsyl family.

[ilayn]

The good ones are here Reference-LAPACK/lapack#651 but they are brand new and we can't jump on to that yet.

[Kai-Striega]

@mainak33 I think this is good to go forward. Do you have enough information to get started? If you still have questions, feel free to ask here or open up a draft PR and I'll take a look.

[mainak33]

@mainak33 I think this is good to go forward. Do you have enough information to get started? If you still have questions, feel free to ask here or open up a draft PR and I'll take a look.

Hi Kai. Thank you @ilayn for your help till now as well. Draft PR at : #18571 . I'm working on it and am almost there. I'm repeating info on the draft PR for convenience below:

I am working on tests right now. In some tests I tried to calculate lhs and rhs of the sylvester equations after solution, but they don't seem to meet the desired accuracy of machine precision*100 that is used in other tests. So, I still have to debug. For now I have these questions:

1. I'm not too sure what the input matrix forms should be. Lapack doc does say (a,b) and (d,e) should be quasi upper triangular and upper triangular respectively. Does quasi upper triangular mean block upper triangular?
2. How do we check the matrix forms like (quasi) triangular inside the template code?
3. I am still not sure how f2py is actually building the *.pyf.src template signature files. Does f2py itself have the mechanism to replace the templates or is there some extra build step somewhere to generate the actual *.pyf from the templates?
4. If you are aware of other sources for tests let me know (NAG is mentioned in code. But not sure where to find it.)

[ilayn]

I'm not too sure what the input matrix forms should be. Lapack doc does say (a,b) and (d,e) should be quasi upper triangular and upper triangular respectively. Does quasi upper triangular mean block upper triangular?

Quasi-upper triangular is triangular matrix and occasional entries sticking out on the first subdiagonal. typical to get for real Schur forms. You can call it block upper triangular but the block size can only be 2 or 1.

How do we check the matrix forms like (quasi) triangular inside the template code?

You don't, it's caller's problem

I am still not sure how f2py is actually building the *.pyf.src template signature files. Does f2py itself have the mechanism to replace the templates or is there some extra build step somewhere to generate the actual *.pyf from the templates?

It's black magic. Don't bother with it. Joking aside it's really complicated code to read so I wouldn't know where to start.

If you are aware of other sources for tests let me know (NAG is mentioned in code. But not sure where to find it.)

https://www.nag.com/numeric/cl/nagdoc_cl26.2/nagdoc_fl26.2/html/f08/f08yhf.html is one example. But you can create problems yourself by taking arbitrary arrays and performing schur on them. There is no point testing for functionality since if something is broken we should take it to LAPACK folks. So please don't overdo it with 1000x1000 arrays.

[Kai-Striega]

It's black magic. Don't bother with it. Joking aside it's really complicated code to read so I wouldn't know where to start.

👍 from me for this, if you're still interested the source code is here: https://github.com/numpy/numpy/tree/main/numpy/f2py

[mainak33]

@ilayn @Kai-Striega I have made the changes requested by ilayn and added tests. See changes in Updated the draft PR branch here: https://github.com/scipy/scipy/pull/18571/files. If it passes the CI, and you are happy with the changes, I can create a real PR from a cleaner branch with a squashed commit. Thanks!

FYI: It seems like the routine doesn't like anything below an absolute tolerance of 1e-2, atleast on my local machine. But on inspection the results look good enough. Also I saw the value in some other test and I think it should be ok. (I don't think anything in scipy is affecting solution precision)

[mainak33]

have made the changes requested by ilayn and added tests. See changes in Updated the draft PR branch here: https://github.com/scipy/scipy/pull/18571/files. If it passes the CI, and you are happy with the changes, I can create a real PR from a cleaner branch with a squashed commit. Thanks!

Hmm looks like some tests are failing on some platforms. Any suggestions?