BLAS, LAPACK and Numerical Libraries

BLAS (Basic Linear Algebra Subprograms) is a specification of low-level routines, originally written in Fortran 77, in three levels. LAPACK (Linear Algebra PACKage) is a library of higher-level algorithms written in terms of BLAS calls, so a fast BLAS makes all of LAPACK fast. Both are interfaces: the netlib distributions are the reference implementations, written for clarity and not for speed, and every serious platform supplies its own optimized implementation of the same interface. A program written against the interface can therefore be moved from the reference library to OpenBLAS, oneMKL or a vendor library by relinking, with no change to the source.

The three BLAS levels differ in what they operate on and, as a consequence, in how much arithmetic they do for each number read from memory. That ratio is what an optimized library exploits.

Level 3 is the one that matters for performance: an optimized dgemm reaches a large fraction of the peak speed of the processor, and LAPACK is organised so that its blocked algorithms spend most of their time in level 3 calls.

dgesv solves a general square system A x = b by LU factorization with partial pivoting. It overwrites a with the factors, b with the solution, and ipiv receives the row interchanges. The example is one implicit heat step on a rod of five points with a hot spot in the middle: with r = alpha dt / dx^2, the system is (I - r L) u_new = u_old, where L is the second-difference matrix and the two end values are held at zero.

Output (identical with GNU Fortran 16.2 and LLVM Flang 23):

The solution is the expected smoothed hot spot, and the residual maxval(abs(a0*u - u0)) is at rounding level: always compute it, because info = 0 only says that the factorization completed, not that the answer is accurate.

Source: LAPACK Users' Guide: Naming Scheme, LAPACK reference documentation

LAPACK reports problems through its info argument and, for illegal arguments, through the routine xerbla. The convention is the same for every routine, and every call should check it:

The info > 0 case is the one that depends on the numbers. Two matrices, one of them singular only on paper:

Output (identical with both compilers):

Source: LAPACK reference documentation (dgesv)

The first matrix has rank 2, and info is still zero: rounding during the elimination leaves a last pivot of order 10^-16 and not exactly zero. The routine cannot know that this number is "really" zero. info > 0 is a report about exact zeros only; for ill-conditioned systems the check to rely on is the condition number (LAPACK offers dgecon, and the expert driver dgesvx returns an estimate) or the residual, as in the previous section.

An illegal argument is a different kind of failure. Here lda = 2 is smaller than the matrix order 3:

With OpenBLAS the library prints a message from xerbla and returns, with info set to the negated argument number:

The message is written to standard error and the program carries on. LLVM Flang 23 prints the same two lines, in the other order, because standard output is buffered when it is a pipe. The behaviour of xerbla belongs to the library. The reference xerbla (in the BLAS directory of the netlib sources) prints the same message and then executes a plain stop. Linking the same program statically against the netlib liblapack.a and libblas.a of the container printed only the line ** On entry to DGESV parameter number 4 had an illegal value and never reached the print after the call; the process exit status was 0, because a plain stop reports success. A program that has to survive a bad argument must check info, and a program that has to detect it everywhere must run its tests against the reference library as well as the optimized one.

Source: LAPACK Users' Guide: Naming Scheme, LAPACK reference documentation (dgesv, dgetrf, dpotrf)

A real symmetric matrix A is positive definite if x^T A x > 0 for every nonzero x. Such a matrix has a unique factorization

where L is lower triangular with a positive diagonal, called the Cholesky factor. It is the fastest and most stable way of solving systems with positive definite matrices (about half the work of an LU factorization, and no pivoting), and the standard test for positive definiteness is to try it and see whether it fails. Positive definite matrices are everywhere in simulation: the implicit heat step above, covariance matrices, stiffness matrices and normal equations are all of this kind. Column by column, for j = 1, …​, n:

If the quantity under the square root is not positive, A is not positive definite. The program below does three things: it codes the formulas with whole-array dot_product calls, it calls LAPACK’s dpotrf on the same matrix and compares the factors, and it uses dpotrs to solve with the factor. The matrix is the implicit rod step with four points.

Output (identical with GNU Fortran 16.2 and LLVM Flang 23):

Source: LAPACK reference documentation (dpotrf, dpotrs), J3/23-007r1: 16.9.73 DOT_PRODUCT

What to take from it:

(The output is identical with LLVM Flang 23, and matches the dgesv solution of the previous section.) Choose the routine for the structure of the matrix: general (ge), symmetric positive definite (po), tridiagonal (gt, pt) or banded (gb, pb). The structure-specific routines are faster and use less memory.

Source: LAPACK Users' Guide: Naming Scheme, LAPACK reference documentation (dpotrf, dpotrs, dptsv)

Several LAPACK routines need a scratch array work of a size that depends on the problem and on the block size the library considers best. The routine tells the size itself: call it with lwork = -1 and it computes the optimal length, stores it in work(1) and does nothing else. Allocate that much and call again.

The test matrix is the one-dimensional second-difference matrix, whose eigenvalues are 2 - 2 cos(k pi / 5). Output (identical with both compilers):

The documented minimum for dsyev is lwork >= max(1, 3*n - 1), which is 11 for n = 4. The optimal value is much larger on OpenBLAS, 136; using the minimum would work and be slower. The query is cheap, and a function in the program that owns the work array can do it once and keep the array. Using an allocatable work array is the modern replacement for the fixed-size scratch arrays and the lwork arithmetic of old code.

Source: LAPACK reference documentation (dsyev)

Linking is where most of the time is lost with numerical libraries, because nothing in the source says which BLAS the program gets. The options, all of which provide the same interface:

Two practical points. First, on a system with the Debian alternatives mechanism, the same reference-style link line -lblas -llapack may end up with OpenBLAS at run time: in the container used for this page libblas.so.3 is a link managed by update-alternatives and points to the OpenBLAS build. That is a feature of the packaging (a program can be switched to another implementation without relinking), and a surprise when timing the reference library. To measure the netlib code the page links libblas.a statically. Second, link order matters on Unix linkers: a library must appear before the libraries it depends on, so -llapack -lblas, never -lblas -llapack.

The language’s own matmul (16.9.134) and dot_product (16.9.73) intrinsics compute the same mathematical products as dgemm and ddot, and they are always available. The standard defines the result and nothing about how it is computed, so how fast matmul is, and whether it uses a BLAS, is a quality-of-implementation matter. The differences that matter:

Reference BLAS, OpenBLAS and matmul were timed on one machine (an Apple-silicon Mac running the aarch64 container, 10 cores available to it, GNU Fortran 16.2 with -O2), multiplying two 1000 by 1000 double precision matrices with this program:

These numbers come from single runs and vary: the multithreaded OpenBLAS dgemm took between 0.014 s and 0.051 s in repeated runs, mostly because of thread start-up, and the matmul figures between 0.049 s and 0.088 s. The ratios are what to read. The reference dgemm was about six times slower than GNU Fortran’s own matmul (the netlib code is not meant for speed); the optimized dgemm was faster than matmul by about 1.2 times with one thread and between 1.7 and about 4 times with all the cores, not by a factor of a hundred; matmul in LLVM Flang 23 was several times slower than in GNU Fortran. These are measurements of one machine and one matrix size. A different processor, size or thread count changes them, so measure on the machine that matters before deciding, as described in Performance and Memory Layout.

GNU Fortran can route matmul to a BLAS without a source change. -fexternal-blas makes the compiler generate calls to BLAS for some matrix operations like matmul when the matrices are larger than a limit, -fblas-matmul-limit=n sets that limit (it compares the size, or the geometric mean of the dimensions if the matrices are not square, with n), and -fexternal-blas64 is the variant for a BLAS with 64-bit integers (it also requires -ffrontend-optimize); the BLAS library has to be named at link time. The table shows the effect: with OpenBLAS matmul became about 2.6 times faster and came close to the dgemm time, and with the reference BLAS it became slower, since the reference dgemm is slower than the compiler’s own code. The option exists to get the benefit of an optimized BLAS in existing code and is only sensible when one is linked. In LLVM Flang 23.1.4 neither flang --help nor --help-hidden lists an option with BLAS in its description, so there is nothing equivalent to check.

Source: GNU Fortran manual: Options for code generation conventions, J3/23-007r1: 16.9.73 DOT_PRODUCT and 16.9.134 MATMUL

Writing interfaces by hand is tedious, so the fortran-lang standard library provides them. stdlib_linalg is a high-level, kind-generic interface to linear algebra built on BLAS and LAPACK: solve, lstsq, det, inv, chol, eigh, eig, svd, qr, schur, norms and helpers such as diag, eye, trace and is_symmetric. The module and its status are documented in the stdlib specification; individual procedures are marked there as stable or experimental, so check the status before depending on one. The linalg page also documents stdlib_linalg_blas and stdlib_linalg_lapack, a Modern Fortran translation of reference LAPACK 3.10.1 (free format, implicit none, explicit intent) with kind-generic names (gemm alongside the original stdlib_dgemm). The documentation says that optimized libraries should be preferred when available and that defining the macros STDLIB_EXTERNAL_BLAS and STDLIB_EXTERNAL_LAPACK when building stdlib wraps the external libraries instead of the built-in translation.

Used from fpm, with the stdlib = "*" registry dependency and an external BLAS to link; the project below is the manifest:

Source: fpm manifest reference: Link external libraries

On 2026-10-19 the registry dependency resolved to stdlib 0.8.1 (the stdlib-fpm branch with pre-processed sources). The program solves the implicit heat step, factorizes it, computes eigenvalues and checks the singular case with the error type linalg_state_type:

Built with GNU Fortran 16.2 through fpm (fpm’s default flags, not the strict -std=f2023 line):

Source: Fortran stdlib: linalg specification

The results agree with the hand-written dgesv/dpotrf examples above (the solution x = 16.9492 76.2712 ..., the same Cholesky residual at rounding level), and the product of the eigenvalues equals the determinant. Things to know about the interface:

stdlib was built and run with GNU Fortran only; the same examples were not built with Flang through fpm. The package itself is described in Fortran stdlib.

Source: Fortran stdlib: linalg specification, stdlib release v0.8.1

Libraries for other jobs follow the same pattern as BLAS and LAPACK: a compiled library, an interface to write or include, and a link line. Most are written in C and give Fortran a thin binding. Those named below are well-known community and vendor projects, not part of the language.

(HDF5, NetCDF and PETSc were not installed or run for this page; the table lists what each project documents, with the project pages in the references.) The Fortran 2003 FFTW interface is the one to look at, because it shows the C interoperability tools doing real work; see Interoperability with C. This program transforms a row of 16 samples containing three waves and prints the bins with energy:

Compiled with gfortran 09_fftw.f90 -I/usr/include -lfftw3 and the same flags for Flang; both print

The bin is 3 and its magnitude is 8, which is n/2 for a unit-amplitude wave. Three details: the include 'fftw3.f03' file is found through -I, the plan is made before the data is written (planning with some flags overwrites the arrays), and the arrays come from fftw_alloc_real so that they are aligned for the vector instructions. With -Wall -Wextra the vendor header produces about forty -Wunused-parameter warnings (one per constant it declares that the program does not use), which is a property of the header, not of the program, and -Wno-unused-parameter silences them.

Source: FFTW manual: Calling FFTW from Modern Fortran

The two Fortran books consulted for this section treat numerical libraries as an application of the language and touch BLAS and LAPACK only through worked examples. Where older teaching material differs from current practice:

Not checked here, because no compiler for them is installed: Intel ifx, NAG, HPE Cray and NVIDIA HPC compilers; and oneMKL, which was not installed or run. See History, Standards and Compilers for the per-compiler support table.

See also Best Practices and Style and Common Gotchas and Pitfalls.