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