# What is the correct way of computing $LL^T$ in packed format with blas/lapack

I have a triangular matrix stored in packed format (ie $$L$$). I need to compute $$LL^T$$ (not the decomposition, just the multiplication). What would be the preferred way of computing this with blas/lapack?.

One obvious way would be to follow the following steps:

1. Convert packed format to full triangular (DTPTTR)
2. Compute triangular multiplication in full format (DLAUU2)
3. Convert full symmetric matrix back into packed format (DTRTTP)

However, this approach seems extremely bad to me. In the first step, it would require doubling memory requirement and the third step is needed just to reverse the first step.

Is there any better/preferred way of carrying out such an operation still using blas/lapack?

• Be on the lookout for performance issues. Some BLAS libraries have highly optimized routines for commonly used functions but leave uncommonly used functions (like functions involving packed storage) unoptimized. Jan 18 at 20:13

The closest routine in BLAS for this would be DSYRK which performs the update C := alpha*A*A**T + beta*C, however C and A use full storage.

Here is a solution with a naive triple loop:

module matmult
implicit none
public
contains

! SLLTPM performs the update
!
!    B := A*A**T
!
! where B is a n by n symmetric matrix, and A is an n by n
! lower-triangular matrix. Packed storage is used for both A and B.
subroutine slltpm(n,a,b)
implicit none
integer, intent(in) :: n
real, intent(in) :: a(*)
real, intent(out) :: b(*)

integer :: i, j, k

do j = 1, n
b(id(j,j):id(n,j)) = 0
do k = 1, j
!$omp simd do i = j, n b(id(i,j)) = b(id(i,j)) + a(id(i,k))*a(id(j,k)) end do end do end do contains integer function id(i,j) integer, intent(in) :: i, j id = i + (2*n-j)*(j-1)/2 end function end subroutine end module program main use matmult, only: slltpm implicit none integer, parameter :: n = 4 real :: ad(n,n), bd(n,n) integer, parameter :: np = n*(n+1)/2 real :: a(np), b(np) ! A := L a = [real :: 1, 2, 3, 4, 5, 6, 7, 8, 9, 10] ! B := AA^T call slltpm(n,a,b) call print_packed(n,a,"a (packed)") call print_packed(n,b,"b (packed)") call unpack(n,a,ad) bd = matmul(ad,transpose(ad)) call print_dense(n,ad,"a (dense)") call print_dense(n,bd,"b (dense)") contains subroutine unpack(n,a,b) integer, intent(in) :: n real, intent(in) :: a(*) real, intent(out) :: b(n,n) integer :: id, i, j id(i,j) = i + (2*n-j)*(j-1)/2 b = 0 do j = 1, n do i = j, n b(i,j) = a(id(i,j)) end do end do end subroutine subroutine print_packed(n,a,str) integer, intent(in) :: n real, intent(in) :: a(*) character(*), intent(in) :: str integer :: id, i, j id(i,j) = i + (2*n-j)*(j-1)/2 print *, str do i = 1, n print *, a( [(id(i,j), j = 1, i)] ) end do end subroutine subroutine print_dense(n,a,str) integer, intent(in) :: n real, intent(in) :: a(n,n) character(*), intent(in) :: str integer :: i, j print *, str do i = 1, n print *, (a(i,j), j = 1, n) end do end subroutine end program  Judging by the output in Compiler Explorer, both gfortran and ifx use an AVX fused multiply-add instruction for the innermost loop: gfortran v13.2 .L5: vmovups ymm0, YMMWORD PTR [rsi+rax] vfmadd213ps ymm0, ymm2, YMMWORD PTR [rdx+rax] vmovups YMMWORD PTR [rdx+rax], ymm0 add rax, 32 cmp r9, rax jne .L5  ifx v2024.0 .LBB0_19: vmovups ymm2, ymmword ptr [r8 + 4*r9] vfmadd213ps ymm2, ymm1, ymmword ptr [rdx + 4*r9] vmovups ymmword ptr [rdx + 4*r9], ymm2 add r9, 8 cmp r9, rcx jb .LBB0_19  The next step toward better performance would be cache blocking to optimize memory access, and multi-threaded execution (with !$omp parallel do schedule(...) on the outer loop).

• external :: llt_packed? I suspect you mean external :: lltpm, but better would be an interface block (and best a module). And you could avoid the need for C preprocessing with the marginally better statement function - at least it's standard Fortran. Jan 18 at 7:45
• Thanks for catching that error. Statement functions are marked as obsolescent when you compile with options such as -std=f95 and higher. Standardization of the preprocessor is being discussed for F202Y. Jan 18 at 9:02
• Are they? I would just calculate the index in place anyway - I just think preprocessing is evil. Jan 18 at 9:39
• An internal procedure can also be used, fully-standard conformant, and gets inlined with -O2. Jan 18 at 10:32
• Agreed - but not guaranteed to be inlined, it's implementation dependent. That's why I'd do the calculation in place if I really care about the performance. But this is all small beer, nice answer. Jan 18 at 11:15