Skip to content

08. LAPACK

LAPACK, which stands for Linear Algebra PACKage, is a highly regarded software library for numerical linear algebra. It provides routines for solving systems of linear equations, linear least squares problems, eigenvalue problems, and singular value decomposition. LAPACK is widely used in scientific computing and engineering because of its efficiency, reliability, and extensive functionality.

LAPACK routines are optimized for performance through the use of the Basic Linear Algebra Subprograms (BLAS) for lower-level operations (such as matrix multiplication).

We wish to use LAPACK to solve a linear system in the form where is a known matrix, a known vector and an unknown vector. One example of a LAPACK solver which achieves this with double precision is dgbsvx. By examining the documentation, this LAPACK routine takes many inputs. Consequently, it's useful to write a specialised subroutine for LAPACK calls that have all of these extra inputs already defined.

Consider the subroutine linear_algebra.f90: solver_banded_double_precision. We begin this subroutine with

   Subroutine solver_banded_double_precision(n_input, nband, sub_diag, sup_diag, l, rhs, soln)

The documentation for dgbsvx describes what data the subroutine requires. The next section of our code maps the data the subroutine receives to these requirements. Importantly, we only want the right hand side to be a vector, not a matrix. Therefore the dimensions of B (the LAPACK right hand side) is (:,1).

The dgbsvx call is then made. This is followed by a check on the solution output. This is discussed in the code. 

 !! X is the output
      Call dgbsvx(FACT, TRANS, N, KL, KU, NRHS, AB, LDAB, AFB, LDAFB, IPIV,&
         &EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, WORK, IWORK, INFO)

We then have an if statement on the value of INFO which is returned from dgbsvx and describes whether the solution was successful. If INFO is not equal to zero, then dgbsvx has failed. Writing messages allows us to determine what has gone wrong.