External Libraries for Scientific Computing

Advantages of using external libraries

Many times, you may find it useful and convenient to install softwares which have been developed and written by others. In this way, you can simply focus more on writing kernel parts of your main algorithms, leaving others to write the rest of numerical subroutines for you. Examples of such routines may vary depending on your needs, but in general, some of the popular choices are numerical linear algebra, I/O, parallel load distributions, mesh generations, etc.

This approach not only saves your efforts and time, but also keeps your kernal algorithms computationally efficient and compatible on various platforms as long as the external libraries are continuosly 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, include and link those libraries to compile with your Fortran subroutines.

Two most popular examples of external libraries, especially for scientific computing, are the BLAS (Basic Linear Algebra Subroutines) and LAPACK (Linear Algebra PACKage) packages.

They are both freely-available softwares which allow users to call various linear algebra subroutines in their programs.

The BLAS package is originally a reference Fortran library. It is a collection of basic linear algebra routines that are standard buidling blocks for basic vector and matrix operations. With this reason, many numerical software packages (including LAPACK, LINPACK, OCTAVE, MATLAB, Mathematica, NumPy, and R) adopts the BLAS and develop their linear algebra softwares as BLAS-compatible libraries.

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

With BLAS and LAPACK, you can call routines in LAPACK from your routines to perform mathematical operations in numerical linear algebra such as (see also the LAPACK name conventions):

  • solving system of linear equations,
  • least-square problems,
  • eigenvalue problems,
  • singular value problems,
  • random number generations, etc.
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 begining: Where do you want to install and run the packages? If you are a Mac user, there is no need to install them because Mac comes with BLAS and LAPACK since Maverick. You can see an already compiled LAPACK library in:

$ cd /usr/lib/
$ ls -all
$ ...
$ 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 in the location.

To further check if you have the BLAS and LAPACK properly installed already, go to lectureNote/chapters/chapt02/codes/lapack/linear and compile solve1.f90 using the LAPACK library flag -llapack (basically you can check out this flag by compiling any Fortran code which doesn’t even call any LAPACK routines ):

$ gfortran -llapack solve1.f90

If you see a mesage such as:

ld: library not found for -llapack

then you’ll need to install the BLAS/LAPACK on your Mac similarly as described below for Linux.

On Linux, as just mentioned, you can quickly check out the availability of BLAS/LAPACK by compiling any Fortran routines using the -llapack flag:

$ gfortran -llapack dummy.f90

Now let’s consider that you find BLAS/LAPACK are not installed on your Linux machine, and you need to install them. In what follows, we consider installing both BLAS and LAPACK in a directory /packages under your home directory, for instance:

$ cd $HOME
$ mkdir packages

1. BLAS installation:

  1. You need to first install the BLAS library. You can download a tarball blas-3.7.1.tgz directly from the BLAS website or use a wget command on your terminal:

    $ cd $HOME/packages/
$ wget http://www.netlib.org/blas/blas-3.7.1.tgz

    On Mac, wget doesn’t come with the general xcode installation. In this case you have two options:

    1. use homebrew or macports as your package manager (see Items for the Class); how to use these package managers to install a software is different from wget, or
    2. install wget manually (see wget-mac)

  2. Untar it using:

    $ tar xvf blas-3.7.1.tgz

  3. Upon untarring, you will see a BLAS source file directory. Go there and compile the source codes to generate a static library:

    $ cd BLAS-3.7.1
$ gfortran -O2 -c *.f   # compiles *.f and generates object files
$ ar cr libblas.a *.o   # archives the object files into a static archive file libblas.a

Now you’re done with the BLAS installation. Let’s now install LAPACK.

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.

2. LAPACK installation:

  1. As before, you can download a tarball lapack-3.7.1.tgz from the LAPACK website or type:

    $ cd $HOME/packages
$ wget http://www.netlib.org/lapack/lapack-3.7.1.tgz
$ tar xvf lapack-3.7.1.tgz

  2. Go to the directory:

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

  3. Edit make.inc as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
####################################################################
#  LAPACK make include file.                                       #
#  LAPACK, Version 3.7.1                                           #
#  June 2017                                                       #
####################################################################

SHELL = /bin/sh

#  CC is the C compiler, normally invoked with options CFLAGS.
#
CC     = gcc
CFLAGS = -O3

#  Modify the FORTRAN and OPTS definitions to refer to the compiler
#  and desired compiler options for your machine.  NOOPT refers to
#  the compiler options desired when NO OPTIMIZATION is selected.
#
#  Note: During a regular execution, LAPACK might create NaN and Inf
#  and handle these quantities appropriately. As a consequence, one
#  should not compile LAPACK with flags such as -ffpe-trap=overflow.
#
FORTRAN = gfortran
OPTS    = -O2 -frecursive
DRVOPTS = $(OPTS)
NOOPT   = -O0 -frecursive

#  Define LOADER and LOADOPTS to refer to the loader and desired
#  load options for your machine.
#
LOADER   = gfortran
LOADOPTS =

#  The archiver and the flag(s) to use when building an archive
#  (library).  If your system has no ranlib, set RANLIB = echo.
#
ARCH      = ar
ARCHFLAGS = cr
RANLIB    = ranlib

#  Timer for the SECOND and DSECND routines
#
#  Default:  SECOND and DSECND will use a call to the
#  EXTERNAL FUNCTION ETIME
#TIMER = EXT_ETIME
#  For RS6K:  SECOND and DSECND will use a call to the
#  EXTERNAL FUNCTION ETIME_
#TIMER = EXT_ETIME_
#  For gfortran compiler:  SECOND and DSECND will use a call to the
#  INTERNAL FUNCTION ETIME
TIMER = INT_ETIME
#  If your Fortran compiler does not provide etime (like Nag Fortran
#  Compiler, etc...) SECOND and DSECND will use a call to the
#  INTERNAL FUNCTION CPU_TIME
#TIMER = INT_CPU_TIME
#  If none of these work, you can use the NONE value.
#  In that case, SECOND and DSECND will always return 0.
#TIMER = NONE

#  Uncomment the following line to include deprecated routines in
#  the LAPACK library.
#
#BUILD_DEPRECATED = Yes

#  LAPACKE has the interface to some routines from tmglib.
#  If LAPACKE_WITH_TMG is defined, add those routines to LAPACKE.
#
#LAPACKE_WITH_TMG = Yes

#  Location of the extended-precision BLAS (XBLAS) Fortran library
#  used for building and testing extended-precision routines.  The
#  relevant routines will be compiled and XBLAS will be linked only
#  if USEXBLAS is defined.
#
#USEXBLAS = Yes
#XBLASLIB = -lxblas

#  The location of the libraries to which you will link.  (The
#  machine-specific, optimized BLAS library should be used whenever
#  possible.)
#
#BLASLIB      = ../../librefblas.a                #Dongwook:  comment out the default one
BLASLIB  = -lblas -L$(HOME)/packages/BLAS-3.7.1   #Dongwook: put your full BLAS path
#BLASLIB  = $(HOME)/packages/BLAS-3.7.1/libblas.a  #Dongwook: equivalent to the above line
CBLASLIB     = ../../libcblas.a
LAPACKLIB    = liblapack.a
TMGLIB       = libtmglib.a
LAPACKELIB   = liblapacke.a

  1. 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 that your machine allows for parallel processings.

Note

There is a master Makefile which calls make.inc in the same directory, $HOME/packages/lapack-3.7.1/. You may want to take a look at the makefile to learn how this make command compiles the LAPACK source codes.

Note

Note that -lblas option does not look for blas.o, but libblas.a in case with a statically linked library, or libblas.so in case with a shared library. Therefore, including -lblas in compilation is equivalente to compiling libblas.a.

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, say, in your .bashrc or .bash_profile. If it’s not the case, however, you need to add -L flags followed by the directory path, in the makefile if you’ve put libraries in nonstandard places.

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

3. Testing:

With these you should be able to run the codes in lectureNote/chapters/chapt02/codes/lapack/linear:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
!! /codes/lapack/solve1.f90
!! This routine solves Ax=b using DGESV routine in LAPACK.
!!
!! Recall that the naming convention follows (see http://www.netlib.org/lapack/lug/node26.html)
!!   D:  double precision
!!   GE: general type of matrix
!!   SV: simple driver by factoring A and overwriting B with the solution x
!!
!! See also:
!! a. https://software.intel.com/sites/products/documentation/doclib/mkl_sa/11/mkl_lapack_examples/dgesv.htm
!! b. http://www.netlib.org/clapack/old/double/dgesv.c
!!


program solve1
  implicit none
  integer, parameter :: n = 3
  real, dimension(n) :: x,b
  real, dimension(n,n) :: a
  integer :: i, info, lda, ldb, nrhs
  integer, dimension(n) :: ipiv

  a = reshape((/1., 2., 2., 0., -4., -6., 0., 0., -1./), (/n,n/))
  b = (/3.,-6.,1./)
  x = b

  nrhs = 1
  lda = n
  ldb = n

  call dgesv(n, nrhs, a, lda, ipiv, x, ldb, info)

  print*, 'solving Ax = b'
  do i = 1,n
     print *, i, x(i)
  end do
  
end program solve1

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# codes/lapack/random/makefile

FC = gfortran
FFLAGS = -O2  -fdefault-real-8 -fdefault-double-8 
LFLAGS = -lblas -llapack
.PHONY: clean test1

%.o : %.f90
	$(FC) $(FFLAGS) -c $<

# this is correct
solve1.exe: solve1.o
	$(FC) solve1.o $(LFLAGS) -o solve1.exe

# this can be a bug
#solve1.exe: solve1.o
#	$(FC) $(LFLAGS) solve1.o  -o solve1.exe

test1: solve1.exe
	./solve1.exe


clean:
	rm -f *.o *.exe

and the codes in lectureNote/chapters/chapt02/codes/lapack/random. Note in the above makefile, it is important to include $(LFLAGS) after solve1.o so that solve1.o can refer to the libraries.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
! Initialize the random number generator using current time,
! so a new sequence of random numbers is generated each 
! execution time.

! Taken from http://gcc.gnu.org/onlinedocs/gfortran/RANDOM_005fSEED.html

SUBROUTINE init_random_seed()
    INTEGER :: i, n, clock
    INTEGER, DIMENSION(:), ALLOCATABLE :: seed
  
    CALL RANDOM_SEED(size = n)
    ALLOCATE(seed(n))
  
    CALL SYSTEM_CLOCK(COUNT=clock)
  
    seed = clock + 37 * (/ (i - 1, i = 1, n) /)
    CALL RANDOM_SEED(PUT = seed)
  
    print *, "Using random seed = ", seed
    print *, " "

    DEALLOCATE(seed)
END SUBROUTINE

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
! /codes/lapack/random/randomsys1.f90

program randomsys1
    implicit none
    integer, parameter :: nmax=1000
    real(kind=8), dimension(nmax) :: b, x
    real(kind=8), dimension(nmax,nmax) :: a
    real(kind=8) :: err
    integer :: i, info, lda, ldb, nrhs, n
    integer, dimension(nmax) :: ipiv

    ! initialize random number generator seed
    ! if you remove this, the same numbers will be generated each
    ! time you run this code.
    call init_random_seed()  

    print *, "Input n ... "
    read *, n
    if (n<1 .or. n>nmax) then
        print *, "n must be positive and no greater than ",nmax
        stop
        endif
    call random_number(a(1:n,1:n))
    call random_number(x(1:n))
    b(1:n) = matmul(a(1:n,1:n),x(1:n)) ! compute RHS
    
    nrhs = 1 ! number of right hand sides in b
    lda = nmax  ! leading dimension of a
    ldb = nmax  ! leading dimension of b

    call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)

    ! Note: the solution is returned in b
    ! and a has been changed.

    ! compare computed solution to original x:
    print *, "         x          computed       rel. error"
    do i=1,n
        err = abs(x(i)-b(i))/abs(x(i))
        print '(3d16.6)', x(i),b(i),err
        enddo

end program randomsys1

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
! /codes/lapack/random/randomsys2.f90

program randomsys2
    implicit none
    real(kind=8), dimension(:),allocatable :: x,b
    real(kind=8), dimension(:,:),allocatable :: a
    real(kind=8) :: err
    integer :: i, info, lda, ldb, nrhs, n
    integer, dimension(:), allocatable :: ipiv

    ! initialize random number generator seed
    ! if you remove this, the same numbers will be generated each
    ! time you run this code.
    call init_random_seed()  

    print *, "Input n ... "
    read *, n

    allocate(a(n,n))
    allocate(b(n))
    allocate(x(n))
    allocate(ipiv(n))

    call random_number(a)
    call random_number(x)
    b = matmul(a,x) ! compute RHS
    
    nrhs = 1 ! number of right hand sides in b
    lda = n  ! leading dimension of a
    ldb = n  ! leading dimension of b

    call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)

    ! Note: the solution is returned in b
    ! and a has been changed.

    ! compare computed solution to original x:
    print *, "         x          computed       rel. error"
    do i=1,n
        err = abs(x(i)-b(i))/abs(x(i))
        print '(3d16.6)', x(i),b(i),err
        enddo

    deallocate(a,b,ipiv)

end program randomsys2

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
! /codes/lapack/random/randomsys3.f90

program randomsys3
    implicit none
    real(kind=8), dimension(:),allocatable :: x,b,work
    real(kind=8), dimension(:,:),allocatable :: a
    real(kind=8) :: errnorm, xnorm, rcond, anorm, colsum
    integer :: i, info, lda, ldb, nrhs, n, j
    integer, dimension(:), allocatable :: ipiv
    integer, allocatable, dimension(:) :: iwork
    character, dimension(1) :: norm

    ! initialize random number generator seed
    ! if you remove this, the same numbers will be generated each
    ! time you run this code.
    call init_random_seed()  

    print *, "Input n ... "
    read *, n

    allocate(a(n,n))
    allocate(b(n))
    allocate(x(n))
    allocate(ipiv(n))

    call random_number(a)
    call random_number(x)

    b = matmul(a,x)    ! compute RHS

    ! compute 1-norm needed for condition number

    anorm = 0.d0
    do j=1,n
        colsum = 0.d0
        do i=1,n
            colsum = colsum + abs(a(i,j))
            enddo
        anorm = max(anorm, colsum)
        enddo

    
    nrhs = 1 ! number of right hand sides in b
    lda = n  ! leading dimension of a
    ldb = n  ! leading dimension of b

    call dgesv(n, nrhs, a, lda, ipiv, b, ldb, info)

    ! compute 1-norm of error
    errnorm = 0.d0
    xnorm = 0.d0
    do i=1,n
        errnorm = errnorm + abs(x(i)-b(i))
        xnorm = xnorm + abs(x(i))
        enddo

    ! relative error in 1-norm:
    errnorm = errnorm / xnorm


    ! compute condition number of matrix:
    ! note: uses A returned from dgesv with L,U factors:

    allocate(work(4*n))
    allocate(iwork(n))
    norm = '1'  ! use 1-norm
    call dgecon(norm,n,a,lda,anorm,rcond,work,iwork,info)

    if (info /= 0) then
        print *, "*** Error in dgecon: info = ",info
        endif

    print 201, n, 1.d0/rcond, errnorm
201 format("For n = ",i4," the approx. condition number is ",e10.3,/, &
           " and the relative error in 1-norm is ",e10.3)    

    deallocate(a,b,ipiv)
    deallocate(work,iwork)

end program randomsys3

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
# /codes/lapack/random/makefile

FC = gfortran
FFLAGS = -O3
LFLAGS = -lblas -llapack
.PHONY: clean test1 test2 test3

%.o : %.f90
	$(FC) $(FFLAGS) -c $<

randomsys1.exe: randomsys1.o init_random_seed.o
	$(FC) randomsys1.o init_random_seed.o $(LFLAGS) -o randomsys1.exe

test1: randomsys1.exe
	./randomsys1.exe

randomsys2.exe: randomsys2.o init_random_seed.o
	$(FC) randomsys2.o init_random_seed.o $(LFLAGS) -o randomsys2.exe

test2: randomsys2.exe
	./randomsys2.exe

randomsys3.exe: randomsys3.o  init_random_seed.o
	$(FC) randomsys3.o init_random_seed.o $(LFLAGS) -o randomsys3.exe

test3: randomsys3.exe
	./randomsys3.exe

clean:
	rm -f *.o *.exe

Some more references: