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_MAJORrefers 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_MAJORconstant when calling LAPACK functions. Conversely,LAPACK_COL_MAJORis 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