Solves the linear system of equations

*Ax = B*

by Gaussion elimination

where `A`

has to be square, and *banded*, i.e. with the only nonzero
elements in bands near the diagonal.

The matrix `A`

is either inputted as a full square matrix or as the non-zero
bands.

uses lapack subroutine dgbsv (FORTRAN)

1 2 | ```
Solve.banded(abd, nup, nlow, B = rep(0, times = ncol(abd)),
full = (nrow(abd) == ncol(abd)))
``` |

`abd ` |
either a matrix containing the (nonzero) bands, rotated row-wise (anti-clockwise) only, or a full square matrix. |

`nup ` |
number of nonzero bands above the diagonal; ignored if |

`nlow ` |
number of nonzero bands below the diagonal; ignored if |

`B ` |
Right-hand side of the equations, a vector with length = number
of rows of |

`full ` |
if |

If the input matrix `abd`

is square, it is assumed that the full,
square A is inputted, unless `full`

is set to `FALSE`

.

If `abd`

is not square, then the number of columns denote the
number of unknowns, while the number of rows equals the nonzero bands,
i.e. nup+nlow+1

matrix with the solution, `X`

, of the banded system of equations A X =B,
the number of columns of this matrix = number of columns of `B`

.

A similar function but that requires a totally different input can now
also be found in the `Matrix`

package

Karline Soetaert <karline.soetaert@nioz.nl>

J.J. Dongarra, J.R. Bunch, C.B. Moler, G.W. Stewart, LINPACK Users' Guide, SIAM, 1979.

`Solve.tridiag`

to solve a tridiagonal system of linear equations.

`Solve`

the generalised inverse solution,

`solve`

the R default

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 | ```
# 1. Generate a banded matrix of random numbers, full format
nup <- 2 # nr nonzero bands above diagonal
ndwn <- 3 # nr nonzero bands below diagonal
nn <- 10 # nr rows and columns of A
A <- matrix(nrow = nn, ncol = nn, data = runif(1 : (nn*nn)))
A [row(A) < col(A) - nup | row(A) > col(A) + ndwn] <- 0
diag(A) <- 1 # 1 on diagonal is easily recognised
# right hand side
B <- runif(nrow(A))
# solve it, using the default solver and banded (inputting full matrix)
Full <- solve(A, B)
Band1 <- Solve.banded(A, nup, ndwn, B)
# 2. create banded form of matrix A
Aext <- rbind(matrix(ncol = ncol(A), nrow = nup, 0),
A,
matrix(ncol = ncol(A), nrow = ndwn, 0))
abd <- matrix(nrow = nup + ndwn + 1, ncol = nn,
data = Aext[col(Aext) <= row(Aext) &
col(Aext) >= row(Aext) - ndwn - nup])
# print both to screen
A
abd
# solve problem with banded version
Band2 <- Solve.banded(abd, nup, ndwn, B)
# compare 3 methods of solution
cbind(Full, Band1, Band2)
# same, now with 3 different right hand sides
B3 <- cbind(B, B*2, B*3)
Solve.banded(abd, nup, ndwn, B3)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.