LAPACK

LAPACK (Linear Algebra PACKage) is a BSD-licenced FORTRAN 77 library that provides routines for solving the most commonly occurring problems in numerical linear algebra, like:

LAPACK supersedes the older EISPACK and LINPACK libraries. Packages are available for all Unix-like operating systems. On FreeBSD, install LAPACK with:

# pkg install math/lapack

Solving Linear Equations

The example program solves the following system of linear equations with LAPACK:

A * x = B

The LAPACK subroutine sgesv() computes the solution to a real system of linear equations A · x = B, where A is an n-by-n matrix and x and B are n-by-nrhs matrices.

call sgesv(n, nrhs, a, lda, ipiv, b, ldb, info)
Argument Type Description
n integer The number of linear equations, i. e., the order of the matrix A (n ≥ 0).
nrhs integer The number of right hand sides, i. e., the number of columns of the matrix B (nrhs ≥ 0).
a real, dimension(lda, *) Array of type real of size lda-by-n.
lda integer The leading dimension of array A (lda ≥ max(1, n)).
ipiv integer, dimension(*) The pivot indices that define the permutation matrix.
b real, dimension(ldb, *) On entry, the n-by-nrhs matrix of right hand side matrix b. On exit, if info = 0, the n-by-nrhs solution matrix.
ldb integer The leading dimension of array B (ldb ≥ max(1, n)).
info integer = 0: successful exit;
< 0: if info = −i, the i-th argument has an illegal value;
> 0: if info = i, U(i, i) is zero. The factorisation has been completed, but the factor U is exactly singular, so the solution could not be computed.

The LAPACK routines have to be imported using the external statement:

! example.f90
program main
    ! Example program that solves the matrix equation A * x = B using LAPACK.
    implicit none
    external :: sgesv

    real    :: a(2, 2)  ! Matrix A.
    real    :: b(2)     ! Matrix B.
    real    :: pivot(2) ! Pivot indices (list of swap operations).
    integer :: rc       ! Return code.

    a = reshape([ 2., 3., 1., 1. ], [ 2, 2 ])
    b = [ 5., 6. ]

    call sgesv(2, 1, a, 2, pivot, b, 2, rc)

    if (rc == 0) then
        print '(a, 2(f0.4, ", "))', 'Solution (x, y): ', b
    else
        print '(a, i0)', 'Error: ', rc
    end if
end program main

Compile, link, and run the example with:

$ gfortran9 -L/usr/local/lib/ -o example example.f90 -llapack -lblas
$ ./example
Solution (x, y): 1.0000, 3.0000

If LAPACK and BLAS are compiled as static libraries, point the library search path -L to their actual location.

LAPACK95

LAPACK95 is a Fortran 95 interface to the FORTRAN 77 LAPACK library. It brings several improvements over the original LAPACK package, by using language features introduced in Fortran 95, like assumed-shape arrays, optional arguments, and generic interfaces. The underlying LAPACK code remains unchanged, as LAPACK95 provides only interfaces to existing routines. The users’ guide and an index of all routines are available online.

Installation

LAPACK95 has been ported to most Unix-like operating systems. On FreeBSD, install the package with:

# pkg install math/lapack95

If you choose to build LAPACK95 from FreeBSD Ports, you can either link the BLAS reference implementation or the optimised ATLAS library. But LAPACK95 can also be compiled simply from source:

$ fetch https://www.netlib.org/lapack95/lapack95.tgz
$ tar xfvz lapack95.tgz
$ cd LAPACK95/

Edit make.inc and set your preferred compiler and flags, for example:

FC    = gfortran9 -ffree-form
FC1   = gfortran9 -ffixed-form
OPTS0 =

Then, compile LAPACK95 with make:

$ cd SRC/
$ make single_double_complex_dcomplex

This outputs the static library lapack95.a and the module files for linking inside directory lapack95_modules/.

Solving Linear Equations

In comparison, calling LAPACK95 routines requires less boilerplate code than LAPACK itself. For example, we can call la_gesv() to solve any real or complex linear system A · x = B.

call la_gesv(a, b[, ipiv][, info])
Argument Type Description
a type(wp), dimension(:, :) On entry, matrix A, an array of type wp (real or complex). On Exit, the factors from the factorisation.
b type(wp), dimension(:, :) On entry, the matrix B. On exit, if info = 0, the solution matrix.
ipiv integer, dimension(:) The pivot indices that define the permutation matrix (optional).
info integer = 0: successful exit;
< 0: if info = −i, the i-th argument has an illegal value;
> 0: if info = i, U(i, i) is zero. The factorisation has been completed, but the factor U is exactly singular, so the solution could not be computed. If info is not present and an error occurs, then the program is terminated with an error message (optional).

The LAPACK95 example solves the linear equation:

A * x = B

We set the working precision wp to the double precision dp provided by LAPACK95.

! example95.f90
program main
    ! Example program that solves the matrix equations A * x = B using LAPACK95.
    use :: la_precision, only: wp => dp
    use :: f95_lapack,   only: la_gesv
    implicit none
    real(kind=wp) :: a(2, 2) = reshape([ 2., 3., 1., 1. ], [ 2, 2 ])
    real(kind=wp) :: b(2)    = [ 5., 6. ]

    call la_gesv(a, b)

    print '(a, 2(f0.4, ", "))', 'Solution (x, y): ', b
end program main

LAPACK95 is linked with -llapack95 -llapack -lblas -ltmglib:

$ gfortran9 -L/usr/local/lib/ -I/usr/local/lib/lapack95_modules/ \
  -o example95 example95.f90 -llapack95 -llapack -lblas -ltmglib
$ ./example95
Solution (x, y): 1.0000, 3.0000

On Linux, change the prefix to /usr instead of /usr/local.

Solving Linear Least Squares Problems

The two-dimensional conformal coordinate transformation is a more complex problem that can be solved by a system of normal equations. The normal equations of a least square adjustment can be expressed as

least square

with the coefficient (or Jacobian) matrix A, the unknown transformation parameters in vector X, the common points in source system in vector L, and the residuals in vector V:

matrices

We may apply the normal matrix N = ATA to get the solution to the normal equations:

least square

But LAPACK95 lets use solve the normal equations directly. The required transformation parameters include the translation factors TX and TY, the rotation angle

rotation angle

and the scale factor

scale factor

We transform additional coordinates given in the source system into the target system by applying the parameters:

transformation

The following example is taken from the book Adjustment Computations by Ch. D. Ghilani. Given are the coordinates A, B, and C both in source (x, y) and target system (E, N). At least three common points are required to solve the normal equations. The additional points 1, 2, 3, and 4 in the source system will then be transformed into the target system.

PointxyEN
A121.622 –128.0661,049,422.4051,089.20
B141.228 187.7181,049,413.9549,659.30
C175.802 135.7281,049,244.9549,884.95
1174.148 –120.262
2513.520 –192.130
3754.444 –67.706
4972.788 120.994

We can call the LAPACK95 routine la_gels() to compute the minimum-norm least squares solution to one or more real or complex linear systems of the form A · x = b, AT · x = b or AH · x = b using a QR or LQ factorisation of A. The routine does not calculate any residuals.

call la_gels(a, b[, trans][, info])
Argument Type Description
a type(wp), dimension(:, :) On entry, matrix A, an rectangular array of shape (:, :) (real or complex). On Exit, the details of its QR or LQ factorisation.
b type(wp), dimension(:, :) On entry, the matrix B of shape (:, :). On exit, if info = 0, the solution matrix.
trans character(len=1) Specifies the form of the system of equations, either N (default) for no transpose, T for transpose, or C for conjugate transpose (optional).
info integer = 0: successful exit;
< 0: if info = −i, the i-th argument has an illegal value (optional).

The source code of the example:

! adjust.f90
program main
    use :: la_precision, only: wp => dp
    use :: f95_lapack
    implicit none
    real(kind=wp), parameter :: PI = acos(-1.0_wp)  ! Pi.
    integer,       parameter :: n  = 3              ! Number of common points.

    real(kind=wp) :: a(n * 2, 4)
    real(kind=wp) :: c(n * 2)
    real(kind=wp) :: l(n * 2)
    real(kind=wp) :: m(8), t(8)
    real(kind=wp) :: s, w
    integer       :: i, j

    print '(/, a)', repeat('-', 64)
    print '(a)',    '-- Two-Dimensional Conformal Coordinate Transformation'
    print '(a)',    '-- (Four-Parameter Helmert Transformation)'
    print '(a)',    repeat('-', 64)

    ! Common points (E, N) in the target system.
    l = [ 1049422.40_wp, 51089.20_wp, &
          1049413.95_wp, 49659.30_wp, &
          1049244.95_wp, 49884.95_wp ]

    ! Common points (x, y) in the source system.
    c = [ 121.622_wp, -128.066_wp, &
          141.228_wp,  187.718_wp, &
          175.802_wp,  135.728_wp ]

    ! Print L vector.
    print '(/, a)', 'L Vector'
    print '(a)',    repeat('-', 64)

    do i = 1, size(l)
        write (*, '(f12.3)') l(i)
    end do

    ! Fill coefficient matrix A.
    print '(/, a)', 'A Matrix'
    print '(a)',    repeat('-', 64)

    do i = 1, (n * 2) - 1, 2
        a(    i, :) = [     c(i), -c(i + 1), 1.0_wp, 0.0_wp ]
        a(i + 1, :) = [ c(i + 1),      c(i), 0.0_wp, 1.0_wp ]
    end do

    ! Print A matrix.
    do i = 1, size(a(:, 1))
        do j = 1, size(a(1, :))
            write (*, '(f12.3, " ")', advance='no') a(i, j)
        end do

        write (*, '(a)')
    end do

    ! Compute minimum-norm least squares solution to the linear system
    ! A * X = L using LAPACK95. The vector `l` contains the solution
    ! afterwards.
    call la_gels(a, l)

    ! Print solution.
    print '(/, a)', 'Transformation Parameters'
    print '(a)',    repeat('-', 64)

    print '(a, f12.4)', ' a = ', l(1)
    print '(a, f12.4)', ' b = ', l(2)
    print '(a, f12.4)', 'TX = ', l(3)
    print '(a, f12.4)', 'TY = ', l(4)

    ! Compute rotation angle `w`.
    w = atan2(l(2), l(1))

    ! Add 360° if angle is negative.
    if (w < 0) w = w + (2 * PI)

    ! Compute scale factor `s`.
    s = l(1) / cos(w)

    print '(a, f12.4, a)', ' w = ', w * (180_wp / PI), ' deg'
    print '(a, f12.4)',    ' s = ', s

    ! Transform coordinates from source to target system.
    m = [ 174.148_wp, -120.262_wp, &
          513.520_wp, -192.130_wp, &
          754.444_wp,  -67.706_wp, &
          972.788_wp,  120.994_wp ]

    do i = 1, size(m), 2
        t(i)     = (s * cos(w)) * m(i) - (s * sin(w)) * m(i + 1) + l(3)
        t(i + 1) = (s * sin(w)) * m(i) + (s * cos(w)) * m(i + 1) + l(4)
    end do

    ! Print transformed coordinates.
    print '(/, a)', 'Transformed Coordinates'
    print '(a)',    repeat('-', 64)

    do i = 1, size(t), 2
        write (*, '(2(f12.3, " "))') t(i), t(i + 1)
    end do
end program main

Compile the program with:

$ gfortran9 -L/usr/local/lib/ -I/usr/local/lib/lapack95_modules/ \
  -o adjust adjust.f90 -llapack95 -llapack -lblas -ltmglib

Then, run:

$ ./adjust
----------------------------------------------------------------
-- Two-Dimensional Conformal Coordinate Transformation
-- (Four-Parameter Helmert Transformation)
----------------------------------------------------------------

A Matrix
----------------------------------------------------------------
    121.622     128.066       1.000       0.000
   -128.066     121.622       0.000       1.000
    141.228    -187.718       1.000       0.000
    187.718     141.228       0.000       1.000
    175.802    -135.728       1.000       0.000
    135.728     175.802       0.000       1.000

L Vector
----------------------------------------------------------------
 1049422.400
   51089.200
 1049413.950
   49659.300
 1049244.950
   49884.950

Transformation Parameters
----------------------------------------------------------------
 a =      -4.5125
 b =      -0.2537
TX = 1050003.7145
TY =   50542.1311
 w =     183.2181 deg
 s =       4.5196

Transformed Coordinates
----------------------------------------------------------------
 1049187.361    51040.629
 1047637.713    51278.829
 1046582.113    50656.241
 1045644.713    49749.336

References