UMBC logo
UMBC High Performance Computing Facility
How to use LAPACK on tara

Introduction

The Linear Algebra PACKage (LAPACK) library is useful for matrix algebra operations like singular value decomposition, solving systems of linear equations, and computing eigenvalues. You may want to consider using it instead of writing your own routines. For lower level operations like matrix multiplication, see BLAS. In this tutorial, we will compile and run a program that uses LAPACK to compute SVDs. Before you begin, make sure to read the tutorial for compiling C programs.

Where to find documentation

On this page, we'll show how to use LAPACK functionality within AMD Core Math Library (ACML). We'll assume that you're using the default "pgi-9.0-mvapich2-1.4rc2" switcher setting - other MPI+compiler combinations may not work with the code below.

For more information about the LAPACK routines in ACML, see

Example

In this example we will create a matrix, then compute its singular value decomposition (SVD) using the LAPACK dgesvd function. The SVD of A (m x n) will consist of U (m x m), S (p x 1, where p = min(m,n)), and VT (n x n) Our matrices will contain doubles, and be stored in column-major order. Here is the code
#include <stdio.h>
#include <acml.h>

#define MATRIX_IDX(n, i, j) j*n + i
#define MATRIX_ELEMENT(A, m, n, i, j) A[ MATRIX_IDX(m, i, j) ]

void init_matrix(double* A, int m, int n)
{
   double element = 0.0;
   for (int j = 0; j < n; j++)
   {
      for (int i = 0; i < m; i++)
      {
         MATRIX_ELEMENT(A, m, n, i, j) = 1.0 / (i + j + 1);
         element += 0.25;
      }
   }
}

void print_matrix(const double* A, int m, int n)
{
   for (int i = 0; i < m; i++)
   {
      for (int j = 0; j < n; j++)
      {
          printf("%8.4f", MATRIX_ELEMENT(A, m, n, i, j));
      }
      printf("\n");
   }
}

int main(int argc, char** argv)
{
   int m = 3;
   int n = 4;

   int p = (m < n ? m : n);

   double A[m * n];
   double U[m * m];
   double VT[n * n];
   double S[p * 1];

   int info;

   init_matrix(A, m, n);

   printf("Matrix A (%d x %d) is:\n", m, n);
   print_matrix(A, m, n);

   dgesvd('A', 'A', m, n, A, m, S, U, m, VT, n, &info);
   if (info != 0)
   {
      fprintf(stderr, "Warning: dgesvd returned with a non-zero status (info = %d)\n", info);
   }

   printf("\nMatrix U (%d x %d) is:\n", m, m);
   print_matrix(U, m, m);

   printf("\nVector S (%d x %d) is:\n", p, 1);
   print_matrix(S, p, 1);

   printf("\nMatrix VT (%d x %d) is:\n", n, n);
   print_matrix(VT, n, n);

   return 0;
}


Download: ../code-2010/lapack-svd/pgi_acml/main.c
The important line to notice, where the actual SVD computation is happening, is this
dgesvd('A', 'A', m, n, A, m, S, U, m, VT, n, &info);
where the 'A' arguments are flags saying that we want the entire U and VT matrices filled in by the procedure. For more information, see the documentation at Netlib, or try "man dgesvd" at the command line. Here is the Makefile
EXECUTABLE := matrix_svd
OBJS := main.o

CC := mpicc
CFLAGS := -O3 -c99 -Minform=warn -fastsse

INCLUDES :=
LIBLOCS :=
LDFLAGS := -lm -lacml -lpgftnrtl

%.o: %.c %.h
    $(CC) $(CFLAGS) $(DEFS) $(INCLUDES) -c $< -o $@

$(EXECUTABLE): $(OBJS)
    $(CC) $(CFLAGS) $(DEFS) $(INCLUDES) $(OBJS) -o $@ $(LIBLOCS) $(LDFLAGS)

clean:
    -rm -f *.o $(EXECUTABLE)


Download: ../code-2010/lapack-svd/pgi_acml/Makefile
The important part is that we need to link to two more libraries: libacml and libpgftnrtl. This is accomplished by adding them to the LDFLAGS variable. Compiling the code should look something like this
[araim1@tara-fe1 pgi_acml]$ make
mpicc -O3 -c99 -Minform=warn -fastsse   -c -o main.o main.c
mpicc -O3 -c99 -Minform=warn -fastsse   main.o -o matrix_svd  -lm -lacml -lpgftnrtl
[araim1@tara-fe1 pgi_acml]$ ls
main.c  main.o  Makefile  matrix_svd
[araim1@tara-fe1 pgi_acml]$ 
Then running the code should look this
[araim1@tara-fe1 pgi_acml]$ ./matrix_svd
Matrix A (3 x 4) is:
  1.0000  0.5000  0.3333  0.2500
  0.5000  0.3333  0.2500  0.2000
  0.3333  0.2500  0.2000  0.1667

Matrix U (3 x 3) is:
 -0.8199  0.5563  0.1349
 -0.4662 -0.5123 -0.7213
 -0.3322 -0.6543  0.6794

Vector S (3 x 1) is:
  1.4519
  0.1433
  0.0042

Matrix VT (4 x 4) is:
 -0.8015 -0.4466 -0.3143 -0.2435
  0.5729 -0.3919 -0.5127 -0.5053
  0.1692 -0.7398  0.1245  0.6392
 -0.0263  0.3157 -0.7892  0.5261
[araim1@tara-fe1 pgi_acml]$ 
Running the code through the batch system should require no special steps. A standard SLURM script like the one below should be sufficient.
#!/bin/bash
#SBATCH --job-name=matrix_svd
#SBATCH --output=slurm.out
#SBATCH --error=slurm.err
#SBATCH --partition=develop
#SBATCH --nodes=1
#SBATCH --ntasks-per-node=1

./matrix_svd

Download: ../code-2010/lapack-svd/pgi_acml/run.slurm