Calling sparse matrix product in Accelerate

Question

I often use Accelerate framework on macOS to speedup my calculations with matrices. The main program is written in Fortran, and the call to dgemm is essentially identical to that of the documentation on netlib.

In my current project, great performance improvement is expected from the use of sparse-dense matrix multiplication. I found that Accelerate contains a sparse_matrix_product_dense_double function (documentation). According to the documentation, the function should work on macOS 10.11+ (I have 13.2). Here is my question: is it possible to directly call this function from fortran? My naive tests with gfortran

program dgemm_demo
  integer, parameter :: m = 3, n = 4, k = 5 

  double precision, allocatable :: A(:,:), B(:,:), C(:,:)

  allocate(A(m,k), B(k, n), C(m, n))

  A = 1.0D0
  B = 2.0D0

  call dgemm ('N', 'N', M, N, K, 1.0D0, A, m, B, k, 0.0D0, C, m)

  call sparse_matrix_product_dense_double('CblasColMajor', 'CblasTrans', n, 1.0D0,&
    sparse_matrix_create_double(m, k), B, n, C, n)
  
end program dgemm_demo

return linking error

Undefined symbols for architecture arm64

Answer 1

It is possible. Accelerate exposes an API very similar to Sparse BLAS as part of vecLib. (The Sparse BLAS API is defined by the BLAS Technical Forum.)

Search for BLAS.h using Finder to get a full overview of the available routines. You can also locate the header using the command: $ find /Library/Developer -name BLAS.h.

As a vendor-agnostic, yet performant, cross-platform Sparse BLAS alternative, you might also consider librsb.

Here's an example of multiplying a sparse matrix A with a dense matrix B using the Sparse BLAS in Accelerate:

! accelerate.f90
!
!   Example of calling Sparse BLAS from Accelerate framework on MacOS
! 
!   Compile using:
! 
!      gfortran -Wall -O2 accelerate.f90 -framework Accelerate
!
program accelerate

use, intrinsic :: iso_c_binding
implicit none

integer, parameter :: dp = c_double

integer, parameter :: spdim = c_int64_t   ! Sparse Dimension (uint64_t)
integer, parameter :: spstr = c_int64_t   ! Sparse Stride
integer, parameter :: spind = c_int64_t   ! Sparse Index

integer, parameter :: cblas_order = c_int
integer, parameter :: cblas_transpose = c_int

! Transpose flags
integer(cblas_transpose), parameter :: NOTRANS   = 111
integer(cblas_transpose), parameter :: TRANS     = 112
integer(cblas_transpose), parameter :: CONJTRANS = 113

! Order flags
integer(cblas_order), parameter :: ROWMAJOR = 101
integer(cblas_order), parameter :: COLMAJOR = 102

! Sparse Status
integer(c_int), parameter :: SPARSE_SUCCESS             =     0
integer(c_int), parameter :: SPARSE_ILLEGAL_PARAMETER   = -1000
integer(c_int), parameter :: SPARSE_CANNOT_SET_PROPERTY = -1001
integer(c_int), parameter :: SPARSE_SYSTEM_ERROR        = -1002

interface

    !> C = alpha * op(A) * B + C, where A is sparse, B and C are dense
    function dspmm(order,transa,n,alpha,A,B,ldb,C,ldc) &
            bind(c,name="sparse_matrix_product_dense_double")
        import c_int, c_int64_t, c_double, c_ptr
        implicit none
        integer(c_int), value :: order
        integer(c_int), value :: transa
        integer(c_int64_t), value :: n
        real(c_double), value :: alpha
        type(c_ptr), value :: A
        real(c_double) :: B(ldb,*)
        integer(c_int64_t), value :: ldb
        real(c_double) :: C(ldc,*)
        integer(c_int64_t), value :: ldc
        integer(c_int) :: dspmm
    end function

    ! Create sparse matrix of doubles
    function dspcreate(m,n) bind(c,name="sparse_matrix_create_double")
        import c_int64_t, c_ptr
        integer(c_int64_t), value :: m, n
        type(c_ptr) :: dspcreate
    end function

    ! Finalize entry phase
    function spcommit(A) bind(c,name="sparse_commit")
        import c_ptr, c_int
        type(c_ptr), value :: A
        integer(c_int) :: spcommit
    end function

    ! Destroy sparse matrix
    function spdestroy(A) bind(c,name="sparse_matrix_destroy")
        import c_ptr, c_int
        type(c_ptr), value :: A
        integer(c_int) :: spdestroy
    end function

    ! Insert into array in triplet format
    function dspinsert_entries(A,n,val,indx,jndx) &
            bind(c,name="sparse_insert_entries_double")
        import c_ptr, c_int64_t, c_double, c_int
        type(c_ptr), value :: A
        integer(c_int64_t), value :: n
        real(c_double), intent(in) :: val(n)
        integer(c_int64_t), intent(in) :: indx(n), jndx(n)
        integer(c_int) :: dspinsert_entries
    end function

end interface

real(dp) :: A(3,3), B(3,3)
integer :: i

real(dp) :: VA(5), C(3,3)
integer(spdim) :: IA(5), JA(5), nnzA
type(c_ptr) :: spA

!
! Dense matrix multiplication
!

A = reshape([1,0,1,1,1,0,0,0,1],[3,3])
B = reshape([(i,i=1,9)],[3,3])

write(*,'(A,/,3(F10.3,2X))') "Dense result = ", matmul(A,B)

!
! Sparse matrix multiplication
!

nnzA = 5
VA = [ 1, 1,    &
          1,    &
       1,    1]

IA = [1,1,2,3,3] - 1    ! Zero-based indexing expected!
JA = [1,2,2,1,3] - 1

! Initialize matrix of size nrows x ncols
spA = dspcreate(3_spdim,3_spdim)
if (.not. c_associated(spA)) error stop "Sparse matrix creation failed."

! Insert entries from a list of triplets
call spcheck( dspinsert_entries(spA,nnzA,VA,IA,JA) )

! Finalize entry phase by committing
call spcheck( spcommit(spA) )

C = 0.0_dp

! Multiply sparse matrix A with dense matrix B
! Result is stored in C
!
call spcheck( dspmm(COLMAJOR,NOTRANS,3_spdim,1.0_dp,spA,B,3_spdim,C,3_spdim) )
call spcheck( spdestroy(spA) )

write(*,'(/,A,/,3(F10.3,2X))') "Sparse result = ", C

contains

    subroutine spcheck(ierr)
        integer(c_int), intent(in) :: ierr
        
        if (ierr == SPARSE_SUCCESS) &
            return

        write(*,'(A,I0)') 'Error encountered: IERR = ', ierr
        select case(ierr)
        case(SPARSE_ILLEGAL_PARAMETER)
            write(*,'(A)') '  Operation was not completed because one or more'
            write(*,'(A)') '  of the arguments had an illegal value.'
        case(SPARSE_CANNOT_SET_PROPERTY)
            write(*,'(A)') '  Matrix properties can only be set before any'
            write(*,'(A)') '  values are inserted into the matrix.'
        case(SPARSE_SYSTEM_ERROR)
            write(*,'(A)') '  An internal error has occured, such as'
            write(*,'(A)') '  not enough memory.'
        case default
            write(*,'(A)') '  Unknown error code.'
            error stop 1
        end select

    end subroutine

end program