Solution of a banded system of linear equations
Solves the linear system of equations
Ax = B
by Gaussion elimination
A has to be square, and banded, i.e. with the only nonzero
elements in bands near the diagonal.
A is either inputted as a full square matrix or as the non-zero
uses lapack subroutine dgbsv (FORTRAN)
either a matrix containing the (nonzero) bands, rotated row-wise (anti-clockwise) only, or a full square matrix.
number of nonzero bands above the diagonal; ignored if
number of nonzero bands below the diagonal; ignored if
Right-hand side of the equations, a vector with length = number
of rows of
If the input matrix
abd is square, it is assumed that the full,
square A is inputted, unless
full is set to
abd is not square, then the number of columns denote the
number of unknowns, while the number of rows equals the nonzero bands,
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
A similar function but that requires a totally different input can now
also be found in the
Karline Soetaert <email@example.com>
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)