LAPACK

LAPACK (Linear Algebra PACKage) is a highly optimized library designed for solving linear algebra problems such as solving systems of linear equations, computing matrix factorizations, and performing matrix inversions. LAPACK is written in Fortran, but there are wrappers and bindings available for C and other languages.

Installing LAPACK

To use LAPACK, you need to have LAPACK and BLAS (Basic Linear Algebra Subprograms) installed.

Linux

You can install them with a package manager, for example:

sudo apt-get install liblapack-dev libblas-dev

MacOS

brew install lapack

This will also install any necessary dependencies, including BLAS.

To work with LAPACK in C, you will want to install the LAPACKe package, which provides the necessary bindings for C:

brew install lapacke

To make sure LAPACK is installed correctly, you can check if lapacke.h and the shared libraries (liblapacke.dylib, liblapack.dylib, libblas.dylib) are available in the default Homebrew installation path (/usr/local/include for headers and /usr/local/lib for libraries). You can also list installed libraries by running:

ls /usr/local/lib | grep lapack
ls /usr/local/lib | grep blas

Compiling

When compiling your program, you need to link LAPACK and BLAS libraries. For example:

gcc -o matrix_inverse matrix_inverse.c -llapacke -llapack -lblas -lm

or

gfortran -o matrix_inverse matrix_inversion.f90 -llapack -lblas

  • -llapacke: Links the LAPACKe C interface.
  • -llapack: Links the LAPACK library.
  • -lblas: Links the BLAS library, which LAPACK depends on.
  • -lm: Links the math library (for operations like square roots or other math functions).

[!WARNING] If the libraries are installed in a non-standard path (e.g., /opt/homebrew/ on newer Macs with M1/M2 chips), you may need to specify the location of the libraries explicitly when compiling, like this:

gcc -o matrix_inverse matrix_inverse.c -I/opt/homebrew/include -L/opt/homebrew/lib -llapacke -llapack -lblas -lm

LAPACK naming schemes

The name of each LAPACK routine is a coded specification of its function. All driver and computational routines have names of the form XYYZZZ, where for some driver routines the 6th character is blank.

The first letter, X, indicates the data type as follows:

S	REAL
D	DOUBLE PRECISION
C	COMPLEX
Z	COMPLEX*16 or DOUBLE COMPLEX

When we wish to refer to an LAPACK routine generically, regardless of data type, we replace the first letter by x. Thus xGESV refers to any or all of the routines SGESV, CGESV, DGESV and ZGESV.

The next two letters, YY, indicate the type of matrix (or of the most significant matrix). Most of these two-letter codes apply to both real and complex matrices; a few apply specifically to one or the other, as indicated in the table below.

Open Table
BD	bidiagonal
DI	diagonal
GB	general band
GE	general (i.e., unsymmetric, in some cases rectangular)
GG	general matrices, generalized problem (i.e., a pair of general matrices)
GT	general tridiagonal
HB	(complex) Hermitian band
HE	(complex) Hermitian
HG	upper Hessenberg matrix, generalized problem (i.e a Hessenberg and a
 	triangular matrix)
HP	(complex) Hermitian, packed storage
HS	upper Hessenberg
OP	(real) orthogonal, packed storage
OR	(real) orthogonal
PB	symmetric or Hermitian positive definite band
PO	symmetric or Hermitian positive definite
PP	symmetric or Hermitian positive definite, packed storage
PT	symmetric or Hermitian positive definite tridiagonal
SB	(real) symmetric band
SP	symmetric, packed storage
ST	(real) symmetric tridiagonal
SY	symmetric
TB	triangular band
TG	triangular matrices, generalized problem (i.e., a pair of triangular matrices)
TP	triangular, packed storage
TR	triangular (or in some cases quasi-triangular)
TZ	trapezoidal
UN	(complex) unitary
UP	(complex) unitary, packed storage

The last three letters ZZZ indicate the computation performed. For example, SGEBRD is a single precision (S) routine that performs, on a general matrix (GE), a bidiagonal reduction (BRD).

Matrix Inversion

Here’s a guide on how to use LAPACK in C to compute the inverse of a matrix.

To compute the inverse of a matrix , you first perform an LU decomposition of the matrix and then use the result to compute the inverse. LAPACK provides the dgetrf function for LU decomposition and dgetri for matrix inversion.

Below is an example of how to invert a matrix using LAPACK in C:

#include <stdio.h>
#include <lapacke.h>

int main() {
    // Define the matrix A (for example, a 3x3 matrix)
    int n = 3;
    double A[9] = {
        4.0, 2.0, 1.0,
        8.0, 7.0, 3.0,
        2.0, 1.0, 5.0
    };

    // Define pivot array and workspace
    int ipiv[3];  // Pivot indices array
    int info;     // Status info

    // Perform LU decomposition
    info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, n, n, A, n, ipiv);
    if (info != 0) {
        printf("LU decomposition failed with error code %d\n", info);
        return -1;
    }

    // Perform matrix inversion using LU decomposition
    info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, A, n, ipiv);
    if (info != 0) {
        printf("Matrix inversion failed with error code %d\n", info);
        return -1;
    }

    // Output the inverse matrix
    printf("Inverse of A:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            printf("%6.2f ", A[i * n + j]);
        }
        printf("\n");
    }

    return 0;
}

LAPACK_ROW_MAJOR is a constant or macro used when calling LAPACK functions to specify the memory layout of the matrix data. It tells LAPACK that the matrix is stored in row-major order, meaning that the matrix elements are stored in a contiguous block of memory row by row.

Consider this matrix:

  • Row-major order: All elements of a row are stored consecutively in memory. In C, matrices are typically stored in row-major order by default. This is what LAPACK_ROW_MAJOR refers to. In row-major order, the elements would be stored in memory like this:

  • Column-major order: All elements of a column are stored consecutively in memory. Fortran, and hence LAPACK’s native format, uses column-major order by default. So in this order, the above matrix would be stored as:

[!NOTE] Why LAPACK_ROW_MAJOR?

Since LAPACK is traditionally written in Fortran (which uses column-major order), C and C++ users need to specify how their matrices are stored. LAPACK functions can handle both row-major and column-major layouts, but you need to tell LAPACK how your data is arranged in memory. To inform LAPACK that you’re passing a matrix stored in row-major order, you’d use the LAPACK_ROW_MAJOR constant when calling LAPACK functions. Conversely, LAPACK_COL_MAJOR is used if the matrix is stored in column-major order (the Fortran default).

Below is an example of how to invert a matrix using LAPACK in FORTRAN:

program matrix_inversion
   implicit none
   real(8), dimension(2,2) :: A
   real(8), dimension(2) :: ipiv
   integer :: info

   ! Matrix A to be inverted
   A = reshape((/ 4.0d0, 7.0d0, 2.0d0, 6.0d0 /), (/ 2, 2 /))

   ! Call LAPACK routine to perform LU decomposition and matrix inversion
   call dgetrf(2, 2, A, 2, ipiv, info) ! LU decomposition
   if (info /= 0) then
      print *, 'Error in LU decomposition, info = ', info
      stop
   endif

   call dgetri(2, A, 2, ipiv, info)    ! Matrix inversion
   if (info /= 0) then
      print *, 'Error in matrix inversion, info = ', info
      stop
   endif

   print *, 'Inverse of matrix A:'
   print *, A
end program matrix_inversion