How can I create a quadratic band matrix, where I give the diagonal and the first diagonal below and above the diagonal? I am looking for a function like
tridiag(upper, lower, main)
where length(upper)==length(lower)==length(main)-1
and returns, for example,
tridiag(1:3, 2:4, 3:6)
[,1] [,2] [,3] [,4]
[1,] 3 1 0 0
[2,] 2 4 2 0
[3,] 0 3 5 3
[4,] 0 0 4 6
Is there an efficient way to do it?
The Thomas algorithm is an efficient way of solving tridiagonal matrix systems. It is based on LU decompo- sition in which the matrix system Mx = r is rewritten as LUx = r where L is a lower triangular matrix and U is an upper triangular matrix.
A tridiagonal system of equations must be solved for every time step, the right-hand side of which is determined by the initial and boundary conditions. The Thomas algorithm developed in Section 9.4 is used to efficiently solve each system.
Reducing a matrix to a tridiagonal form is an iterative process which uses Jacobi rotations to reduce matrix el- ements to zero. The purpose of this research project is to implement an existing algo- rithm for tridiagonal reduction using CUDA, thus leveraging the parallelism present in GPUs to accelerate the process.
A tridiagonal matrix is a matrix that has non-zero elements only at the main diagonal, diagonal below and above it. All other elements are zero. For this reason tridiagonal matrices of dimension smaller than or equal to 3 seem meaningless. Example 1: [a11, a22, 0 , 0 , 0 , 0 ]
This function will do what you want:
tridiag <- function(upper, lower, main){
out <- matrix(0,length(main),length(main))
diag(out) <- main
indx <- seq.int(length(upper))
out[cbind(indx+1,indx)] <- lower
out[cbind(indx,indx+1)] <- upper
return(out)
}
Note that when the index to a matrix is a 2 column matrix, each row in that index is interpreted as the row and column index for a single value in the vector being assigned.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With