ENH: linalg: Add wrapper for ?tgsyl

[mainak33]

Reference issue

Closes #18404

What does this implement/fix?

This adds a low-level wrapper for Lapack functions <s/d>tgsyl

Additional information:

• Not final code. Help and review welcome
• 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]

Thanks @mainak33 please see my comments sorry again for not looking at it earlier about lwork.

[ilayn]

Looks almost done to me. However the test need a full PEP8 pass then this is getting ready for merge.

[mainak33]

Looks almost done to me. However the test need a full PEP8 pass then this is getting ready for merge.

Fixed PEP8 styling issues with new tests in latest checkin

[Kai-Striega]

Looks much better after the PEP8 changes. Thank you

[Kai-Striega]

I'm happy with your work, and I believe the failing tests are unrelated to your changes. Well done @mainak33! Let's wait for a final pass by @ilayn and then it should be fine

[ilayn]

If you separate the spectrum a bit further then float32 should also get ~1e-4 level precision. So something like this should make it fly. Nevermind the scale comment above.

    atol = 1e-4 if dtype == np.float32 else 1e-10

    m, n = 10, 15

    a, _ = schur(generate_random_dtype_array([m, m], dtype=dtype))
    d = np.triu(generate_random_dtype_array([m, m], dtype=dtype))
    b, _ = schur(generate_random_dtype_array([n, n], dtype=dtype))
    e = np.triu(generate_random_dtype_array([n, n], dtype=dtype))
    c = generate_random_dtype_array([m, n], dtype=dtype)
    f = generate_random_dtype_array([m, n], dtype=dtype)

    # Adjust the spectrum of the matrices to get a more precise solution
    adjustment = np.eye(m, dtype=dtype)
    a += adjustment
    d -= 5 * adjustment
    b += adjustment
    e -= 5 * adjustment

[mainak33]

If you separate the spectrum a bit further then float32 should also get ~1e-4 level precision. So something like this should make it fly...

I had to make a few minor changes to your suggested code due to matrix sizes but it worked on my local machine. I pushed it. Let's see if it passes all CI tests.

[mainak33]

If you separate the spectrum a bit further then float32 should also get ~1e-4 level precision. So something like this should make it fly...

I had to make a few minor changes to your suggested code due to matrix sizes but it worked on my local machine. I pushed it. Let's see if it passes all CI tests.

Tests are failing on some platforms for atol=1e-4 for single precision. Trying with atol =1e-3

[ilayn]

Well these are tricky problems anyways so if it works with 1e-3 it is what it is, no need to dwell on it. If things pass, let's get it in instead of wasting your time.

[ilayn]

Ah OK, the transposed problem is a different one to separate, hence the artificial push is making the 'T' option solution worse. That's why we get that silly result in float32 probably. Almost there but I'll check this tomorrow.

[mainak33]

Ah OK, the transposed problem is a different one to separate, hence the artificial push is making the 'T' option solution worse. That's why we get that silly result in float32 probably. Almost there but I'll check this tomorrow.

Thanks for looking into it. I'm afraid I don't really know the specifics of this algorithm, so won't be much help figuring out the tolerances beyond trial and error. I'll leave it upto you to fine tune the test and tolerances. Let me know if I can help otherwise or if I need to take any action once we are ready to merge the PR.

[mainak33]

The one failing test on macOS tests (meson) / Meson build (3.11) (https://github.com/scipy/scipy/actions/runs/5162854702/jobs/9300833616?pr=18571) looks like its related to #18601 ?

[ilayn]

Ah dammit. I'm stupid. That's why indeed. Anyways let's wait for that one. This one is indeed ready with smaller tolerances

[mainak33]

Ah dammit. I'm stupid. That's why indeed. Anyways let's wait for that one. This one is indeed ready with smaller tolerances

No worries. I agree. This one looks ready pending the build fix for macOS. Thanks for all the help!

[mainak33]

@ilayn @Kai-Striega I just rebased and pushed and it is now passing all tests with tighter tolerances that Ilayn worked on (Looks like the BLAS package updates for #18601 must have worked to fix the OSX tests for this as well). I think this is now ready for merge.

[Kai-Striega]

I'm happy with this. Waiting on @ilayn for a final review and merge.

[ilayn]

Looks good, thank you very much @mainak33 for your patience and @Kai-Striega for the reviews. Let's click it