Skip to content

LAPACK

LAPACK - Fortran Interface

Linear Algebra PACKage (LAPACK) provides Fortran 90 routines for solving systems of simultaneous linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems. The associated matrix factorizations (LU, Cholesky, QR, SVD, Schur, generalized Schur) are also provided, as are related computations such as reordering of the Schur factorizations and estimating condition numbers. Dense and banded matrices are handled, but not general sparse matrices. In all areas, similar functionality is provided for real and complex matrices, in both single and double precision.

Useful documentation can be found from the following sources:

LAPACKE - C Interface

LAPACK includes the LAPACKE package, a C language API for LAPACK. You can find info on how to use the package by reading the LAPACKE webpage: The LAPACKE C Interface to LAPACK.

Optimized Versions

HPE Cray provides an optimized LAPACK version in its scientific math library called LibSci. When building a code with the HPE Cray compiler wrappers (that is, ftn, cc and CC), this library is linked automatically and you don’t have to explicitly add additional flag to use LAPACK.

The Intel MKL (Math Kernel Library) also provides an optimized LAPACK version, too. Unlike Cray Libsci, however, you have to add a special flag in order to use MKL. You can find info on how to do this from the MKL page.

Examples

The following example code solves a linear equation, $A \ \mathbf{x} = \mathbf{b}$, by factorizing the matrix $A$ first and then solving for the unknown $\mathbf{x}$. It uses the general matrix type version of the LAPACK routines.

The code employs a special matrix type described in https://www.scribd.com/document/397967102/LAPACK-and-ScaLAPACK-Examples. When such a matrix is multiplied by the vector $\mathbf{x} = \left(1\ 2\ 3\ \cdots n\right)^T$ where $n$ is the dimension of the matrix, the result vector $\mathbf{b}$ is given by $(n + 1) \left(1\ 1\ 1 \cdots 1\ 2\right)^T$. This property holds for any $n$ greater than 1, although the author of the contents in the aforementioned webpage restricts to the even number cases only.

```fortran program lapacktest_gen

! Solve A x = b for x, using LAPACK functions. ! ! NERSC

implicit none integer n double precision, allocatable :: a(:,:), x(:) integer, allocatable :: ipiv(:) integer info integer i

! Set up; assemble A and RHS

print , ‘Enter the matrix dimension, n (n > 1)’ read , n

allocate(a(n,n), x(n)) allocate(ipiv(n))

a = 0.d0 do i = 1, n if (i > 1) then a(i-1, i) = a(i-1, i) - 1.d0 end if a(i,i) = a(i,i) + 3.d0 if (i + 1 <= n) then a(i+1, i) = a(i+1, i) - 1.d0 end if a(n+1-i, i) = a(n+1-i, i) + 1.d0 end do

x(:n-1) = (n + 1.d0) x(n) = (n + 1.d0) * 2.d0

! Factorize

call dgetrf(n, n, a, n, ipiv, info)

if (info /= 0) then print *, ‘Error with factorization:’, info stop end if

! Solve

call dgetrs(‘N’, n, 1, a, n, ipiv, x, n, info)

if (info /= 0) then print *, ‘Error with solver:’, info stop end if

print *, (x(i), i=1,n)

deallocate(a, x) deallocate(ipiv)

end ```

To build and run using HPE Cray’s Libsci:

```console $ ftn -o lapacktest_gen lapacktest_gen.f90

$ ./lapacktest_gen Enter the matrix dimension, n (n > 1) 10 1.00000000000000 2.00000000000000 3.00000000000000 4.00000000000000 5.00000000000000 6.00000000000000 7.00000000000000 8.00000000000000 9.00000000000000 10.0000000000000 ```

With a little bit of algebra, it can be shown that the matrix is actually positive definite. So, we can use the positive-definite version, instead:

```fortran program lapacktest_pos

! Solve A x = b for x, using LAPACK functions. ! The matrix is positive definite, and positive definite functions ! are used here. ! ! NERSC

implicit none integer n double precision, allocatable :: a(:,:), x(:) integer info integer i

! Set up; assemble A and RHS

print , ‘Enter the matrix dimension, n (n > 1)’ read , n

allocate(a(n,n), x(n))

a = 0.d0 do i = 1, n if (i > 1) then a(i-1, i) = a(i-1, i) - 1.d0 end if a(i,i) = a(i,i) + 3.d0 if (i + 1 <= n) then a(i+1, i) = a(i+1, i) - 1.d0 end if a(n+1-i, i) = a(n+1-i, i) + 1.d0 end do

x(:n-1) = (n + 1.d0) x(n) = (n + 1.d0) * 2.d0

! Factorize

call dpotrf(‘U’, n, a, n, info)

if (info /= 0) then print *, ‘Error with factorization:’, info stop end if

! Solve

call dpotrs(‘U’, n, 1, a, n, x, n, info)

if (info /= 0) then print *, ‘Error with solver:’, info stop end if

print *, (x(i), i=1,n)

deallocate(a, x)

end ```

You can build and run as before:

```console $ ftn -o lapacktest_pos lapacktest_pos.f90

$ ./lapacktest_pos Enter the matrix dimension, n (n > 1) 10 0.999999999999999 2.00000000000000 3.00000000000000 4.00000000000000 5.00000000000000 6.00000000000000 7.00000000000000 8.00000000000001 9.00000000000000 10.0000000000000 ```

An instruction on how to build the code with the MKL library can be found in the MKL webpage. When you build with the Cray compiler wrapper, you can see a warning message as shown below, but you can ignore it. The numerical result should be basically the same.

text Warning: Headers and libraries from cray-libsci/XX.XX.X will be ignored because they conflict with -mkl.

Below is a C code for solving the same linear equation using LAPACKE functions.

```C /* Solve A x = b for x, using LAPACKE functions. The matrix is positive definite, and positive definite functions are used here.

NERSC */

include

include

include

int main (int argc, const char * argv[]) { double a, x; lapack_int info, n; int i;

/ Set up; assemble A and RHS /

printf(“Enter the matrix dimension, n (n > 1)\n”); scanf(“%d”, &n);

a = calloc(n*n, sizeof(double)); x = calloc(n, sizeof(double));

for (i=0; i 0) a[i - 1 + i * n] -= 1.; a[i + i * n] += 3.; if (i + 1 < n) a[i + 1 + i * n] -= 1.; a[n - 1 - i + i * n] += 1.; x[i] = n + 1.; } x[n-1] = (n + 1.) * 2.;

/ Factorize /

info = LAPACKE_dpotrf(LAPACK_COL_MAJOR, ‘U’, n, a, n);

if (info != 0) { printf(“Error with factorization: %d\n”, info); exit (1); }

/ Solve /

info = LAPACKE_dpotrs(LAPACK_COL_MAJOR, ‘U’, n, 1, a, n, x, n);

if (info != 0) { printf(“Error with solver: %d\n”, info); exit (1); }

for (i=0; i<n; i++) { printf(” %f\n”, x[i]); }

free(a); free(x);

return 0; } ```

You can build and run as follows:

```console $ cc -o lapacketest_pos lapacketest_pos.c

$ ./lapacketest_pos Enter the matrix dimension, n (n > 1) 10 1.000000 2.000000 3.000000 4.000000 5.000000 6.000000 7.000000 8.000000 9.000000 10.000000 ```

To use MKL instead of LibSci, you need to use mkl_lapacke.h instead of lapacke.h - just replace the line #include <lapacke.h> with #include <mkl_lapacke.h>. Note also that you can see a warning message when you build using the Cray compiler wrapper cc. Just ignore it.

text Warning: Headers and libraries from cray-libsci/XX.XX.X will be ignored because they conflict with -mkl.