Solving a system of linear equations with a block tridiagonal symmetric positive definite coefficient matrix
Goal
Solve a system of linear equations with a Cholesky-factored symmetric positive definite block tridiagonal coefficient matrix.
Solution
Given a coefficient symmetric positive definite block tridiagonal matrix (with square blocks each of the same NB-by-NB size) is LLT factored, the solving stage consists of:
- Solve the system of linear equations with a lower bidiagonal coefficient matrix which is composed of N by N blocks of size NB by NB and with diagonal blocks which are lower triangular matrices:
Solve the N local systems of equations with lower triangular diagonal blocks of size NB by NB which are used as coefficient matrices and respective parts of the right hand side vectors.
Update the local right hand sides.
Solve the system of linear equations with an upper bidiagonal coefficient matrix which is composted of block size N by N blocks of size NB by NB and with diagonal blocks which are upper triangular matrices.
Solve the local systems of equations.
Update the local right hand sides.
Source code: see the BlockTDS_SPD/source/dpbltrs.f file in the samples archive available at https://www.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook-samples-120115.zip.
Cholesky factorization of a symmetric positive definite block tridiagonal matrix
…
…
CALL DTRSM('L', 'L', 'N', 'N', NB, NRHS, 1D0, D, LDD, F, LDF)
DO K = 2, N
CALL DGEMM('N', 'N', NB, NRHS, NB, -1D0, B(1,(K-2)*NB+1), LDB, F((K-2)*NB+1,1), LDF, 1D0, F((K-1)*NB+1,1), LDF)
CALL DTRSM('L','L', 'N', 'N', NB, NRHS, 1D0, D(1,(K-1)*NB+1), LDD, F((K-1)*NB+1,1), LDF)
END DO
CALL DTRSM('L', 'L', 'T', 'N', NB, NRHS, 1D0, D(1,(N-1)*NB+1), LDD, F((N-1)*NB+1,1), LDF)
DO K = N-1, 1, -1
CALL DGEMM('T', 'N', NB, NRHS, NB, -1D0, B(1,(K-1)*NB+1), LDB, F(K*NB+1,1), LDF, 1D0, F((K-1)*NB+1,1), LDF)
CALL DTRSM('L','L', 'T', 'N', NB, NRHS, 1D0, D(1,(K-1)*NB+1), LDD, F((K-1)*NB+1,1), LDF)
END DO
…
Routines Used
Task |
Routine |
Description |
---|---|---|
Solve a local system of linear equations |
DTRSM |
Solves a triangular matrix equation. |
Update the local right hand sides |
DGEMM |
Computes a matrix-matrix product with general matrices. |
Discussion
Consider a system of linear equations described by:
![](/content/dam/docs/us/en/cookbook/2024-0/F8A364CA-94CE-4612-B306-2611BF395A451.png)
Assume that matrix A is a symmetric positive definite block tridiagonal coefficient matrix with all blocks of size NB by NB. A can be factored as described in Factoring block tridiagonal symmetric positive definite matrices to give:
![](/content/dam/docs/us/en/cookbook/2024-0/F8A364CA-94CE-4612-B306-2611BF395A452.png)
Then the algorithm to solve the system of equations is:
Solve the system of linear equations with a lower bidiagonal coefficient matrix in which the diagonal blocks are lower triangular matrices:
Y1=L1-1 F1 //trsm() do i=2,N Gi=Fi - Ci - 1 Yi - 1 //gemm() Yi=Li-1 Gi //trsm() end do
Solve the system of linear equations with an upper bidiagonal coefficient matrix in which the diagonal blocks are upper triangular matrices:
XN=LN-T YN //trsm() do i=N-1,1,-1 Zi=Fi-CiT Xi + 1 //gemm() Xi=Li-T Zi //trsm() end do