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, such as:

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

Link your Fortran application against LAPACK with -llapack -lblas.

Solving Linear Equations

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

AX = B

The LAPACK subroutine sgesv() computes the solution to a real system of linear equations AX = 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.

LAPACK routines have to be imported individually using the external statement if implicit none (type, external) is set (Fortran 2018). Otherwise, all LAPACK routines are visible by default.

! example.f90
! Example program that solves the matrix equation AX = B using LAPACK.
program main
    implicit none (type, external)
    external :: sgesv
    real     :: a(2, 2)  ! Matrix A.
    real     :: b(2)     ! Vector b/x.
    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, i0)', 'Error: ', rc
        stop
    end if

    print '("Solution (x1, x2): ", f0.4, ", ", f0.4)', b
end program main

Assumed, the LAPACK libraries and modules are located in /usr/local/lib/, compile, dynamically link, and run the example with:

$ gfortran13 -L/usr/local/lib -o example example.f90 -llapack -lblas
$ ./example
Solution (x1, x2): 1.0000, 3.0000

To link LAPACK and BLAS statically, run instead:

$ gfortran13 -o example example.f90 /usr/local/lib/liblapack.a /usr/local/lib/libblas.a

You may have to change the paths to the actual locations.

LAPACK95

LAPACK95 is a modern interface to the FORTRAN 77 LAPACK library that brings several improvements over the original 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 wrappers 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 we can also compile LAPACK95 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    = gfortran13 -ffree-form
FC1   = gfortran13 -ffixed-form
OPTS0 =

Then, compile LAPACK95 with make and the desired precision:

$ cd SRC/
$ make single_double_complex_dcomplex

The Makefile outputs the static library lapack95.a and the Fortran module files inside directory lapack95_modules/.

Solving Linear Equations

In comparison to LAPACK, less boilerplate code is required when calling LAPACK95 routines. For example, the la_gesv() routine takes a minimum of two parameters to solves any real or complex linear system AX = B of arbitrary precision.

call la_gesv(a, b[, ipiv][, info])
Argument Type Description
a type(kind=wp), dimension(:, :) On entry, matrix A, an array of type wp [real, complex]. On exit, the factors from the factorisation.
b type(kind=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 example program calls the LAPACK95 routine to solve the system of linear equations

AX = B

The precision of the variables has to be declared in advance. We set our working precision wp to double precision dp provided by the la_precision module. For single precision, select sp instead.

! example95.f90
program main
    ! Example program that solves the matrix equations AX = B using LAPACK95.
    use :: la_precision, only: wp => dp
    use :: f95_lapack,   only: la_gesv
    implicit none
    real(kind=wp) :: a(2, 2)
    real(kind=wp) :: b(2)

    a = reshape([ 2.0_wp, 3.0_wp, 1.0_wp, 1.0_wp ], [ 2, 2 ])
    b = [ 5.0_wp, 6.0_wp ]

    call la_gesv(a, b)
    print '("Solution (x1, x2): ", f0.4, ", ", f0.4)', b
end program main

LAPACK95 is linked dynamically with -llapack95 -llapack -lblas -ltmglib. The include search path has to be pointed to the directory containing the LAPACK95 module files, e. g., /usr/local/lib/lapack95_modules/:

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

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 squares adjustment are expressed as

least square

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

matrices

We then would have to calculate the normal matrix N = ATA to get the solution to the normal equations

least square

In this particular case, the coefficient matrix has to be filled with the following partial differentiations of the given model equations:

coefficient matrix

The solution vector X contains the transformation parameters a, b, and the translations TE, TN. The parameters are used to determine the rotation angle

rotation angle

and the scale factor

scale factor

With the parameters given, we can easily transform additional coordinates from source (x, y) to target (E, N) system by applying the parameters in the model equations:

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 (E, N) system. 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

The LAPACK95 routine la_gels() returns the minimum-norm least squares solution to one or more real or complex linear systems of the form AX = B, ATX = B, or AHX = B, by using either QR or LQ factorisation of matrix A.

call la_gels(a, b[, trans][, info])
Argument Type Description
a type(kind=wp), dimension(:, :) On entry, matrix A, an rectangular array of shape (:, :) [real, complex]. On Exit, the details of its QR or LQ factorisation.
b type(kind=wp), dimension(:, :) On entry, the matrix B of shape (:) or (:, :). 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 matrix arguments A and B of the LAPACK95 subroutine are passed by reference. In this particular case, with four degrees of freedom given, B will contain the solution and the residual sum of squares of the observations afterwards:

solution

As the input matrices are overwritten, we have to copy them to work arrays first, before calling la_gels(). Once the transformation parameters are determined, the measured points in array m are transformed from source to target system by applying the coordinates and the parameters on the functional model, then storing the result in array t. The complete source code of the example program is listed below.

! adjust.f90
program main
    use :: la_precision, only: wp => dp
    use :: f95_lapack
    implicit none
    integer,       parameter :: N    = 6
    integer,       parameter :: NRHS = 4
    integer,       parameter :: NM   = 8
    real(kind=wp), parameter :: PI   = acos(-1.0_wp)

    real(kind=wp), allocatable :: a(:, :)
    real(kind=wp), allocatable :: c(:), l(:), v(:), x(:), m(:), t(:)
    real(kind=wp)              :: w, r, s
    integer                    :: i

    ! Allocate arrays.
    allocate (a(N, NRHS), c(N), l(N), v(N), x(NRHS), m(NM), t(NM))

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

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

    ! Measured points (x, y) in source 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 ]

    ! Fill coefficient matrix A.
    do i = 1, N - 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

    ! Get transformation parameters `x`, `w`, and `s`.
    call solve(a, l, x, w, s)

    ! Calculate residuals of observations.
    v = matmul(a, x) - l

    ! Calculate adjustment's reference variance, given the
    ! 6 observations and 4 degrees of freedom.
    r = dot_product(v, v) / (size(a, 1) - size(a, 2))

    ! Transform measured points from source (x, y) to target (E, N) system.
    do i = 1, NM, 2
        t(i)     = (s * cos(w)) * m(i) - (s * sin(w)) * m(i + 1) + x(3)
        t(i + 1) = (s * sin(w)) * m(i) + (s * cos(w)) * m(i + 1) + x(4)
    end do

    call output(a, l, x, m, t, v, r, w, s)
contains
    subroutine solve(a, l, x, w, s)
        !! Solves AX = L, then calculates rotation angle and scale factor.
        real(kind=wp), intent(inout) :: a(:, :)  ! Input:  Coefficient matrix A.
        real(kind=wp), intent(inout) :: l(:)     ! Input:  Observations.
        real(kind=wp), intent(inout) :: x(:)     ! Output: Solution [ a, b, Te, Tn ].
        real(kind=wp), intent(out)   :: w, s     ! Output: Rotation angle and scale.
        real(kind=wp), allocatable   :: wa(:, :) ! Work array of `a`.
        real(kind=wp), allocatable   :: wl(:)    ! Work array of `l`.

        allocate (wa(size(a, 1), size(a, 2)), wl(size(l)))
        wa = a
        wl = l

        ! Compute minimum-norm least squares solution to AX = L using LAPACK95.
        ! Work array `wl` contains both the solution and the residuals afterwards.
        ! The solution is stored in `wl(1:size(a, 2))`, and the sum of squares
        ! of the residuals in `wl(size(a, 2) + 1:size(a, 1))`.
        call la_gels(wa, wl)

        x = wl(1:size(a, 2))  ! Solution [ a, b, Te, Tn ].
        w = atan2(x(1), x(2)) ! Calculate rotation angle.
        w = modulo(w, 2 * PI) ! Add 360 deg if angle is negative.
        s = x(1) / cos(w)     ! Scale factor.
    end subroutine solve

    subroutine output(a, l, x, m, t, v, r, w, s)
        real(kind=wp), intent(inout) :: a(:, :)
        real(kind=wp), intent(inout) :: l(:), x(:), m(:), t(:), v(:)
        real(kind=wp), intent(in)    :: r, w, s
        integer                      :: i

        print '(a)', repeat('-', 64)
        print '("-- Two-Dimensional Conformal Coordinate Transformation")'
        print '(a)', repeat('-', 64)

        ! Output L matrix.
        print '(/, "L Matrix")'
        print '(a)', repeat('-', 64)
        print '(f12.3)', [ (l(i), i = 1, size(l)) ]

        ! Output A matrix.
        print '(/, "A Matrix")'
        print '(a)', repeat('-', 64)
        print '(4(f12.3, x))', [ (a(i, :), i = 1, size(a, 1)) ]

        ! Output residuals.
        print '(/, "Transformed Control Points")'
        print '(6x, "E", 12x, "N", 12x, "VE", 11x, "VN")'
        print '(a)', repeat('-', 64)
        print '(4(f12.3, x))', [ (l(i), l(i + 1), v(i), v(i + 1), i = 1, size(l) - 1, 2) ]

        ! Output solution.
        print '(/, "Transformation Parameters")'
        print '(" a = ", f12.3)', x(1)
        print '(" b = ", f12.3)', x(2)
        print '("TE = ", f12.3)', x(3)
        print '("TN = ", f12.3)', x(4)

        print '(/, "Rotation", 11x, "= ", f8.4, " deg")', w * (180 / PI)
        print '("Scale", 14x, "= ", f8.4)', s
        print '("Reference Variance = ", f8.4)', r

        ! Output transformed coordinates.
        print '(/, "Transformed Points")'
        print '(6x, "x", 12x, "y", 12x, "E", 12x, "N")'
        print '(a)', repeat('-', 64)

        print '(4(f12.3, x))', [ ( m(i), m(i + 1), t(i), t(i + 1), i = 1, size(m) - 1, 2) ]
    end subroutine output
end program main

Compile and dynamically link the LAPACK95 example:

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

Then, run the adjustment computation:

$ ./adjust
----------------------------------------------------------------
-- Two-Dimensional Conformal Coordinate Transformation
----------------------------------------------------------------

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

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

Transformed Control Points
      E            N            VE           VN
----------------------------------------------------------------
 1049422.400    51089.200        0.004       -0.029
 1049413.950    49659.300        0.101       -0.077
 1049244.950    49884.950       -0.105        0.106

Transformation Parameters
 a =       -4.512
 b =       -0.254
TE =  1050003.715
TN =    50542.131

Rotation           = 183.2181 deg
Scale              =   4.5196
Reference Variance =   0.0195

Transformed Points
      x            y            E            N
----------------------------------------------------------------
     174.148     -120.262  1049187.361    51040.629
     513.520     -192.130  1047637.713    51278.829
     754.444      -67.706  1046582.113    50656.241
     972.788      120.994  1045644.713    49749.336

The program first outputs the observations and the coefficient matrix, then the transformation parameters, the residuals of the transformed observations, the variance, and the transformed coordinates of points 1, 2, 3, and 4.

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.2621,049,187.3651,040.63
2513.520 –192.1301,047,637.7151,278.83
3754.444 –67.7061,046,582.1150,656.24
4972.788 120.9941,045,644.7149,749.34

Fortran Libraries

References

< BLAS [Index] ANSI Escape Sequences >