Rewrite LAPACK so it's not in F77

@ilayn

From @ilayn, who can hopefully expand:

We must (full emphasis on must) port LAPACK out of F77 for fast computing as scientific community (of all languages). LAPACK codebase is treated as the scripture, and for a good reason because the folks who wrote the algorithms have immense insight and expertise. But that does not mean they coded the most optimal code and also F77 does not lend itself to anything modern hardware offers. You get only what gfortran offers for the loop optimization. While the algorithms are rock solid, they have no way other than writing for loops as if there is no tomorrow. I know this because I spent reading similar code for the good chunk of last year due to META: FORTRAN Code inventory · Issue #18566 · scipy/scipy · GitHub. LAPACK routines are often older than SciPy’s *PACK libraries.

Indeed this needs some explanation because 99% of the responses will be "yeah right, sure" which is not the audience I care about so I don't feel discouraged by any of such responses hence feel free if you feel like it.

Acknowledgements

I think current and past authors and contributors of BLAS and LAPACK should be knighted and placed in chateaus and mansions with 4 lambos and 2 yachts each. There is no way to measure their contributions to the scientific community. It is in fact make me feel a bit uneasy given how much money is poured into AI and not to this fundamental library.

So I will make some remarks below but this is purely from future-proof perspective and nothing related to past decisions. It is what it is. And in fact I know maintainers of Reference LAPACK already contemplate about similar points. But if I summarize them here, it is going to open up a different discussion which I want to skip completely.

Meat of the argument

The goal is to store LAPACK in a format that is both

Accessible to be compiled in any platform with no significant toolchain problems,
Readable to be ported to any platform with no significant preliminary knowledge in that language.

and currently it is neither with fortran77.

Very brief and opinionated personal summary;

Currently, LAPACK is the fundamental reason why we need fortran compilers or even we need to care about fortran in general if we are not using fortran.

BLAS has been ported to C/Assembly multiple times in various forms and except the reference implementation, nobody uses it in its native form. If BLAS is to NumPy, LAPACK is to SciPy in our case; LAPACK heavily uses BLAS operations but essentially it is just a never-ending stream of for loops and index olympics to do clever stuff with arrays. Thus it has no special property that blocks such porting to whichever language we would be using.

Therefore fundamentally, F77 is the gatekeeper of taking that code and implementing anywhere that has a suitable compiler of any other language. In other words, if the reference implementation was in C, it would have been substantially easier to port it to almost all platforms under the blue sky. And I would speculate that it would have been translated to many languages.

In fact the totality of the numerical community is forced to either find a fortran compiler or use a compiled library but still forced to have a Fortran layer on top to prevent ABI compatibility issues. As every SciPy maintainer by now knows, it is just an endless pain point. Plus, initiatives such as Pyodide are always getting stuck at wrapping fortran code and anytime SciPy rewrites anything from F77 to C they remove a good chunk of unnecessary code that was meant for handling f77.

Besides, Fortran community is dwarfed by C community and also when it comes to low-level optimizations, lot of folks who can perform SSE optimizations cannot do anything about the LAPACK code since it is in Fortran. So it is a loss of capacity and skill. Anything to be done is already being done to BLAS such as BLIS or GotoBLAS/OpenBLAS, since it is a smaller codebase and also simple to follow the functionality goals, but for LAPACK very little has been done, say, OpenBLAS having a custom parallelized ?getrf written in C and so on. Therefore to be able to increase the attention that LAPACK gets from the international community it has to be archived in a proper common and accessible language. And I argue here that F77 is not that. Therefore a port is sorely missed. And similarly that's why modern fortran is not relevant here. Since if a port will take place, we don't need to modernize fortran we need to get it out to a more common language.

So far there have been only two, battle tested, and proven, software languages that has true archival properties. First one is obvious, Fortran77 and the other one is C (I'm talking about multi-decades). Hence we are left with C. Now, it would be unfair to call C as a better fortran alternative when it comes to numerical array computing. Because it absolutely is not. It has less unergonomic parts however it is just tedious to have straightforward code. But a C array is just a pacifier to hide pointer arithmetic. And we are talking about a huge array library. Some might say well what is wrong with A[z + m*y + m*n*x] or A[x][y][z] instead of A[x, y, z] and we don't need to have a debate. I take these syntax as self-evident. You would too if you spent your days in writing array code. And this is the reason C++ is not even considered. It is very difficult to write unabstracted and sound numerical code in C++ without using some strange new features that have questionable future such as C++2017 std::vector all the way to C++ 2023 std::mdspan. It seems like every year another fancy container is discovered in the standard. In C++ 2026 even LAPACK is supposed to be coming in the std::linalg https://www.youtube.com/watch?v=-UXHMlAMXNk. No offense to C++ devs but even the examples look worse than F77 versions. Hence I have a strong opinion that this is not something we should port to.

The most problematic issue for C is Microsoft and its relentless insistence on not conforming to C99 standards. Therefore, it is always some stupid boilerplate code that needs to be included at the top of the project structure and can hurt in the long run. However, it has been done countless times and can stay as a necessary evil. Or can be left for the folks who want to use MSVC. If open source C compilers become widely adopted in Windows ecosystem then probably the best strategy might be is to drop MSVC completely and be done with it. But until then, it would hurt majority of the users. This might sound bad but it is much worse for Windows if Fortran is used already which is why we are here.

An ode to fortran though

Fortran, old or new, regardless of its age, has superior array ergonomics, hands down, no competition. It is fantastic in that regard, so much so that, NumPy borrowed almost all identically. But that's where good things about F77 end for us. Because it is typically 1-indexed and column major memory layout, it does not play well with anything else (except julia and matlab). It is possible to change the indexing in fortran however LAPACK code did not. Hence one of the biggest annoyances of such rewrite will be 2D array index bonanza if C to be selected but the advantages so far outweigh this very important annoyance.

Language designers completely miss the number crunching community and we are always piggybacking on regular langs such as Python's @ operator and ndarray silliness of all C inspired languages and VLAs and countless other ergonomics killing details and so on.

In summary, fortran is a very good language to stick with, if you don't interface with other languages (ctypes are not a thing in this context) or you are not packaging anything to be run on all platforms. Unfortunately, the majority of the contemporary community does and are. Also its future is bound to a very particular audience and to a very small committee skewed towards HPC use cases.

Obligatory Rust and Zig section

Why not Zig or Rust? Rust already has a couple of attempts, a very nice one is https://github.com/sarah-ek/faer-rs and a there are few others via ndalgebra crate. But its impossible looking syntax is killing the starting rationale of this rewrite.

Zig has multidimensional arrays but it is no different than C and also does not come with complex numbers which is in my opinion two very important trains missed. It is now available in std::math but I can't make heads or tails from this issue ziglang/zig#16278. Typically, folks who don't use complex numbers dismiss this as "It's just a struct what about it?!?" and an endless debate ensues.

Why not others langs?

I don't know them. Put it below if you want to make a case. But we might as well use C to get all the middle man out despite all the problems given above.

`#IFNDEF HORROR`, it's here

So what kind of a beast are we talking about? Well, as of 20240813 using cloc tool

-------------------------------------------------------------------------------
Language                     files          blank        comment           code
-------------------------------------------------------------------------------
Fortran 77                    3548           3236         922596         622695
C                             2870          15422          90590         166352
C/C++ Header                    12           3144            393          51340
Fortran 90                      30            364           8594           6040
CMake                           40            520            537           5594
make                            19            445            492           4837
TeX                              1            242            178           1264
YAML                             5             66             44            413
Python                           1             33             23            271
Markdown                         7             53              0            266
XML                              2             17              8            251
-------------------------------------------------------------------------------
SUM:                          6535          23542        1023455         859323
-------------------------------------------------------------------------------

This is a utterly spectacularly annoying problem. However before you close the tab, notice that this is the totality of BLAS/LAPACK tests and other stuff multiple versions etc. If we just look at the SRC directory of LAPACK (no BLAS no nothing else)

-------------------------------------------------------------------------------
Language                     files          blank        comment           code
-------------------------------------------------------------------------------
Fortran 77                    2111           2641         573661         314089
Fortran 90                      18             38           6919           3302
C/C++ Header                     1              2              0           2275
make                             2             55             72            584
CMake                            1             27             39            526
-------------------------------------------------------------------------------
SUM:                          2133           2763         580691         320776
-------------------------------------------------------------------------------

Now it is still at least spectacularly annoying but halved. However note that most LAPACK code is in 4 flavors. That is the same routine for <s, d, c, z> (single, double, single complex, double complex) hence this is actually around 90k SLOC of code or original code. It does not reduce the pain but it definitely reduce the suffering.

Now that is just very annoying but it is precisely how much fortran77 code we had in SciPy too and we are currently passed the halfway point nearing 70%. It took me 16 months to get to here mostly 1/8 FTE or 5 hours/week on average. To be honest that's pretty much doable in my book given my track record in SciPy rewrites.

Also note that not all LAPACK code is essential. Many parts of it has deprecated things or historical warts.

Last points

There is no more important code in the scientific ecosystem than LAPACK. Your blackhole research or cure for cancer can resume in a much much better state if we collectively do this and be done with it because then we can in fact handle this unspoken embarrassing situation.

As I have written in my acknowledgements, we must support and lift the original reference library to the next level.

In case it will avoid pointless debates for the future, I don't like working with C especially after using C3, Zig and bunch of other langs which fix "a different part of C and not the parts I need". I am just forced into it because of the lack of another specialized array language. My unicorn wish is to have a truly "arrays are at the center" language. Think of the ever so wonderful Lua and tables, but instead <insert lang> and arrays.

A few small point to summarize the pages of thoughts I could have written about this topic (not enough time, sorry):

(1) For SciPy, once we get rid of all our own Fortran code, we have no fundamental problem with Reference LAPACK being written in F77. The problem we have with BLAS and LAPACK is an ABI problem - multiple problems actually, from LAPACKE not being well-enough supported to the routines like cdotc with two incompatible ABIs, to ILP64 symbol names. That the reference implementation is in Fortran doesn't matter so much, since we never want to use it at runtime (except maybe for debugging).

(2) I'm honestly not sure if it would help for LAPACK to be rewritten. Whether some vendor or open source effort writes an accelerated implementation is also based on need. The payoff for an individual routine to be written with SIMD/assembly/etc. is much smaller for LAPACK than it is for BLAS, so the chances of that effort being done are also much lower.

A few small point to summarize the pages of thoughts I could have written about this topic (not enough time, sorry):

We should get to a retreat in a cabin somewhere so we can compare pages.

(1) I can also go the other way, the reason why you have ABI problems is actually f77. The very reason for LAPACKE existing is also f77. I can give you another word that hides a saga; LSAME. Or you can think of NaN/INF handling, if you think of all the ABI issues and which language they are originating from, there is a pattern there. The main issue is the absence of a formal fortran ABI. I understand your view since this is a SciPy pain and in your view if the ABI was OK, why bother. However, this is not a SciPy problem and much larger. And f77 has blocked all the possible progress with the exact argument so far spelled out who knows how many times.

(2) This one I don't follow yet. Whatever makes BLIS fast can be applied to intermediate steps in LAPACK too which I think what https://github.com/flame/libflame was supposed to be in the first place. It is just not as tidy as BLAS routines. If you can unroll things in BLAS, you can still unroll in LAPACK too. Also this assumes that there is no room for improvement in LAPACK which is quite not true. Even LU can do much better. That's what Julia/Rust folks are trying to generalize in many different codebases. But the code cannot be modified as is and has to be rewritten. The main problem is that it is just not feasible to decode all that old crust and reorganize the new code and squeeze it back into F77 without causing mayhem. After quite some thousands of line of code translation, I can guarantee you, there are low-hanging fruits left and right everywhere. You can do the same exercise, just look at our translation of minpack or quadpack C code. There are many places I know that it can be 10 lines long and 10 times faster instead of bare looping with insane-level index tracking across nested functions.

The point I wanted to make is that it is possible to do optimizations without knowing the actual algorithm if you are using C if you encounter a for loop. In f77 you need a surgery. Much easier in modern fortran but not in f77. Also it is not "practically" possible to do it in F77 unless you know everything about the algorithm because the algorithm is baked into indices such as if if (mod(j, m).eq.1) GOTO 888. Why 1? Why mod? Because f77 does not have proper data structures so every flow is guided by counting with single letter indices hiding massive opportunities of performance.

We were holding the same idea about our fortran codebase being not feasible to translate, needs unreasonable human hours, they are battle-tested, crazily optimized yada yada. And I am not that blind or that stupid, I saw them everywhere too. We had all the reasons not to do it but now we are where we are. If we do not do this step, regardless of the bugs that will seep through, we are not going to go anywhere and still be discussing ABI issues in the decades to come which is a big waste of time to treat the symptom. To be honest, I am starting to think that everybody is making too much of a deal about these old algorithms. They are not magical code written in Ocaml or Assembly but just flat out, for loops. There are just too many of them.

@ilayn

@ilayn What do you think should be the next step?

A couple more questions:

How good is the test coverage in LAPACK?
What is the coverage of "functions" that we actually need from LAPACK which is already implemented in the rust libraries (just as an example) - it might be simpler to translate from rust to C

I don't have a concrete plan yet to start the work but the following are something that would help understanding the situation;

Opinion of actual maintainers of these libs would be really nice. If they say this is stupid, that's also a data point. I just want to know what the current feeling is. Because this is always thought as impossible task and I tend to agree. I'm just trying to understand how much impossible it is.
If C is the right solution or not.
Should the lwork query be eliminated or not.
Most algorithm is written with column-major in mind. So there will be some changes needed (loop orders to start with).

and many more questions. However, before all this a few examples should be done I think such as ?gesv, ?getrf, ?potrf then by swapping the underlying BLAS vendor such OpenBLAS, BLIS and so on would give some idea of what we are dealing with.

How good is the test coverage in LAPACK?

Not sure, but there is some. They are quite strict (up to a few ULP) hence when there is a test I'd expect it to have good coverage. But this is never guaranteed.

What is the coverage of "functions" that we actually need from LAPACK which is already implemented in the rust libraries

I don't know any direct Rust versions. But many rewrites are providing the same functionality by re-coding (not translating). But I am not sure, I did not really sit down for a comprehensive research.

Note that NumPy used to avoid LAPACK by using the so-called "LAPACK LITE" in a similar fashion but it was just f2c'd LAPACK.

Keeping in mind that I am not a linear algebra person: I tend to think of C as legacy language at this point, and more broadly any language that lacks memory safety by default. So personally I'd go with Rust for any new infrastructure projects; syntax can be learnt, a C ABI can be exposed... My impression is that https://docs.rs/num/latest/num/complex/struct.Complex.html is the standardized Rust-specific API/ABI for complex numbers, it's used by ndarray, nalgebra, faer, and sprs crates. The matrixmultiply crate (inspired by BLIS) appears to use an experimental respresentation that is ABI-compatible with the previous one. (C ABI compat for complexs seems like more of a a pain though, because it sounds like there's no cross-architecture ABI: rust-lang/rfcs#793).

A random thought: https://github.com/immunant/c2rust transpiler is apparently successful enough to be used in some real world projects, e.g. as the starting point for a video codec (https://www.memorysafety.org/blog/porting-c-to-rust-for-av1/). One can imagine the equivalent for Fortran 77.

As I tried to argue above, C is indeed the legacy language however memory safety is not an issue on these decades old code. Hence it is a matter of whether all the necessary details are included in the language. The complex number issue is one, like you point out, but also the 2D array libraries are still competing so it seems to me that there is a very high cost of choosing the wrong one.

One can imagine the equivalent for Fortran 77.

This only exists for F77->Modern Fortran with partial success so far and I don't think there is any help coming anytime soon.

Have you guys seen that there is a Modern Fortran implementation? https://github.com/perazz/fortran-lapack

It is already included on the Fortran stdlib with preprocessing for different types: https://github.com/fortran-lang/stdlib/blob/master/src/stdlib_linalg_lapack.fypp

I do and have discussed this on fortran-lang discourse. It is a reformulation of the modern function declaration with little or no code rewrite (their parsing python script is also available in case interested) hence it is not solving the problems I have mentioned above.

@ilayn

Responding to @ilayn:

I don't have a concrete plan yet to start the work but the following are something that would help understanding the situation;

Opinion of actual maintainers of these libs would be really nice. If they say this is stupid, that's also a data point. I just want to know what the current feeling is. Because this is always thought as impossible task and I tend to agree. I'm just trying to understand how much impossible it is.

If C is the right solution or not.

Should the lwork query be eliminated or not.

Most algorithm is written with column-major in mind. So there will be some changes needed (loop orders to start with).

and many more questions. However, before all this a few examples should be done I think such as ?gesv, ?getrf, ?potrf then by swapping the underlying BLAS vendor such OpenBLAS, BLIS and so on would give some idea of what we are dealing with.

I think the LAPACK maintainers would balk at the idea of rewriting it in its entirety. However, it's already been pointed out earlier in this thread that SciPy only needs a small subset of LAPACK. It's conceivable that that specific subset of LAPACK could be rewritten to fit NumPy+SciPy needs. The two major challenges on that front would be (1) the shear scale of LAPACK's code for eigenvalue problems, and (2) the fact that LAPACK's approach to blocking algorithms is actually very very important for decent speed. So while it wouldn't be impossible for SciPy to implement LAPACK-level functions according to its needs, it would take a ton of effort to do so without dramatic impacts to performance.

The LAPACK maintainers generally feel that C++ is (and/or will be) the successor to Fortran for numerical linear algebra. Various subsets of them have two active efforts in view of this.

TLAPACK implements LAPACK-level algorithms (from scratch!) while templating on the matrix datatype. Its algorithms can run on plain buffers, on Eigen matrices, or on C++/Kokkos mdspans. Development of TLAPACK is quite active, but it's still early days.
LAPACK++ is a portability layer for vendor LAPACK that makes very modest use of features specific to C++. Its API templates on numerical precision (s/c/d/z) and dispatches the vendor function appropriately. It does away with the concept of lwork, handling the workspace query call for the user. It requires that matrices be stored in column-major order just like you'd expect with classic LAPACK. Finally, it has proper exceptions rather than only returning error information in status codes. Development of LAPACK++ has slowed, but that's not much of an issue for CPU-only code, where their coverage is quite thorough--- more thorough than would be needed for NumPy+SciPy.

I agree with the sentiment but the issue is much more convoluted than yeah let's switch to C++ for variety of reasons. TLAPACK in particular this is pretty good summary of the situation Reference-LAPACK/lapack#988 (comment) For LAPACK++ I would only mention that templating is not equal to generic programming. Using C++ for generating 4 flavors of the same code is not meant to be a modernization but a mere compactification with a lot of gotchas and really really bad language practices.

As I tried to argue above, just because it exists should not make C++ a valid candidate. We have been dealing with it for decades now and just because Rust exists seriously 99% of the Cpp conferences are now with talks that mention oh we can make C++ memory safe. I am not saying that memory safety is important in numerical computing (it probably is not) but people have a concrete example to point to and show the benefits of a proper memory-safety benefits in driver code or this or that. And this makes a difference a lot. Because C++ fans finally cannot just say the same things that they were mentioning for decades now. Rust actually makes a difference where it matters. Like I mentioned above this is not about language wars but finally settling for a proper numerical language that does not feel like piggybacking. When Rust exists we don't need to talk about C++ (or C to be perfectly honest).

Just because a better language does not exist yet, should not invalidate, in numerical computing in my opinion, argument's validity. Frankly the only other language that insists on going the wrong direction is the Fortran folks in this regard insisting that we have a better language at home we don't need a new one.

In particular, both committees have zero respect for any kind of improvement in ergonomics or removing legacy behavior and so on. While it might seem like this is a SciPy specific subset, it is in fact the entire LAPACK that I have in mind. Not because it is easy but we suffered sufficiently enough with oh what can you do mindset. This is not magical code, and the reason I know it isn't is because I translated already a lot of code that was deemed magical. But it is just bunch of annoying for loops. Original LAPACK maintainers might indeed balk at it but this is not a blocker of any kind. This code is too important to be left alone as is with bandaids or templates on top of it.

@ilayn

Thanks for your quick and detailed response, @ilayn! I have some follow up points from your message.

For LAPACK++ I would only mention that templating is not equal to generic programming. Using C++ for generating 4 flavors of the same code is not meant to be a modernization but a mere compactification with a lot of gotchas and really really bad language practices.

I get that templating isn't equal to generic programming. However, I think you're selling LAPACK++'s templating short by saying that it's a "mere compactification." I would like to hear more about the "with a lot of gotchas and really really bad language practices," if you care to elaborate. I'm always looking to learn.

While it might seem like this is a SciPy specific subset, it is in fact the entire LAPACK that I have in mind. Not because it is easy but we suffered sufficiently enough with oh what can you do mindset. This is not magical code, and the reason I know it isn't is because I translated already a lot of code that was deemed magical.

Translating all of LAPACK would be silly. For example, GEQRP is utterly obsolete and is completely superseded by GEQP3, and even that is agonizingly slow in the reference implementation. (Think: slower than an SVD.) At the same time, I appreciate this ambition to rewrite more of LAPACK than just the LAPACK functions used in SciPy.

I think that doing a rewrite that preserved the column-major ordering would be relatively straightforward, although still a huge effort. The problem with switching to row-major is that with LAPACK-level operations you can't easily switch from row-major to column-major just by fiddling with transposition flags like we do in GEMM, GEMV, etc... If you try to do Householder QR on a row-major matrix, you're in for a massive performance penalty from cache effects. I know as much because I've used the trick to reduce Householder QR on row-major data to Householder LQ on column-major data (using LAPACK's GELQF), and that's practically as slow as GEQP3! To my knowledge, fixing this isn't a matter of changing the order of some loops.

I would like to hear more about the "with a lot of gotchas and really really bad language practices," if you care to elaborate. I'm always looking to learn.

What I mean is that LAPACK++ is for C++ workflows. If you are working on C++, that's great, you do you (you = whoever is reading). But if you are not then you deal with C++ unwillingly. I don't mean to start bashing C++ since it is not relevant but I think a lot of people had done the bashing already sufficiently well both within and outside of C++ circles. I don't have a beef with it; I just don't think it is the right tool for the job. But like I said, this is a dead horse we don't need to keep beating. I have a control embedded background and in that tiny sliver of a scope, I don't it is even productive. But, say, for gamedevs, somehow C++ is a must so something is done right in that scope apparently, I just don't know what it is.

What is relevant here though is that it is not an archival language for storing algorithms as it continuously getting "modern", without really breaking or fixing anything. Hence here I am questioning both its archival quality and moreover it's interoperability. I hope you agree that in that limited scope the current king is C and C alone. And I repeat this very often I'm not a big fan of C either but it is a simple enough language with relatively sufficient tooling and quite decent support in all platforms. Moreover, almost all languages have functioning C API. Thus my fundamental thesis here is that we need to get LAPACK to a point where it does not need extra tooling that is very scarce and only supports some of the platforms.

Translating all of LAPACK would be silly.

Fully agree. What I could not articulate there is all = full non-deprecated suite. And not something like Eigen that only offers a limited subset not necessarily compatible with LAPACK features.

I think that doing a rewrite that preserved the column-major ordering would be relatively straightforward, although still a huge effort.

I agree. But 80-20% also applies.

If you try to do Householder QR on a row-major matrix, you're in for a massive performance penalty from cache effects.

This is a very good point. However I can also bring certain functions that boost performance for row-major layout (think of row swaps or ??trf workloads. The crux of the matter when it comes to memory-layout favoring algorithms is possibly to work on hardware accelerated cache-oblivious transpose algorithms. Or in the case of BLAS-like behavior, microkernels that can process different layout reads. I think BLIS is doing a great job at that. Then we can keep the algorithms intact. But again these require nontrivial amount of design and thinking. So I totally agree on that however in its current F77 form, nobody will even bother working on these problems due to the reasons I mentioned above. Hence a mechanical translation would not immediately return the investment but to enable these advanced optimization topics, it is a necessary evil in my opinion.

Note that a lot of optimization opportunities already being discovered by many others. The faer-rs project is doing really nice stuff without a hand-optimized special BLAS library https://www.youtube.com/watch?v=LLG3AxRyfTw I think this is because they don't need to operate within the classical shortcomings. Hence there is now prior art that we have quite a good probability that this type of work will do some good. Rust being the language here is not relevant; it is the fact that they, a tiny number of people, can hit these performances is the impressive part.

In the other side of things, I think you can take any C code these days and convert it to any other language due to its simplicity and coding-agent assistance availability. Hence if the reference is a simple straightforward implementation all platforms can take the reference or machine-translate. This is not possible if the reference lives in an obscure flavor of a not so common language.

I try my best to keep the discussion away from the language flames but all I want to emphasize is that "if this is not the correct language then what?" answer is probably not C++ in my opinion. Other than that there are lots and lots of smart folks working on these things so I hope I don't disrespect their efforts with these comments.

Now I'll reply your email before the plane takes off :)

Those are all great points. I certainly agree that C is king from an archival and portability perspective. I also agree that there's more potential for interesting optimizations/innovations if a project has a lower barrier to new contributors.

Going back to your questions in this comment, I think it would be good to keep the workspace query. If you try to get rid of it then you have to deal with dynamic memory allocation and you make it harder to recycle workspace in complicated algorithms. I'm sure formal code generation tools and/or LLMs could make wrapper functions that hide the workspace query.

What are your thoughts on how a C implementation would handle separate implementations for (s,c,d,z) scalar types, and the desire for other precisions? Maybe some light macro magic could get the job done. My understanding of LAPACK development (which could be wrong) is that they only write the z implementation from scratch and then use code generation to get the s,c,d implementations.

If you try to get rid of it then you have to deal with dynamic memory allocation and you make it harder to recycle workspace in complicated algorithms.

True. Workspace queries are typically dependent on a system-dependent parameter NB and once it is known we can compute the exact amount of each lwork number algorithmically. The main issue is that NB parameter detection. Either some heuristics like ATLAS or just hard-coding a known number for different architectures can be done. Then instead of lwork we can ask for a specific size of array like 3n*m + nb*n**2 or so.

What are your thoughts on how a C implementation would handle separate implementations for (s,c,d,z) scalar types, and the desire for other precisions?

No idea yet. I wish we don't do this in C but I don't see anything else that is portable enough. Some folks from Jane Street, a few years ago, pushed me into OCaml similar to what FFTW folks do; "Write it on OCaml and transpile to C". That's also how I found out that I am not alone, see item 2.10 here I don't know any OCaml and I think it would limit the contributor number to fingers. But maybe for the initial brunt of the work it might be the right way to do. Then once enough momentum is there, it can continue with C if there is nothing else than C.

LAPACK folks' strategy is also how I translate SciPy F77 code too; write one version and run an involved regex scheme on the file to have the others.

Here is where I put my money where my mouth is (such a weird idiom).

https://github.com/ilayn/semicolon-lapack

It is in very alpha stage. The mechanical translations of four flavors are pretty much done (99ish%). It runs but there is some way to stability. I still need to port all the tests, hence it's a work in progress.

Probably the most controversial one ; it is all 0-indexed. Everything; from arguments to internals to ipiv-like outputs.

Performance is good and better in some cases, if you turn on the machine specific optimizations even better (only tested very few things and only against OpenBLAS). I don't know why maybe only on my own machine or a fluke.

It can detect the integer size at compile time and then be linked to the following (tested with all of them)

OpenBLAS LP64 (conda, libopenblas.so) "LP64"
OpenBLAS ILP64 (local build, INTERFACE64=1) "ILP64"
MKL LP64 (conda, libmkl_intel_lp64) "LP64"
MKL ILP64 (conda, libmkl_intel_ilp64) "ILP64"
Reference BLAS LP64 (local build, libcblas.so) "LP64"
Reference BLAS ILP64 (local build, BUILD_INDEX64, libcblas64.so) "ILP64"
BLIS LP64 (local build, --enable-cblas) "LP64"
BLIS ILP64 (local build, --int-size=64 --blas-int-size=64) "ILP64"

The integer agnostic-version is not pushed yet (hopefully this week) but I did not get any issues. Controls its own CBLAS header so no more lower case, upper case, CBLAS_FUNC(...) etc.

Generated an entirely different documentation. It takes a while but it looks normal and not doxygen.

I don't know, I did a lot of things; there are many conveniences entirely for my own taste. Maybe good maybe not. Have a go at it if you want. All feedback, criticism, insults etc. are very welcome.

Resulting C code (I wasn't super off)

━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 Language              Files        Lines         Code     Comments       Blanks
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 C                      1933       451241       290056       116939        44246
 C Header                  6         2383         2218          104           61
 Meson                     1         2052         1956           87            9
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 Total                  1940       455676       294230       117130        44316
━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

@ilayn

Very interesting, thanks for sharing @ilayn!

A couple of questions that came up after a quick read of this comment and skimming the repo:

Are ILP64 libraries using the now-standardized _64 symbol suffix handled?
Are "dual builds" that contain both LP64 and ILP64 symbols in the same library handled, and is that an either/or handling, or is it possible for semicolon-lapack to provde both LP64 and ILP64 routines in a single shared library?
It sounds like one still needs to make a build-time choice, correct? Did you think about or play with ideas to make the BLAS library selectable at runtime? Or is that out of scope for this project, and could be built separately if desired?

No, I did not make any choices yet; it's pretty much up to us.

But because we control the compilation, we can generate, rename symbols and so on. One thing I am trying to figure out is the two-pass compilation (that has both symbols), essentially a fat binary. Then runtime switching is not an issue. However, then we need to link to two CBLAS (LP and ILP) separately. It's not impossible but I'm not sure if it is better than having two separate binaries.

I am mostly pushing to finish the test porting and correctness to reach to an acceptable level, so I did not spend too much time on the toolchain. Plus, not really an expert on multi-platform architecture. Let me know if I am doing something wrong, it's not really meant to be a personal project, it's better if we shape it towards a generic use.

However, then we need to link to two CBLAS (LP and ILP) separately. It's not impossible but I'm not sure if it is better than having two separate binaries.

Yeah I agree, that's not better. If it's done at all, it should only be done if the shared library being linked against contains both sets of symbols (i.e., a "BLAS fat binary") in my opinion.

I am mostly pushing to finish the test porting and correctness to reach to an acceptable level, so I did not spend too much time on the toolchain.

Makes perfect sense, first things first. I can see value in a pure C11 translation of LAPACK; my mind was already starting to imagine possibilities / use cases:)

I'll try to play with what you have now when I get some time.

Performance is good and better in some cases, if you turn on the machine specific optimizations even better (only tested very few things and only against OpenBLAS). I don't know why maybe only on my own machine or a fluke.

Can you add some benchmarking infrastructure to the repo? I'd love to try this out.

I'm not sure if this qualifies as a benchmark but here is an attempt.

https://github.com/ilayn/semicolon-lapack/tree/main/benchmark

I will try to shave off some time from work, to try it on a clean env.

@mmelnich

I had some more time to work on this;

all LAPACK test suite for double and single (d, s) precision is ported, complex is being ported
ILP64 builds are now possible with custom symbol mangling mechanism in place. So you can name your ILP64 symbols whatever you want (more here https://ilayn.github.io/semicolon-lapack/building.html#symbol-naming)
Correctness-wise at least for the double precision the confidence is increased, complex versions still in the limbo. However the library as a binary can be experimented with. I'll try to finish the complex tests this week.
The coverage of LAPACK testing suite is high but not comprehensive. Hence some smoke tests are added to hit certain workflows.
BLIS CBLAS interface turned out to be incomplete, which is a pity, and OpenBLAS is going to get cblas_?gemmtr in the latest version 0.3.32 (so I have some hacks to accomodate for it).
I also placed a poison-header for cblas.h which means nothing explicitly depends on any vendor CBLAS header.
OpenBLAS/Reference BLAS/MKL both LP64/ILP64 are passing the tests.
Thanks to @mmelnich we have a better benchmarking mechanism in place.

I guess from now on, once the complex tests are ported, it is just correctness fixes and Reference BLAS/LAPACK repo tracking since the majority of the infra is in place.

I also finally got the numbers I wanted to see which is linking to existing CBLAS interfaces;

Configuration	Input `.a`	Output `.so`	Runtime deps	Multi-arch?
Dynamic linking (default)	—	6.9 MB	libopenblas.so.0	Yes (vendor)
CBLAS-only, single target	12 MB	11 MB	libm, libc only	No
CBLAS-only, `DYNAMIC_ARCH=1`	47 MB	33 MB	libm, libc only	Yes
Full vendored OpenBLAS static	70 MB	33 MB	libm, libc only	Yes

The first one compiles as a shared library and works with the system/user-supplied CBLAS binary. Input here means the CBLAS vendor size. Single target is what user compiles for their own architecture (mine is Zen/Cooperlake).

The first two shows when the user brings their own binary we are around 7MBs and jumps to 11MBs if we subsume the vendor objects as a static library.

If we use the dynamic architecture fat binaries, we inherit all, as expected, since CBLAS is highly architecture optimized.

The last two rows shows the cost of Fortran interface which has no effect on this library.

@ilayn

@ilayn I'm very interested in your C port of LAPACK! Just bringing to your attention the proposed LAPACK (and BLAS) registry that would capture semantic information about the routines and their parameters, as discussed on Reference-LAPACK/lapack#1349, which I mention since it seems like it might help your effort.

Navigation Menu

Repository navigation

Rewrite LAPACK so it's not in F77 #7

Acknowledgements

Meat of the argument

Very brief and opinionated personal summary;

An ode to fortran though

Obligatory Rust and Zig section

Why not others langs?

`#IFNDEF HORROR`, it's here

Last points

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Rewrite LAPACK so it's not in F77 #7

Description

Activity

ilayn commented on Aug 13, 2024

Acknowledgements

Meat of the argument

Very brief and opinionated personal summary;

An ode to fortran though

Obligatory Rust and Zig section

Why not others langs?

#IFNDEF HORROR, it's here

Last points

ilayn commented on Aug 13, 2024

rgommers commented on Aug 16, 2024

ilayn commented on Aug 16, 2024

alonme commented on Sep 23, 2024

ilayn commented on Sep 23, 2024

itamarst commented on Sep 23, 2024

ilayn commented on Sep 24, 2024

fedebenelli commented on Sep 24, 2024

ilayn commented on Sep 24, 2024

rileyjmurray commented on Jul 10, 2025

ilayn commented on Jul 10, 2025

rileyjmurray commented on Jul 11, 2025

ilayn commented on Jul 11, 2025

rileyjmurray commented on Jul 12, 2025

ilayn commented on Jul 15, 2025

ilayn commented on Feb 16, 2026

rgommers commented on Feb 16, 2026

ilayn commented on Feb 16, 2026

rgommers commented on Feb 16, 2026

rileyjmurray commented on Feb 16, 2026

ilayn commented on Feb 17, 2026

ilayn commented on Feb 23, 2026

mohawk2 commented on Oct 25, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

`#IFNDEF HORROR`, it's here