Factoring block tridiagonal symmetric positive definite matrices
Goal
Perform Cholesky factorization of a symmetric positive definite block tridiagonal matrix.
Solution
To perform Cholesky factorization of a symmetric positive definite block tridiagonal matrix, with N square blocks of size NB by NB:
Perform Cholesky factorization of the first diagonal block.
Repeat N - 1 times moving down along the diagonal:
Compute the off-diagonal block of the triangular factor.
Update the diagonal block with newly computed off-diagonal block.
Perform Cholesky factorization of a diagonal block.
Source code: see the BlockTDS_SPD/source/dpbltrf.f file in the samples archive available at https://software.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 DPOTRF('L', NB, D, LDD, INFO)
…
DO K=1,N-1
CALL DTRSM('R', 'L', 'T', 'N', NB, NB, 1D0, D(1,(K-1)*NB+1), LDD, B(1,(K-1)*NB+1), LDB)
CALL DSYRK('L', 'N', NB, NB, -1D0, B(1,(K-1)*NB+1), LDB, 1D0, D(1,K*NB+1), LDD)
CALL DPOTRF('L', NB, D(1,K*NB+1), LDD, INFO)
…
END DO
Routines Used
Task |
Routine |
Description |
---|---|---|
Perform Cholesky factorization of diagonal blocks |
DPOTRF |
Computes the Cholesky factorization of a symmetric (Hermitian) positive-definite matrix. |
Compute off-diagonal blocks of the triangular factor |
DTRSM |
Solves a triangular matrix equation. |
Update the diagonal blocks |
DSYRK |
Performs a symmetric rank-k update. |
Discussion
A symmetric positive definite block tridiagonal matrix, with N diagonal blocks Di and N - 1 sub-diagonal blocks Bi of size NB by NB is factored as:
Multiplying the blocks of the matrices on the right gives:
Equating the elements of the original block tridiagonal matrix to the elements of the multiplied factors yields:
Solving for Ci and LiLiT:
Note that the right-hand sides of the equations for LiLiT is a Cholesky factorization. Therefore a routine chol() for performing Cholesky factorization can be applied to this problem using code such as:
L1=chol(D1) do i=1,N-1 Ci=Bi∙Li-T //trsm() Di + 1:=Di + 1 - Ci∙CiT //syrk() Li + 1=chol(Di + 1) end do