Linking against libraries

We saw in the Fortran chapter that we could expand the functionality of our program greatly by using external libraries, like Lapack and BLAS. Similarly, in Python we imported other modules like NumPy and Matplotlib. Of course, we can do the same in C. In fact, we’ve already done this by including the math library.

Recall our first little sample program:

/* File: pythagoras.c
 * Author: Ian May
 * Purpose: Demonstrate overall structure of C program as well as
 *          how to call gcc
 */

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

double pythag(double a, double b)
{
  return sqrt(a*a + b*b);
}

int main(int argc, char **argv)
{
  /* Compute several hypo. lengths */
  double a1 = 1.2,b1 = 3.3;
  double c1 = pythag(a1,b1);
  
  /* Note that unlike Fortran 90, we can define variables anywhere */
  double a2 = -3.0,b2 = 4.0;
  double c2 = pythag(a2,b2);

  // Similar to Python, we can do formatted output using specifiers
  printf("Triangle with sides %f and %f has hypotenuse %f\n",a1,b1,c1);
  printf("Triangle with sides %f and %f has hypotenuse %f\n",a2,b2,c2);

  /* You can also nest function calls together */
  double a3 = 5.0,b3 = 12.0;
  printf("Triangle with sides %f and %f has hypotenuse %f\n",a3,b3,pythag(a3,b3));

  /* We return 0 to tell the OS that everything ran normally */
  return 0;
}

Download this code

Which we compiled by calling:

gcc pythagoras.c -o pythagoras.ex -lm

Notably, the flag -lm tells gcc to find a file called libm.a. The .a extension means that this is a library that can be linked statically. For all intents and purposes, you can think of static libraries as additional .o files, that provide definitions for symbols referred to your source files.

If the library is also written in C, there will usually be a header file associated to the library that fills in the interfaces to the supplied functions, and defines any necessary data structures. This is essentially the pattern we used for our column-major matrix example. If the library is not written in C, e.g. Lapack, then there generally will not be an associated header file. We’ll deal with this later.

Next, we’ll look at a few commonly used C libraries for scientific computing.

A few useful libraries

• The math library. You can link against this using -lm, and the math.h header.

• The gcc quadruple precision library. This is linked by using -lquadmath, and the quadmath.h header.

• The GNU scientific library provides a huge amount of capability.

  • This is organized into a group of separate libraries, so that you only need to link the relevant items. Compare this to SciPy.

  • We aren’t in a position to talk much about software licenses in this class, but note that it is something to be aware of. Fortunately, most licenses are totally compatible with academic goals.

There are many other useful libraries, and we’ll revisit this topic towards the end of the quarter, but these provide just a couple.

Creating your own (static) libraries

Consider this fairly common circumstance: You’ve written a big code spread over several source files, and perhaps some of the files depend on each other. This code base exists to solve one larger problem, or class of problems. You may find that you want to write multiple programs that all use these same source files. Shuffling .o files around isn’t a very good solution, and neither is copying all of the sources to every new project that uses them.

This is of course the exact reason libraries exist. Let’s revisit our column major matrix example, and add two more files.

A header file:

/* File: Hilbert.h
 * Author: Ian May
 * Purpose: Specify interface for a Hilbert matrix constructor
 */

#include "ColMajorMat.h"

#ifndef HILBERT_H
#define HILBERT_H

ColMajorMat CreateHilbertMatrix(int,int);

#endif

Download this code

and a matching source file:

/* File: Hilbert.c
 * Author: Ian May
 * Purpose: Add a function that creates a Hilbert matrix in the CM format
 */

#include <stdlib.h>

#include "ColMajorMat.h"

/* Matrix creation handles allocation */
ColMajorMat CreateHilbertMatrix(int m, int n)
{
  /* Create the matrix */
  ColMajorMat mat = CreateMatrix(m,n);

  /* Fill this matrix like so */
  for (int j=0; j<n; j++) {
    for (int i=0; i<m; i++) {
      mat->data[j*m + i] = 1./(i+j+1); /* Note zero indexing */
    }
  }
  
  /* Pass the matrix back out */
  return mat;
}

Download this code

The remaining files are nearly the same as before. The column major format header is:

/* File: ColMajorMat.h
 * Author: Ian May
 * Purpose: Another iteration on the previous example, now modeling the OOC
 *          approach
 */

#ifndef COLMAJORMAT_H
#define COLMAJORMAT_H

/* We will only use a pointer to the struct, so make that the typedef */
typedef struct _p_ColMajorMat *ColMajorMat;

/* Define a structure that holds the matrix dimensions and a flattened array */
struct _p_ColMajorMat{
  int m, n;
  double *data;
};

/* Functions to create and destroy the matrix struct */
ColMajorMat CreateMatrix(int,int);
void DestroyMatrix(ColMajorMat);

/* Functions to interact with this matrix */
void MatVecMult(ColMajorMat,double[*],double[*]);

#endif

Download this code

The matching source file is:

/* File: ColMajorMat.c
 * Author: Ian May
 * Purpose: Another iteration on the previous example, now modeling the OOC
 *          approach
 */

#include <stdlib.h>

#include "ColMajorMat.h"

/* Matrix creation handles allocation */
ColMajorMat CreateMatrix(int m, int n)
{
  /* Allocate and fill struct */
  ColMajorMat mat = malloc(sizeof(*mat));
  mat->m = m;
  mat->n = n;
  mat->data = malloc(m*n*sizeof(*mat->data));
  /* Pass this struct back out */
  return mat;
}

/* Destruction frees everything creation allocated */
/* Can you see a possible bug here? */
void DestroyMatrix(ColMajorMat A)
{
  /* Free the data array inside the struct */
  free(A->data);
  /* Free the struct */
  free(A);
}

void MatVecMult(ColMajorMat A, double x[A->m], double b[A->m])
{
  for (int i=0; i<A->m; i++) {
    b[i] = 0;
    for (int j=0; j<A->n; j++) {
      b[i] += A->data[j*A->m + i]*x[j];
    }
  }
}

Download this code

and the main file is now:

/* File: Main.c
 * Author: Ian May
 * Purpose: Slight modification to prior example to now use a Hilbert matrix
 */

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

#include "Hilbert.h"

int main(int argc, char **argv)
{
  /* Matrix sizes, easily set from command line */
  int m = argc>1 ? atoi(argv[1]) : 4;
  int n = argc>2 ? atoi(argv[2]) : 3;
  
  /* Arrays to hold vectors we are using */
  double b[m];
  double x[n];
  for (int j=0; j<n; j++) {
    x[j] = 1.0/(1.0 + j*j);
  }

  /* Create and fill a Hilbert matrix */
  ColMajorMat A = CreateHilbertMatrix(m,n);

  /* Store A*x into b, interpreted as matrix-vector product */
  MatVecMult(A,x,b);

  /* Print out result */
  printf("b =");
  for (int i=0; i<m; i++) {
    printf(" %1.2f |",b[i]);
  }
  printf("\n");

  /* Don't forget to clean up! */
  DestroyMatrix(A);
  
  return 0;
}

Download this code

Download these to the same directory. Then you can create a library from the Hilbert.c and ColMajorMat.c files like this:

gcc -c ColMajorMat.c
gcc -c Hilbert.c
ar -crs libhilbert.a Hilbert.o ColMajorMat.o

You should now see the file libhilbert.a. The ar command creates archives (libraries) out of given object files. You can now compile the main program like this:

gcc Main.c -o Hilbert.ex -L. -lhilbert

Note that this is nearly the same as what we saw above when linking against libm.a. The only difference is the the presence of the flag -L.. This tells gcc to look for the library at the location ., which is the current directory. Does this remind you of anything from the #include directive?

Exercise: Try moving the files for the library into another directory and building the main program against them. This is a more typical usage pattern. What do you need to change in the compile command for the main program?

Remark: Of course, all of this can be written into a makefile.

Question: What does the flag -crs do in the ar command?

Interoperability with Fortran

The above pattern of including a header file and linking against a library file works well, and is standard, for external libraries written in C. However, if the library was written in another language, e.g. Fortran, there will not be any header file for you to include. How then can we tell the compiler what functions will be available for use later and what their signatures are?

The extern keyword denotes a function whose definition will later be filled in by the linker. Of course we still need to know the signature of each function we intend to use from the library.

Question: How do other languages annotate functions (or other symbols exported into the library) and how do we know what to call them inside our C code?

This might seem like a silly question, the name of the function is the name of the function right? It turns out that this is slightly different between each language. For instance, to call Fortran functions and subroutines from C code you need to use all lower case letters and append a trailing underscore. Name mangling in C++ can also lead to lots of complication, though we won’t talk about that here.

All of this is easier to see with an example. By far the most common use case for this is to use Lapack from C. Let’s use our above library to make an arbitrary size Hilbert matrix stored in column major format. Then using Lapack, let’s try to find the spectrum of it. Consider this simple driver code:

/* File: HilbertSpectrum.c
 * Author: Ian May
 * Purpose: Demonstrate interoperatbility with Fortran by finding the
 *          spectrum of the Hilbert matrix created in the previous
 *          example
 */

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

#include "Hilbert.h"

/* Declare dsyev as external */
/* Note the trailing underscore, and that all arguments are pointers */
extern void dsyev_(char*,char*,int*,double*,int*,double*,double*,int*,int*);

int main(int argc, char **argv)
{
  /* Matrix sizes, easily set from command line */
  int n = argc>1 ? atoi(argv[1]) : 3;

  /* Make space to store the eigenvalues */
  double evals[n];
  
  /* Create and fill a Hilbert matrix */
  ColMajorMat A = CreateHilbertMatrix(n,n);

  /* Call dsyev */
  char jobz = 'N',uplo = 'U';
  int lwork = 3*n,info = 0;
  double work[lwork];
  dsyev_(&jobz,&uplo,&n,A->data,&n,evals,work,&lwork,&info);

  /* Print out result */
  printf("evals =");
  for (int i=0; i<n; i++) {
    printf(" %e |",evals[i]);
  }
  printf("\n");

  /* Don't forget to clean up! */
  DestroyMatrix(A);
  
  return 0;
}

Download this code

Then, assuming that the previous library code is nested under a subdirectory called ColMajorLib, you can compile this example by:

gcc -I ColMajorLib HilbertSpectrum.c -o HilbertSpectrum.ex -L ColMajorLib/ -lhilbert -llapack

Then running this code should produce the outputs:

./HilbertSpectrum
evals = 2.687340e-03 | 1.223271e-01 | 1.408319e+00 |
./HilbertSpectrum 6
evals = 1.082799e-07 | 1.257076e-05 | 6.157484e-04 | 1.632152e-02 | 2.423609e-01 | 1.618900e+00 |

We should make some observations about this example code:

• The -I flag to gcc works just like the -L flag, but now for included files.

• The -L flag doesn’t interfere with the linker finding Lapack in the standard location.

• We linked Lapack just like in the Fortran chapter, that is we just used -llapack.

• The external declaration for DSYEV is kind of messy.

  • The trailing underscore is needed for the linker to understand the symbol name.

  • All arguments have to be pointers since Fortran always passes values by reference. Even scalar quantities.

• Since all arguments have to be pointers we need to create several dummy variables for the sole purpose of passing a value into dsyev.

These last points make the utilized code somewhat hard to read. It is generally easier to write a small wrapper function that behaves like C expects which calls dsyev (or whatever function) for you.

Exercise: Write this wrapper function and tidy up the call site in main.