In case anyone finds it helpful, I’ve put together a small example using the general banded solvers from various Fortran libraries.

Origin | Factorization | Solve | Driver | Pivoting | Comment |
---|---|---|---|---|---|

LAPACK | DGBTRF | DGBTRS | DGBSV | Y | diagonals stored as rows |

LINPACK | DGBFA, DGBCO | DGBSL | / | Y | diagonals stored as rows |

SLATEC | DNBFA, DNBCO | DNBSL | DNBFS | Y | diagonals stored as columns |

PPPACK | BANFAC | BANSLV | / | N | diagonals stored as rows; only suitable for diagonally-dominant matrices |

BODE [1] | LTRI | TSOL | / | Y | diagonals stored as columns; only for systems with `ml == mu` ; based on the procedures `bandet1` and `bansol1` [3] (in Algol) |

Thorson [2] | band | solve | / | Y | diagonals stored as columns; only for systems with `ml == mu` |

Engeln-Müllges & Uhlig [4] | BANDP | BANDS | BAND | Y | diagonals stored as columns |

Numerical Recipes [5] | bandec | banbks | / | Y | based on the procedures `bandet1` and `bansol1` [3] |

Are there any other banded solvers available? If yes, please share them here.

Literature:

- [1] Hopkins, T. R., & Wait, R. (1978). A comparison of Galerkin, collocation and the method of lines for partial differential equations. International Journal for Numerical Methods in Engineering, 12(7), 1081-1107. https://doi.org/10.1002/nme.1620120703
- [2] Jeff Thorson, Gaussian Elimination on a Banded Matrix
- [3] Martin, R. S., & Wilkinson, J. H. (1967). Solution of symmetric and unsymmetric band equations and the calculation of eigenvectors of band matrices.
*Numerische Mathematik*,*9*, 279-301. Solution of symmetric and unsymmetric band equations and the calculation of eigenvectors of band matrices | Numerische Mathematik - [4] Engeln-Müllges, G., & Uhlig, F. (1996).
*Numerical Algorithms with Fortran*. Springer, Berlin. Numerical Algorithms with Fortran | SpringerLink - [5] Press, W. H., et al. (1992). Numerical Recipes in Fortran 77: the art of scientific computing, Cambridge University Press, http://s3.amazonaws.com/nrbook.com/book_F210.html

Demo:

```
program example
use, intrinsic :: iso_fortran_env, only: error_unit
implicit none
integer, parameter :: wp = kind(1.0d0)
! Constants
integer, parameter :: n = 5
integer, parameter :: ml = 1, mu = 2
integer, parameter :: bw = ml + mu + 1
! Dense matrix and vectors
real(wp) :: ad(n,n), b(n), xla(n), xsl(n)
real(wp) :: wrk(n)
integer :: ipiv(n)
! SLATEC format (extra column needed for pivoting)
real(wp) :: as(n,bw+ml)
integer :: itask, ind
! LAPACK format (additional ml rows needed by factorization)
real(wp) :: ab(bw+ml,5)
integer :: info
! Procedures
external :: dgbsv, dnbfs, banfac, banslv
! Other variables
integer :: m1, m21, m3
! Dense matrix
ad(1,:) = [1, 2, 3, 0, 0]
ad(2,:) = [4, 5, 6, 7, 0]
ad(3,:) = [0, 8, 9, 1, 2]
ad(4,:) = [0, 0, 3, 4, 5]
ad(5,:) = [0, 0, 0, 6, 7]
! Right-hand side
! A * x = [6, 22, 20, 12, 13]^T, x = [1]_n
b = sum(ad,dim=2)
! Convert dense to banded (extra storage is provided for the factorization)
ab = dense2bandedLA(ad,ml,mu)
as = dense2bandedSLATEC(ad,ml,mu)
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
!
! Solve banded system A x = b with LAPACK
!
xla = b
call dgbsv(n,ml,mu,1,ab,n,ipiv,xla,n,info)
if (info /= 0) then
write(error_unit,*) "[dgbsv] an error occured"
error stop info
end if
write(*,'(A,*(/,G0))') "x (LAPACK) = ", xla
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
!
! Solve banded system A x = b with SLATEC
!
xsl = b
itask = 1
call dnbfs(as,n,n,ml,mu,xsl,itask,ind,wrk,ipiv)
if (ind < 0) then
write(error_unit,*) "[dnbfs] an error occured"
error stop ind
end if
write(*,'(A,*(/,G0))') "x (SLATEC) = ", xsl
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
!
! Solve banded system A x = b with PPPACK
! (not intended for arbitrary matrices as it doesn't use pivoting,
! also watch out for single/double precision!)
!
ab = dense2bandedLA(ad,ml,mu)
b = sum(ad,dim=2)
xla = b
! Shift into first row (PPPACK doesn't need extra storage)
ab(1:bw,:) = ab(ml+1:ml+bw,:)
call banfac(ab,size(ab,1),n,ml,mu,info)
if (info /= 1) then
write(error_unit,*) '[banfac] Factorization failed.'
error stop info
end if
call banslv(ab,size(ab,1),n,ml,mu,xla)
write(*,'(A,*(/,G0))') "x (BANSLV) = ", xla
contains
!
! Convert to banded format using LAPACK storage
!
! Description found here:
! - https://www.netlib.org/lapack/lug/node124.html
function dense2bandedLA(A, ml, mu) result(AB)
integer, intent(in) :: ml, mu
real(wp), intent(in) :: A(:,:)
real(wp) :: AB(2*ml+mu+1, size(A,1))
integer :: j, i
ab = 0
do j = 1, size(A,1)
do i = max(1, j - mu), min(size(A,1), j+ml)
AB(ml+mu+1+i-j,j) = A(i,j)
end do
end do
end function
!
! Convert to banded using SLATEC storage
! (diagonals are stored as columns)
!
function dense2bandedSLATEC(A, ml, mu) result(ab)
integer, intent(in) :: ml, mu
real(wp), intent(in) :: A(:,:)
real(wp) :: ab(size(A,1),2*ml+mu+1)
integer :: i, j1, j2, j, k
ab = 0
do i = 1, n
j1 = max(1,i-ml)
j2 = min(n,i+mu)
do j = j1, j2
k = j - i + ml + 1
ab(i,k) = a(i,j)
end do
end do
end function
end program
```