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

• 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``

### 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 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: 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 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: We may apply the normal matrix N = ATA to get the solution to the normal equations: But LAPACK95 lets use solve the normal equations directly. The required transformation parameters include the translation factors TX and TY, the rotation angle and the scale factor We transform additional coordinates given in the source system into the target system by applying the parameters: 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/ \

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``````
• E. Anderson, Z. Bai, C. Bischof, et. al: LAPACK Users’ Guide. Third Edition. 1999
• V. A. Barker, L. S. Blackford, J. J. Dongarra, et. al: LAPACK95 Users’ Guide. 2001
• Ch. D. Ghilani: Adjustment Computations: Spatial Data Analysis. Fifth Edition. 2010