# 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.

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)'

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:

$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:

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)'

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:

$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.

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.

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

NERSC
*/

#include <stdio.h>
#include <stdlib.h>

#include <lapacke.h>

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<n; i++) {
if (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:

$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.

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