Add low-level wrappers for BLAS and LAPACK

[amontoison]

Superseed #1386
cc @ViralBShah @langou

[ViralBShah]

The @ccall macro is bit annoying. I suppose this would have been simpler to do with just ccall. @amontoison Did you auto-generate these from some spec file?

[amontoison]

Hello @ViralBShah, I autogenerated them from a custom Julia script that directly parse the Fortran files in Reference-LAPACK.
I will regenerate them with only ccalldirectly.

Note that that I will visit @langou next month during ten days and we plan to work on that.

[ViralBShah]

There are a few lapack calls with one argument that need fixing too. Perhaps the generator script can be updated? It would be nice to have a gen directory where the generator scripts could be checked in.

In theory, we can use Clang.jl to automatically generate the entire API for CBLAS.

[amontoison]

@ViralBShah I can upload my Fortran parser of F77 but I need to clean it a little bit.
We could the wrappers for CBLAS but we don't use any routine of it in Julia so it will just add more functions in LinearAlgebra.

What would be great to have the "buffer" version of the LAPACK routines with LAPACKE.
But it means a quite big revamp of LinearAlgebra and I am scared that a few people complain.

[ViralBShah]

The reason to do cblas would be then to wrap that into the low level Julia interface and have Linear Algebra use that and delete the fortran bindings altogether.

The LAPACKE style overhaul can easily be done in a new external package.

[codecov]

Codecov Report

❌ Patch coverage is 65.65097% with 124 lines in your changes missing coverage. Please review.
✅ Project coverage is 78.07%. Comparing base (5d05175) to head (fe39ab5).
⚠️ Report is 3 commits behind head on master.

Files with missing lines	Patch %	Lines
src/libblas.jl	61.72%	124 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           master    #1403       +/-   ##
===========================================
- Coverage   93.89%   78.07%   -15.82%     
===========================================
  Files          34       36        +2     
  Lines       15893    20249     +4356     
===========================================
+ Hits        14922    15809      +887     
- Misses        971     4440     +3469     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ViralBShah]

We'll also need to see the impact this has on system image size.

[amontoison]

All wrappers are removed in blas.jl and liblapack.jl so it is more modular now if we want to rely on CBLAS / LAPACKE.
Yes, the system image size could be an issue.

[ViralBShah]

The tests are all passing now. I will check the system size impact. Hopefully it isn’t much.

Given that it is all working now, perhaps not much incentive to change to cblas now, but this puts us in a really good position to do it whenever we need to.

[ViralBShah]

These methods are being generated twice. Should we just edit these by hand or would you like to regenerate the files?

WARNING: Method definition csrot(Any, Any, Any, Any, Any, Any, Any) in module BLAS at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:82 overwritten at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:343.
WARNING: Method definition drot(Any, Any, Any, Any, Any, Any, Any) in module BLAS at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:61 overwritten at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:506.
WARNING: Method definition dsymv(Any, Any, Any, Any, Any, Any, Any, Any, Any, Any) in module BLAS at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:9 overwritten at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:576.
WARNING: Method definition dsyr(Any, Any, Any, Any, Any, Any, Any) in module BLAS at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:40 overwritten at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:583.
WARNING: Method definition srot(Any, Any, Any, Any, Any, Any, Any) in module BLAS at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:68 overwritten at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:770.
WARNING: Method definition ssymv(Any, Any, Any, Any, Any, Any, Any, Any, Any, Any) in module BLAS at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:1 overwritten at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:834.
WARNING: Method definition ssyr(Any, Any, Any, Any, Any, Any, Any) in module BLAS at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:33 overwritten at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:841.
WARNING: Method definition zdrot(Any, Any, Any, Any, Any, Any, Any) in module BLAS at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:75 overwritten at /home/viralbshah/julia/usr/share/julia/stdlib/v1.13/LinearAlgebra/src/libblas.jl:951.

[ViralBShah]

This PR adds 5MB to the system image on Linux. This seems acceptable to me.

My hope is that we can separate out the bits for *, \ and / and relevant BLAS/LAPACK calls to be included in the system image in LinearAlgebraBase, and then not build LinearAlgebra itself into the system image.

[ViralBShah]

Is it possible for the wrapper generator script to annotate the types for all the function calls? Currently these are all Any.

[amontoison]

Hello @ViralBShah, in general we don't type the arguments of the wrappers because it is low-level routines for advance users. I can add types but I am scared to be too much restrictive.

With @langou, we discussed about the Julia wrappers and he thinks that we should interface CBLAS / LAPACKE instead in LinearAlgebra.
It will be also easier to maintain because we can regenerate the wrappers with Clang.jl directly.

Viral, if you want to meet, we are available this week with Julien. We are working together this week on reference-LAPACK.

[ViralBShah]

I am happy to meet, althoughI do not have much more to contribute than what I already have. I was wondering myself if we should use the C interfaces and auto generate.

I think it should be ok to add types to the methods and widen them for things like Integer. The rest have very little room and will be enforced by ccall anyways. Perhaps we get this for free in Clang.jl.

[ViralBShah]

It seems like Apple Accelerate does not use LAPACKE.

[ViralBShah]

The other question is LBT support. However it should be easy to add if C API support is missing.

Perhaps worth considering doing this in an external package?

[stevengj]

If we have these, shouldn't they be called unsafe_foo, since they could crash if passed inconsistent arguments?

Who is this PR for, exactly?

[ViralBShah]

People generally keep asking for a way to call the entire LAPACK and BLAS from Julia - not just our prepackaged factorizations. This adds the full API surface. While it is definitely useful, it is worthwhile to discuss where it should live.

[stevengj]

In my mind, exposing a bunch of unsafe functions like this is not much better than telling them to use ccall. In both cases you need to look up a bunch of low-level arguments and understand low-level types etc. or it will crash.

One alternative would just be to add a @blas_ccall macro that automates calling the correct trampoline, along with documenting a couple of blas types like the correct integer type.

[amontoison]

The complaint of @langou when he tried Julia is that he doesn’t understand why some arguments are exposed and others are not, and I agree with him. The design of the routines / high-level interfaces in LinearAlgebra.LAPACK routines also doesn’t allow reusing buffers and, in the end, we need to do our own ccall in our Julia packages.

Example: https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/main/src/block_krylov_utils.jl#L210-L292

What would be nice is to expose documented high-level wrappers without the prefixes s, d, c, z, with type arguments (like we have right now) but with all arguments. That way, it would be consistent with all calls to reference BLAS/LAPACK.

If we want something with very limited arguments, then we interface them in high-level Julia functions like mul!, qr, lu, etc. But currently we have a hybrid with LinearAlgebra.BLAS and LinearAlgebra.LAPACK, and the solution proposed is adding the low-level wrappers now so that we don’t break Julia users.

Another possible approach could be to keep the old methods (in LinearAlgebra.LAPACK), potentially with a limited number of arguments, to avoid breaking existing code, while providing a clean API with all arguments fully documented and typed.

This would "just" involve duplicating some methods like potrf or geqrf that are already interfaced.

[ViralBShah]

@amontoison How do you feel about a separate package? That package could even use Preferences to pick a BLAS and LAPACK and do other nice things.

[stevengj]

The complaint of @langou when he tried Julia is that he doesn’t understand why some arguments are exposed and others are not, and I agree with him. The design of the routines / high-level interfaces in LinearAlgebra.LAPACK routines also doesn’t allow reusing buffers and, in the end, we need to do our own ccall in our Julia packages.

That's a good argument for adding more high-level methods of these functions that expose more arguments, I agree. I'm not sure what it has to do with the present PR, which is not much better than doing your own ccall?