## 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:

• linear equations,
• least squares solutions of linear systems of equations,
• eigenvalue problems, and
• singular value problems.

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:

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:

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

``\$ gfortran11 -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    = gfortran11 -ffree-form
FC1   = gfortran11 -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

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) = reshape([ 2.0_wp, 3.0_wp, 1.0_wp, 1.0_wp ], [ 2, 2 ])
real(kind=wp) :: b(2)    = [ 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/`:

``````\$ gfortran11 -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

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:

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

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

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

and the 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 on the model equations:

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:

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 = [ real(kind=wp) :: 121.622, -128.066, &
141.228,  187.718, &
175.802,  135.728 ]

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

! Measured points (x, y) in source system.
m = [ real(kind=wp) :: 174.148, -120.262, &
513.520, -192.130, &
754.444,  -67.706, &
972.788,  120.994 ]

! 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:

``````\$ gfortran11 -I/usr/local/lib/lapack95_modules/ -L/usr/local/lib/ \

``````\$ ./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