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:

LAPACK webpage: https://www.netlib.org/lapack/
Online documentation for LAPACK Fortran functions: http://www.netlib.org/lapack/explore-html/
LAPACK Users' Guide: http://www.netlib.org/lapack/lug/
LAPACK GitHub: https://github.com/Reference-LAPACK/lapack
Man pages for LAPACK functions on Perlmutter

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

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

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