Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Efficient creation of tridiagonal matrices

Tags:

r

matrix

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?

like image 988
Karsten W. Avatar asked Mar 10 '15 21:03

Karsten W.


People also ask

Which method is used for tridiagonal system of matrix?

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.

Which method is very efficient method to solve a tridiagonal system of linear equations?

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.

How do you reduce a matrix to tridiagonal form?

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.

What is tridiagonal matrix with example?

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 ]


1 Answers

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.

like image 148
Jthorpe Avatar answered Sep 18 '22 10:09

Jthorpe