External Libraries for Scientific Computing

Advantages of using external libraries

Very rarely will you write software in total isolation, and frequently you will want to rely on existing implementations for some components of your code. This lets you focus on the overall purpose of the code while letting others write and maintain some general purpose routines. Common examples are numerical linear algebra routines, (binary format) I/O, parallel load distribution and communication, etc.

This approach not only saves you effort and time, but also keeps your code computationally efficient and compatible across various platforms (as long as the external libraries are continuously supported and maintained).

About BLAS and LAPACK

In this section we learn how to install external software packages. We are going to search for them, compile them to produce libraries, then include and link those libraries so that you may call them from your code.

Two of the most popular examples of external libraries, especially for scientific computing, are BLAS (Basic Linear Algebra Subroutines) and LAPACK (Linear Algebra PACKage). Linear algebra is a fundamental building block of numerical methods, and hence these libraries are fundamental building blocks for many scientific codes.

They are both freely-available and allow users to call various linear algebra subroutines in their programs. Indeed, linear algebra powers a great deal of scientific computing.

The BLAS package is a reference Fortran library. It is a collection of fundamental linear algebra routines that comprise standard building blocks for vector and matrix operations. These basic operations have been heavily optimized and are incredibly fast. For this reason, many numerical software packages (including LAPACK, LINPACK, OCTAVE, MATLAB, Mathematica, NumPy, and R) adopt BLAS and develop their linear algebra software as BLAS-compatible libraries.

LAPACK is an improved version of LINPACK, which provides standard numerical linear algebra routines using BLAS as a building block for its implementation. That is to say, routines in LAPACK perform computations by calling the routines from BLAS.

With BLAS and LAPACK installed you will have access to all of the implemented linear algebra routines from your own code, including:

  • solving systems of linear equations,

  • solving least squares problems,

  • eigenvalue decompositions,

  • singular value decompositions,

  • random number generation

  • and more…

The naming convention for the subroutines is a little old-school. The LAPACK name conventions) page has more information on the topic.

See also:

Installing BLAS and LAPACK

Since LAPACK routines call the BLAS routines, you have to install the BLAS library first before installing LAPACK.

There is one important question to ask before beginning: Where do you want to install and run the packages? Do you want to install them system-wide? For only one user?

If you are a Linux or MacOS user then there is chance that these libraries are already present on your machine. Otherwise they can easily be obtained from your system’s package manager, which might look like:

$ # On Ubuntu or Debian
$ sudo apt install libblas-dev liblapack-dev
$ # Or on Arch derivatives
$ sudo pacman -S lapack blas

Remark You can often install the LAPACK documentation directly to your machine. This will provide man-pages for the LAPACK routines which can be quite helpful

If you are a Mac user, then it is quite likely that BLAS and LAPACK are already present on your machine. Try looking in the /usr/lib directory for entries containing lapack in their name:

$ cd /usr/lib/
$ ls -al
$ ...
$ liblapack.dylib -> ../../System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libLAPACK.dylib

The -> means that the default provided library file liblapack.dylib in /usr/lib/ is a symbolic link to the file called libLAPACK.dylib buried under that long path.

To further check if you have the BLAS and LAPACK properly installed already, download the code below:

 1!! /codes/lapack/solve1.f90
 2!! This routine solves Ax=b using DGESV routine in LAPACK.
 3!!
 4!! Recall that the naming convention follows (see http://www.netlib.org/lapack/lug/node26.html)
 5!!   D:  double precision
 6!!   GE: general type of matrix
 7!!   SV: simple driver by factoring A and overwriting B with the solution x
 8!!
 9!! See also:
10!! a. https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgesv.htm
11!! b. http://www.netlib.org/clapack/old/double/dgesv.c
12!!
13
14
15program solve1
16  implicit none
17  integer, parameter :: n=3, fp = selected_real_kind(15)
18  real (fp), dimension(n) :: x,b
19  real (fp), dimension(n,n) :: a
20  integer :: i, info, lda, ldb, nrhs
21  integer, dimension(n) :: ipiv
22
23  a = reshape((/1._fp, 2._fp, 2._fp, 0._fp, -4._fp, -6._fp, 0._fp, 0._fp, -1._fp/), (/n,n/))
24  b = (/3._fp,-6._fp,1._fp/)
25  x = b
26
27  nrhs = 1
28  lda = n
29  ldb = n
30
31  call dgesv(n, nrhs, a, lda, ipiv, x, ldb, info)
32
33  print*, 'solving Ax = b'
34  do i = 1,n
35    print *, i, x(i)
36  end do
37
38end program solve1

Download this code

and compile solve1.f90 using the LAPACK library flag -llapack (actually you can check out this flag by compiling any Fortran code which doesn’t even call any LAPACK routines)

$ gfortran solve1.f90 -o solve1.ex -llapack

If you see a message such as

ld: library not found for -llapack

then you’ll need to install the BLAS/LAPACK on your machine.

Note

I strongly recommend installing BLAS and LAPACK through your package manager.

If you can not install BLAS and LAPACK from your package manager, then proceed below for installing from source. The ubiquity of BLAS and LAPACK make them commonly available, however other libraries that you may want to use will only be available as source. The instructions below for LAPACK illustrate a fairly common procedure for this type of installation.

BLAS installation from source:

Note

The current Lapack source distribution includes BLAS. This step is retained in case you need it, but most people should skip ahead to the next section.

  1. You need to first install the BLAS library. You can download a tarball blas-3.10.0.tgz directly from the BLAS website.

  2. Untar it using:

    $ tar -xvf blas-3.10.0.tgz

  3. Upon untarring, you will see a BLAS source file directory. cd into it and verify that the contents of make.inc match your system. Then compile the full library by calling

    $ make all

You have now created a static BLAS library that you can link to your code.

Note

If you wish to learn more about the command ar on Linux/Unix, please check out the following articles: article-ar-fortran, article-ar-C, lib-name.a convention.

Installing LAPACK from source:

  1. As before, you can download a tarball lapack-3.10.0.tgz from the LAPACK website.

    $ cd $HOME/packages
$ tar -xvzf lapack-3.10.0.tar.gz

  2. Go to the directory

    $ cd lapack-3.10.0
$ cp make.inc.example make.inc

  3. Edit make.inc if needed. Here you can specify what compiler to use, what flags to use, and whether or not the included BLAS library should also be built.

  4. You’re now ready to compile LAPACK by typing

    $ make

    or you can run make in a parallel mode

    $ make -j <N>

    where N = 2, 4, or any number up to the number of threads that your machine has.

Note

There is a primary Makefile which calls make.inc in the same directory. You may want to take a look at the makefile to see how it works internally.

Note

Note that -lblas option does not look for blas.o, but libblas.a for a statically linked version of the library, or libblas.so for a shared library.

Note

As just mentioned, the correct way to link a static library (i.e., *.a files) is using -l, but that only works if the library can be found on the search path, $PATH. On the other hand, you can direct the compiler to correct location by using the -L flag followed by the directory.

On the other hand, if you’re interested in linking a shared library (i.e., *.so files), you define an environment variable, LD_LIBRARY_PATH which include the directory paths to the target shared library (see more howto-sharedLib). Also see an article on their differences article-shared-static.

Note

Also the order matters. The linker processes its input files in the order in which they appear on the command line – left to right. See article.

Also read:

  1. article 1

  2. article 2

  3. article 3

Testing your installation

With these you should be able to run the codes below:

 1!! /codes/lapack/solve1.f90
 2!! This routine solves Ax=b using DGESV routine in LAPACK.
 3!!
 4!! Recall that the naming convention follows (see http://www.netlib.org/lapack/lug/node26.html)
 5!!   D:  double precision
 6!!   GE: general type of matrix
 7!!   SV: simple driver by factoring A and overwriting B with the solution x
 8!!
 9!! See also:
10!! a. https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgesv.htm
11!! b. http://www.netlib.org/clapack/old/double/dgesv.c
12!!
13
14
15program solve1
16  implicit none
17  integer, parameter :: n=3, fp = selected_real_kind(15)
18  real (fp), dimension(n) :: x,b
19  real (fp), dimension(n,n) :: a
20  integer :: i, info, lda, ldb, nrhs
21  integer, dimension(n) :: ipiv
22
23  a = reshape((/1._fp, 2._fp, 2._fp, 0._fp, -4._fp, -6._fp, 0._fp, 0._fp, -1._fp/), (/n,n/))
24  b = (/3._fp,-6._fp,1._fp/)
25  x = b
26
27  nrhs = 1
28  lda = n
29  ldb = n
30
31  call dgesv(n, nrhs, a, lda, ipiv, x, ldb, info)
32
33  print*, 'solving Ax = b'
34  do i = 1,n
35    print *, i, x(i)
36  end do
37
38end program solve1

Download this code

 1# codes/lapack/linear/makefile
 2
 3FC = gfortran
 4FFLAGS = -O2 -g -Wall -Wextra
 5LFLAGS = -llapack
 6#LFLAGS = -L ../lapack-3.10.0 -llapack
 7
 8.PHONY: clean all
 9
10all: solve1.ex randomsys1.ex randomsys2.ex randomsys3.ex
11
12# this is correct
13solve1.ex: solve1.o
14	$(FC) solve1.o $(LFLAGS) -o solve1.ex
15
16# this can be a bug depending on your gfortran version and your platform
17# solve1.ex: solve1.o
18# 	$(FC) $(LFLAGS) solve1.o  -o solve1.ex
19
20randomsys1.ex: randomsys1.o init_random_seed.o
21	$(FC) randomsys1.o init_random_seed.o $(LFLAGS) -o randomsys1.ex
22
23randomsys2.ex: randomsys2.o init_random_seed.o
24	$(FC) randomsys2.o init_random_seed.o $(LFLAGS) -o randomsys2.ex
25
26randomsys3.ex: randomsys3.o  init_random_seed.o
27	$(FC) randomsys3.o init_random_seed.o $(LFLAGS) -o randomsys3.ex
28
29%.o : %.f90
30	$(FC) $(FFLAGS) -c $<
31
32clean:
33	rm -f *.o *.ex

Download this code

Note in the above Makefile, it is important to include $(LFLAGS) after solve1.o so that solve1.o can refer to the libraries. Not doing so can cause issues on some systems.

The following programs construct and solve random linear systems of various sizes.

 1! Initialize the random number generator using current time,
 2! so a new sequence of random numbers is generated each 
 3! execution time.
 4
 5! Taken from http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html
 6
 7SUBROUTINE init_random_seed()
 8  INTEGER :: i, n, clock
 9  INTEGER, DIMENSION(:), ALLOCATABLE :: seed
10
11  CALL RANDOM_SEED(size = n)
12  ALLOCATE(seed(n))
13
14  CALL SYSTEM_CLOCK(COUNT=clock)
15
16  seed = clock + 37 * (/ (i - 1, i = 1, n) /)
17  CALL RANDOM_SEED(PUT = seed)
18
19  print *, "Using random seed = ", seed
20  print *, " "
21
22  DEALLOCATE(seed)
23END SUBROUTINE init_random_seed

Download this code

 1! /codes/lapack/random/randomsys1.f90
 2
 3program randomsys1
 4  implicit none
 5  integer, parameter :: nmax=1000, fp = selected_real_kind(15)
 6  real (fp), dimension(nmax) :: b, x
 7  real (fp), dimension(nmax,nmax) :: a
 8  real (fp) :: err
 9  integer :: i, info, lda, ldb, nrhs, n
10  integer, dimension(nmax) :: ipiv
11
12  ! initialize random number generator seed
13  ! if you remove this, the same numbers will be generated each
14  ! time you run this code.
15  call init_random_seed()  
16
17  print *, "Input n ... "
18  read *, n
19  if (n<1 .or. n>nmax) then
20    print *, "n must be positive and no greater than ",nmax
21    stop
22  endif
23  call random_number(a(1:n,1:n))
24  call random_number(x(1:n))
25  b(1:n) = matmul(a(1:n,1:n),x(1:n)) ! compute RHS
26
27  nrhs = 1 ! number of right hand sides in b
28  lda = nmax  ! leading dimension of a
29  ldb = nmax  ! leading dimension of b
30
31  call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
32
33  ! Note: the solution is returned in b
34  ! and a has been changed.
35
36  ! compare computed solution to original x:
37  print *, "          x       computed     rel. error"
38  do i=1,n
39    err = abs(x(i)-b(i))/abs(x(i))
40    print "(ES12.6,' | ',ES12.6,' | ',ES12.6)", x(i),b(i),err
41  end do
42
43end program randomsys1

Download this code

 1! /codes/lapack/random/randomsys2.f90
 2
 3program randomsys2
 4  implicit none
 5  integer, parameter :: fp = selected_real_kind(15)
 6  real (fp), dimension(:),allocatable :: x,b
 7  real (fp), dimension(:,:),allocatable :: a
 8  real (fp) :: err
 9  integer :: i, info, lda, ldb, nrhs, n
10  integer, dimension(:), allocatable :: ipiv
11
12  ! initialize random number generator seed
13  ! if you remove this, the same numbers will be generated each
14  ! time you run this code.
15  call init_random_seed()  
16
17  print *, "Input n ... "
18  read *, n
19
20  allocate(a(n,n))
21  allocate(b(n))
22  allocate(x(n))
23  allocate(ipiv(n))
24
25  call random_number(a)
26  call random_number(x)
27  b = matmul(a,x) ! compute RHS
28
29  nrhs = 1 ! number of right hand sides in b
30  lda = n  ! leading dimension of a
31  ldb = n  ! leading dimension of b
32
33  call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
34
35  ! Note: the solution is returned in b
36  ! and a has been changed.
37
38  ! compare computed solution to original x:
39  print *, "          x       computed     rel. error"
40  do i=1,n
41    err = abs(x(i)-b(i))/abs(x(i))
42    print "(ES12.6,' | ',ES12.6,' | ',ES12.6)", x(i),b(i),err
43  end do
44
45  deallocate(a,b,ipiv)
46
47end program randomsys2

Download this code

 1! /codes/lapack/random/randomsys3.f90
 2
 3program randomsys3
 4  implicit none
 5  integer, parameter :: fp = selected_real_kind(15)
 6  real (fp), dimension(:),allocatable :: x,b,work
 7  real (fp), dimension(:,:),allocatable :: a
 8  real (fp) :: errnorm, xnorm, rcond, anorm, colsum
 9  integer :: i, info, lda, ldb, nrhs, n, j
10  integer, dimension(:), allocatable :: ipiv
11  integer, allocatable, dimension(:) :: iwork
12  character, dimension(1) :: norm
13
14  ! initialize random number generator seed
15  ! if you remove this, the same numbers will be generated each
16  ! time you run this code.
17  call init_random_seed()  
18
19  print *, "Input n ... "
20  read *, n
21
22  allocate(a(n,n))
23  allocate(b(n))
24  allocate(x(n))
25  allocate(ipiv(n))
26
27  call random_number(a)
28  call random_number(x)
29
30  b = matmul(a,x)    ! compute RHS
31
32  ! compute 1-norm needed for condition number
33
34  anorm = 0.d0
35  do j=1,n
36    colsum = 0.d0
37    do i=1,n
38      colsum = colsum + abs(a(i,j))
39    enddo
40    anorm = max(anorm, colsum)
41  enddo
42
43
44  nrhs = 1 ! number of right hand sides in b
45  lda = n  ! leading dimension of a
46  ldb = n  ! leading dimension of b
47
48  call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
49
50  ! compute 1-norm of error
51  errnorm = 0.d0
52  xnorm = 0.d0
53  do i=1,n
54    errnorm = errnorm + abs(x(i)-b(i))
55    xnorm = xnorm + abs(x(i))
56  enddo
57
58  ! relative error in 1-norm:
59  errnorm = errnorm / xnorm
60
61
62  ! compute condition number of matrix:
63  ! note: uses A returned from dgesv with L,U factors:
64
65  allocate(work(4*n))
66  allocate(iwork(n))
67  norm = '1'  ! use 1-norm
68  call dgecon(norm,n,a,lda,anorm,rcond,work,iwork,info)
69
70  if (info /= 0) then
71    print *, "*** Error in dgecon: info = ",info
72  endif
73
74  print 201, n, 1.d0/rcond, errnorm
75201 format("For n = ",i4," the approx. condition number is ",ES10.3,/,&
76      " and the relative error in 1-norm is ",ES10.3)    
77
78  deallocate(a,b,ipiv)
79  deallocate(work,iwork)
80
81end program randomsys3

Download this code

Some more references: