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
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
|
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 -fdefault-real-8 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:¶
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.tgzOn Mac,
wget
doesn’t come with the general xcode installation. In this case you have two options:
- 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
- install
wget
manually (see wget-mac)Untar it using:
$ tar xvf blas-3.7.1.tgz
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.aNow 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:¶
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.tgzGo to the directory
$ cd lapack-3.7.1 $ cp make.inc.example make.incEdit
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
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 callsmake.inc
in the same directory,$HOME/packages/lapack-3.7.1/
. You may want to take a look at the makefile to learn how thismake
command compiles the LAPACK source codes.Note
Note that
-lblas
option does not look forblas.o
, butlibblas.a
in case with a statically linked library, orlibblas.so
in case with a shared library. Therefore, including-lblas
in compilation is equivalente to compilinglibblas.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:
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: