FortranLegacyRunbook Editorial·2026-08-30·9 min read

The Reference Implementation: How 47-Year-Old Fortran Still Anchors NumPy, MATLAB, and R

BLAS shipped as a Fortran library in 1979. LAPACK moved to Fortran 90 in 2008. Every major numerical computing platform in daily use still calls into that lineage for its actual linear algebra — and the 2021 Turing Award said so explicitly.

Most conversation about Fortran, including a lot of the reference material on this site, is understandably about maintaining it directly: IOSTAT codes, compiler-specific runtime errors, the gap between what a fixed-size array catches and what an allocatable one lets slide. That's real, and it matters to the people staring at a segfault in a thirty-year-old solver. But there's a second, far larger Fortran footprint that almost nobody maintaining that solver ever stops to think about: the reference numerical libraries sitting underneath nearly every modern data-science and scientific-computing tool in daily use, written in Fortran, and never replaced.

This isn't a story about a niche corner of the language's install base either. It's the layer that runs every time someone calls `.dot()` on a NumPy array, fits a linear model in R, solves a system of equations in MATLAB, or trains a model on a CPU fallback path. Almost none of the people doing that know, or need to know, that the actual arithmetic underneath their call is still, at its core, Fortran. That invisibility is itself worth noticing: legacy Fortran usually gets discussed as a liability someone has to manage. This particular Fortran is different. It's infrastructure so reliable that an entire generation of scientific software was built directly on top of it without a second thought.

The library nobody had to choose to depend on

BLAS, the Basic Linear Algebra Subprograms, is the standard low-level specification for vector and matrix arithmetic: dot products, matrix-vector multiplication, matrix-matrix multiplication, the actual arithmetic kernels almost every numerical algorithm eventually bottoms out in. It didn't start as an abstract specification. It started as a Fortran library.

1979The original BLAS was published as "Basic Linear Algebra Subprograms for FORTRAN usage." The reference implementation is still a public-domain Fortran library today, and it's what vendor-optimized versions are built to match. (Wikipedia: BLAS)

The specification grew in stages that map directly onto how much work each level of operation does: Level 1, from 1979, covered vector-vector operations. Level 2, added between 1984 and 1988, covered matrix-vector operations. Level 3, finalized around 1990, covered matrix-matrix operations, the most computationally expensive tier and the one modern hardware vendors spend the most engineering effort optimizing. Sparse-matrix extensions followed in 2002. Every one of those additions preserved the same calling convention the original 1979 Fortran routines defined, which is why "most libraries that offer linear algebra routines conform to the BLAS interface" remains true nearly half a century later: Intel's MKL, OpenBLAS, Apple's Accelerate framework, and AMD's AOCL are all different, independently engineered implementations of the exact same interface a Fortran library defined when Jimmy Carter was in office.

That's worth sitting with for a second, because it's an unusual kind of longevity. Hardware vendors compete fiercely on how fast their BLAS implementation runs on their own chips, hand-tuning assembly, exploiting specific cache sizes, adding vector instructions the original authors never imagined. What they don't compete on, and haven't in nearly fifty years, is the shape of the interface itself: the argument order, the naming conventions, the way a matrix gets described to the routine that operates on it. That interface is still the one a Fortran calling convention defined, because rewriting it would break every program built on top of it, and by the 1990s that list of programs was already too long to touch.

LAPACK: designed to replace the Fortran of the 1970s, in Fortran

LAPACK, the Linear Algebra Package, was built as the direct successor to two earlier packages: LINPACK's linear-equation and least-squares routines, and EISPACK's eigenvalue routines. Both predecessors were tuned for the vector supercomputers of the 1970s and 1980s and ran poorly on the cache-based architectures that came after. LAPACK was designed specifically to exploit those caches, and runs orders of magnitude faster than LINPACK on the same modern hardware. It has never left Fortran; it simply moved to a newer dialect of it.

Fortran 90 since version 3.2, in 2008LAPACK was originally written in FORTRAN 77 and moved to Fortran 90 in version 3.2, released in 2008. R, MATLAB, and SciPy are all documented as being built on top of it. (Wikipedia: LAPACK)

The one real concession LAPACK made to a non-Fortran world came in the form of LAPACKE, a standardized C interface added as part of LAPACK starting with version 3.4.0. It's worth being precise about what that concession actually was: not a rewrite of the underlying math, but a wrapper around it. The routines other languages call into are still the same Fortran routines; LAPACKE just gave C programs a cleaner way to reach them.

What's actually calling into it, right now

It would be easy to assume this is a historical footnote that only matters to people who still write raw Fortran. It isn't. NumPy, the array library underneath most of the Python data-science and machine-learning stack, says exactly what it depends on for the actual math, in its own documentation.

The NumPy linear algebra functions rely on BLAS and LAPACK to provide efficient low level implementations of standard linear algebra algorithms.

NumPy documentation, Linear algebra (numpy.linalg)

NumPy isn't the exception; it's the norm. The same BLAS interface is documented as underlying SciPy, GNU Octave, Mathematica, MATLAB, R, and Julia, alongside LAPACK, LINPACK, and Armadillo. A data scientist running a regression in R, an engineer solving a system of equations in MATLAB, and a machine-learning researcher calling `numpy.linalg.solve` are, underneath a few layers of abstraction, all reaching the same family of Fortran-defined routines. On the deep-learning side specifically, PyTorch and TensorFlow's CPU tensor operations for linear algebra commonly route through the same BLAS implementations, OpenBLAS or Intel MKL among them, that everything else in this chain calls into; the GPU path runs through NVIDIA's cuBLAS instead, a separate implementation that still honors the same calling convention. A stack that looks, from the outside, like it's built entirely on Python and C++ is, at the layer that actually does the arithmetic, still leaning on an interface a Fortran library defined in 1979.

The Turing Award committee said this part explicitly

Jack Dongarra wrote or co-wrote LINPACK, BLAS, LAPACK, and ScaLAPACK, and in 2021 the Association for Computing Machinery gave him its highest honor for exactly that body of work.

for pioneering contributions to numerical algorithms and libraries that enabled high performance computational software to keep pace with exponential hardware improvements for over four decades.

2021 ACM A.M. Turing Award citation, awarded to Jack Dongarra

That citation isn't nostalgia for old software; it's computing's most prestigious award naming, specifically, decades-old Fortran numerical libraries as work that let high-performance computing keep up with hardware at all. Dongarra's software is documented as incorporated into MATLAB, Maple, Wolfram Mathematica, GNU Octave, R, and SciPy directly. When people talk about which programming language "won," they're usually counting new projects, job postings, or forum chatter. None of those metrics capture a routine that gets called billions of times a day from inside software nobody thinks of as Fortran at all.

The Turing Award isn't a one-off either. Dongarra has been elected a Fellow of the American Association for the Advancement of Science, the Association for Computing Machinery, the Society for Industrial and Applied Mathematics, and the Institute of Electrical and Electronics Engineers, a foreign member of the Royal Society, elected to the U.S. National Academy of Engineering in 2001 specifically for contributions to numerical software, and elected to the U.S. National Academy of Sciences in 2023. That's four decades of the field's own institutions, repeatedly, independently, choosing to formally recognize work built on top of a Fortran calling convention from 1979. It's difficult to read that pattern as anyone treating the underlying code as obsolete.

The other real bridge: F2PY

There's also a direct, named tool that formalizes exactly the kind of interoperability a legacy Fortran shop would actually use, and it's been part of NumPy's own toolchain for years: F2PY, the Fortran-to-Python interface generator, originally created by Pearu Peterson. It lets Python code call Fortran 77, 90, and 95 subroutines and module procedures directly, and reach into Fortran 77 COMMON blocks and Fortran 90/95 module data, including allocatable arrays, without hand-writing a C binding layer in between. That's not a hypothetical migration path invented for this article. It's a documented, currently maintained piece of NumPy's own build tooling, which means the exact pattern of "wrap the existing Fortran instead of rewriting it" isn't a theory about what a modernization vendor might suggest, it's the literal mechanism the most widely used numerical Python library ships with by default.

Why this matters if your Fortran has nothing to do with any of this

If you maintain an in-house Fortran codebase, a structural solver, a reservoir simulator, a legacy actuarial model, there's a real chance it already links against BLAS or LAPACK directly, since both predate most application code that calls them. That means two practical things. First, when correctness of a legacy solver is in question, the reference BLAS and LAPACK implementations are the actual industry ground truth vendors validate their own optimized libraries against, not an arbitrary comparison point; testing a suspect routine's output against the reference implementation is a legitimate, well-precedented way to isolate whether a bug is in your code or in how it's being called. Second, if a modernization project involves exposing Fortran numerics to Python or another language, that path, wrapping Fortran routines for a higher-level caller, is not a novel or risky idea invented for this one project. It's the exact architecture NumPy itself is built on, at global scale, for going on three decades.

That has a direct bearing on how a modernization decision should actually get scoped. "Rewrite the Fortran" and "wrap the Fortran" get treated as opposite ends of a spectrum, but the BLAS and LAPACK lineage shows a third option that's been running in production the whole time: leave the numerical core exactly as it is, because it's already correct and already fast, and only replace the layers around it that actually need modernizing, the file I/O, the interface, the reporting. A team staring down a full rewrite of a Fortran solver is often better served asking which specific pieces genuinely need to change, rather than assuming the numerics do just because the era they were written in did.

Where this leaves you

  • BLAS shipped as a Fortran library in 1979 and its reference implementation still is one; the interface it defined is what Intel MKL, OpenBLAS, Apple Accelerate, and AMD AOCL all independently implement today.
  • LAPACK, the direct successor to LINPACK and EISPACK, moved from FORTRAN 77 to Fortran 90 in 2008 and has stayed there; its one concession to other languages, the LAPACKE C interface, wraps the same Fortran routines rather than replacing them.
  • NumPy says directly, in its own documentation, that it relies on BLAS and LAPACK for its actual linear algebra — and SciPy, MATLAB, R, GNU Octave, Julia, and Mathematica are documented as built on the same foundation.
  • The 2021 ACM Turing Award went to Jack Dongarra explicitly for this body of numerical software, the field's clearest possible statement that this specific decades-old Fortran work remains foundational rather than historical.
  • A legacy Fortran codebase linking BLAS or LAPACK directly isn't isolated from the modern data-science stack; it shares a dependency with it, which makes both correctness-testing against the reference implementation and interop with Python well-trodden paths rather than experiments.
See the Fortran modernization guideSee the Fortran reading list