Write a Solver for Systems of Equations Using LAPACK on Mac OS

Updated September 2026: a new example with a non-zero answer, a note on column-major order, the -framework Accelerate flag, the new LAPACK interface in recent macOS versions, and the same solver in Python with NumPy and SciPy.

Apple provides the Accelerate library which includes the linear algebra library (LAPACK). In this post, we will go through writing a simple C++ program to use this library on Mac OS.

The Equations¶

Let’s solve the following equations:

We can rewrite the equations as:

The answer is $x_0=1$ and $x_1=2$, so we can check the program’s output.

Column-Major Order¶

LAPACK comes from Fortran, and it expects matrices stored column by column. C and C++ arrays are usually written row by row. So the matrix above goes into the array as its first column (2, 3) followed by its second column (1, -1):

double a[2*2] = {2, 3,
                 1, -1};

If you write it row by row as {2, 1, 3, -1}, LAPACK solves the transposed system and gives you a wrong answer without any error.

The Code¶

Now, let’s write our code. dgetrf_ computes the LU factorization of the matrix, and dgetrs_ uses that factorization to solve for the right-hand side. dgetrs_ writes the solution into b:

#include <iostream>
#include <Accelerate/Accelerate.h>

int main()
{
    int number_of_rows = 2;
    int number_of_cols = 2;
    int number_of_right_hand_side_cols = 1;
    int LDA = 2;
    int LDB = 2;
    int IPIV[2];
    int INFO = 0;
    char TRANS = 'N'; // Non transpose

    // column-major: first column (2, 3), second column (1, -1)
    double a[2*2] = {2, 3,
                     1, -1};
    double b[2] = {4,
                   1};

    dgetrf_(&number_of_rows, &number_of_cols, a, &LDA, IPIV, &INFO);
    if (INFO != 0) {
        std::cout << "LU factorization failed, INFO = " << INFO << std::endl;
        return 1;
    }
    std::cout << "LU factorization executed without errors." << std::endl;

    dgetrs_(&TRANS, &number_of_rows, &number_of_right_hand_side_cols,
            a, &LDA, IPIV, b, &LDB, &INFO);
    if (INFO != 0) {
        std::cout << "Solver failed, INFO = " << INFO << std::endl;
        return 1;
    }

    std::cout << "Result: x_0 = " << b[0] << ", x_1 = " << b[1] << std::endl;
    return 0;
}

IPIV holds the row swaps from the factorization, and it needs one entry per row. INFO is 0 when a call succeeds. A positive INFO from dgetrf_ means the matrix is singular, so there is no unique solution.

Compile and Run¶

Save the code as lapack_test.cpp, then compile it against the Accelerate framework and run it:

clang++ lapack_test.cpp -framework Accelerate -o ltest
./ltest

The output is:

LU factorization executed without errors.
Result: x_0 = 1, x_1 = 2

g++ works too. On macOS it is the same Apple Clang compiler.

The New LAPACK Interface¶

On macOS 13.3 and newer, the compiler warns that these LAPACK functions are deprecated. Apple updated the LAPACK in Accelerate to version 3.9.1 and kept the old version for existing programs. To use the new one, define ACCELERATE_NEW_LAPACK when you compile:

clang++ lapack_test.cpp -DACCELERATE_NEW_LAPACK -framework Accelerate -o ltest

The code stays the same. If you work with very large matrices, also add -DACCELERATE_LAPACK_ILP64 to use 64-bit integers for sizes and indices. In that case, declare the integer variables in the code as __LAPACK_int instead of int.

The Same Solver in Python¶

If you don’t need C++, you can do the same thing in Python with NumPy and SciPy. They don’t have their own linear algebra code for this. Under the hood, they call the same LAPACK routines. On macOS 14 and newer, the NumPy and SciPy packages from pip are built against Apple’s Accelerate, so they use the same LAPACK as the C++ program above. On other systems, they usually use OpenBLAS.

Install them with pip:

pip install numpy scipy

NumPy¶

numpy.linalg.solve calls LAPACK’s dgesv, which does the LU factorization and the solve in one step:

import numpy as np

A = np.array([[2.0, 1.0],
              [3.0, -1.0]])
b = np.array([4.0, 1.0])

x = np.linalg.solve(A, b)
print(x)  # [1. 2.]

Notice that the matrix is written row by row here. NumPy takes care of the column-major order before it calls LAPACK.

SciPy with LU Factorization¶

To keep the two steps separate like in the C++ code, use lu_factor and lu_solve. They call dgetrf and dgetrs. This is useful when you solve for many right-hand sides with the same matrix, because the factorization runs only once:

from scipy.linalg import lu_factor, lu_solve

lu, piv = lu_factor(A)       # dgetrf
x = lu_solve((lu, piv), b)   # dgetrs
print(x)  # [1. 2.]

Calling LAPACK Directly¶

SciPy can also give you the LAPACK functions themselves. get_lapack_funcs picks the right version of a routine for your data type: dgetrf for float64, sgetrf for float32, and zgetrf for complex numbers. The calls look almost the same as in C++:

from scipy.linalg import get_lapack_funcs

getrf, getrs = get_lapack_funcs(("getrf", "getrs"), (A,))
print(getrf.typecode)  # d, so these are dgetrf and dgetrs

lu, piv, info = getrf(A)
if info != 0:
    raise RuntimeError(f"LU factorization failed, INFO = {info}")

x, info = getrs(lu, piv, b)
if info != 0:
    raise RuntimeError(f"Solver failed, INFO = {info}")

print(x)  # [1. 2.]

INFO means the same thing as in the C++ version. One difference: SciPy returns piv with 0-based row indices, while Fortran LAPACK uses 1-based indices. You don’t need to change anything as long as you pass piv from getrf straight to getrs.

If you know the exact routine you want, you can also import it directly with from scipy.linalg.lapack import dgetrf, dgetrs.

To see which LAPACK your NumPy uses, run:

import numpy as np
np.show_config()

Look for accelerate or openblas under the BLAS and LAPACK entries.