2
$\begingroup$

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?

$\endgroup$
1
  • 1
    $\begingroup$ 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. $\endgroup$ Jan 18 at 20:13

1 Answer 1

4
$\begingroup$

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).

$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$
    – Ian Bush
    Jan 18 at 7:45
  • $\begingroup$ 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. $\endgroup$
    – IPribec
    Jan 18 at 9:02
  • $\begingroup$ Are they? I would just calculate the index in place anyway - I just think preprocessing is evil. $\endgroup$
    – Ian Bush
    Jan 18 at 9:39
  • $\begingroup$ An internal procedure can also be used, fully-standard conformant, and gets inlined with -O2. $\endgroup$
    – IPribec
    Jan 18 at 10:32
  • 1
    $\begingroup$ 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. $\endgroup$
    – Ian Bush
    Jan 18 at 11:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.