Introduction

LAPACK and BLAS are originally written in Fortran and meant to be used in Fortran programs. Many vendors supply an optimised version of the LAPACK and BLAS libraries. Naturally, people who are programming in C or C++ and want to make use of the efficient implementation of the LAPACK/BLAS libraries. Here we will demonstrate how this can be done.

At the end you find a complete working example, together with a script to run the same program using various interfaces and LAPACK libraries.

Fortran example of dsyev

The examples wil be based on the call of dsyev, the fortran code to do this is:

      character*1 jobz, uplo
      integer n,lda,lwork,info
      double precision, allocatable a(:,:)
      double precision, allocatable w(:)
      double precision, allocatable, work(:)
      double precision wwork

      n=1000
      lda=1000
      jobz='V'
      uplo='U'
      allocate(a(lda,n))
      ... code to fill matrix a ...
      allocate(w(n))
      lwork = -1       ! workspace query
      call dsyev(jobz, uplo, n, a, lda, w, wwork, lwork, info)
      lwork = wwork+0.1
      ! now the real work:
      allocate(work(lwork))
      call dsyev(jobz, uplo, n, a, lda, w, work, lwork, info)
      deallocate(work)

Using LAPACK as-is

You can use the LAPACK library as-is, by mimicking the way Fortran calls Fortran subroutines. Three points are important here:

• use dsyev_ in stead of dsyev
• every parameter must be a pointer
• be aware of strings of varying length: Fortran has it's own way to deal with strings in parameter lists. In general it is safe to use a pointer to a character if the LAPACK subroutine expects only one character.
• matrices in C and Fortran differ: C is row-major, while Fortran is colomn-major.

So, the C-version of above Fortran fragment would become:

void dsyev_(char *jobz, char *uplo, int *n, double *a, int *lda, 
           double *w, double *work, int *lwork, int *info);
 ... 
char       jobz, uplo;
int  n,lda,lwork,info;
double    *a;
double    *w;
double *work;
double wwork;

n     = 1000;
lda   = 1000;
jobz  =  'V';
uplo  =  'U';
a     = malloc(sizeof(double)*n*lda);
w     = malloc(sizeof(double)*n);
lwork = -1;  // workspace query
dsyev_(&jobz, &uplo, &n, a, &lda, w, &wwork, &lwork, &info);
lwork = wwork+0.1;
work  = malloc(sizeof(double)*lwork);
// now the real thing:
dsyev_(&jobz, &uplo, &n, a, &lda, w, work, &lwork, &info);
free(work);

Using LAPACKE: an interface between C and Fortran

The newer versions of LAPACK from www./netlib.org/lapack , contain an interface to use LAPACK from C. How to use this interface is described in detail in www.netlib.org/lapack/lapacke , we will give here an example, based on the Fortran code above.

Important things are:

• The name of the subroutine becomes LAPACKE_dsyev
• The function returns the info argument, the info argument is not in the parameter list
• All parameters are used in a C-like fashion: only the arrays are pointers
• In general, there are two versions of the function: for example LAPACKE_dsyev and LAPACKE_dsyev_work. The first one determines the optimal workspace itself, the second one expects that the workspace is provided by the caller. The first one must be called without workspace arrays and dimensions. We will demonstrate the first one here.
• If one or more parameters are 2-dimension matrices, an extra parameter is needed: LAPACK_COL_MAJOR or LAPACK_ROW_MAJOR. The first one for Fortran matrix layout, the second for C matrix layout. This parameter is placed before all other parameters.
• LAPACKE can also be used to call LAPACK subroutines from Intel's MKL library, see below.

Here the example:

#include <lapacke.h>
 ... 
char jobz, uplo;
int  n,lda,info;
double    *a;
double    *w;

n     = 1000;
lda   = 1000;
jobz  =  'V';
uplo  =  'U';
a     = malloc(sizeof(double)*n*lda);
w     = malloc(sizeof(double)*n);
info  = LAPACKE_dsyev(LAPACKE_COL_MAJOR,jobz, uplo, n, a, lda, w);

An advantage is, that an header file is provided which prevents trivial mistakes in calling the function.

Working example codes

Here we provide you with working example codes that you can run on the Lisa system. The examples include the usage of MKL from C programs.

Using copy/paste, copy both files (testdsyev.c and doit.txt, see 'attached files' below) to Lisa and give the following commands:

chmod +x doit.txt
./doit.txt

In the output you can see which method is used (LAPACK_PLAIN or LAPACKE), the order of the matrix used; the value of info after the function call and the time it took. In the script doit.txt you can study the linking flags that are used.