能否在Fortran中直接调用Accelerate的sparse_matrix_product_dense_double函数？

我在macOS系统上常用Accelerate框架加速矩阵运算，主程序用Fortran编写，其中dgemm的调用方式和netlib文档完全一致。

当前项目中，使用稀疏-稠密矩阵乘法有望大幅提升性能。我发现Accelerate框架包含sparse_matrix_product_dense_double函数，该函数支持macOS 10.11及以上版本（我的系统是13.2）。

我的问题是：能否直接从Fortran中调用该函数？我用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

测试时出现链接错误：

arm64架构下存在未定义符号

解决思路 #

1. 函数命名与绑定问题 #

Accelerate的稀疏矩阵函数是C语言接口，Fortran直接调用会遇到名字改编问题：C函数名在链接时原样保留，但Fortran默认会将函数名转为小写并添加下划线（比如dgemm会变成dgemm_），两者命名规则不匹配。

必须用Fortran的bind(C)特性显式声明函数接口，告诉编译器这是C风格函数：

interface
  subroutine sparse_matrix_product_dense_double(order, transa, n, alpha, A, B, ldb, C, ldc) bind(C, name='sparse_matrix_product_dense_double')
    use, intrinsic :: iso_c_binding
    integer(c_int), value :: order, transa, n, ldb, ldc
    real(c_double), value :: alpha
    type(c_ptr), value :: A  ! 稀疏矩阵以C指针类型传递
    real(c_double), intent(in) :: B(ldb, *)
    real(c_double), intent(inout) :: C(ldc, *)
  end subroutine
end interface

2. 稀疏矩阵创建逻辑错误 #

你代码中直接调用sparse_matrix_create_double(m, k)是错误的：该函数返回C语言的sparse_matrix_t类型，Fortran需用c_ptr接收；且创建稀疏矩阵必须传入非零元的具体数据（行索引、列指针、值），不能仅传入行列数。正确步骤：

• 按CSC/CSR格式准备稀疏矩阵的非零元数据
• 调用sparse_matrix_create_double创建矩阵对象（用c_ptr存储）
• 使用完毕后调用sparse_matrix_destroy释放资源

3. 编译链接参数 #

编译时必须显式链接Accelerate框架，gfortran的编译命令需添加-framework Accelerate：

gfortran your_code.f90 -o your_program -framework Accelerate

完整修正示例片段 #

program sparse_dense_demo
  use, intrinsic :: iso_c_binding
  implicit none

  ! 声明Accelerate稀疏矩阵相关函数的C接口
  interface
    ! 创建双精度稀疏矩阵（CSC格式）
    function sparse_matrix_create_double(m, n, col_ptrs, row_indices, values, flags) bind(C, name='sparse_matrix_create_double')
      import :: c_int, c_double, c_ptr
      integer(c_int), value :: m, n, flags
      integer(c_int), intent(in) :: col_ptrs(*), row_indices(*)
      real(c_double), intent(in) :: values(*)
      type(c_ptr) :: sparse_matrix_create_double
    end function

    ! 稀疏-稠密矩阵乘法
    subroutine sparse_matrix_product_dense_double(order, transa, n, alpha, A, B, ldb, C, ldc) bind(C, name='sparse_matrix_product_dense_double')
      import :: c_int, c_double, c_ptr
      integer(c_int), value :: order, transa, n, ldb, ldc
      real(c_double), value :: alpha
      type(c_ptr), value :: A
      real(c_double), intent(in) :: B(ldb, *)
      real(c_double), intent(inout) :: C(ldc, *)
    end subroutine

    ! 销毁稀疏矩阵
    subroutine sparse_matrix_destroy(A) bind(C, name='sparse_matrix_destroy')
      import :: c_ptr
      type(c_ptr), value :: A
    end subroutine
  end interface

  integer, parameter :: m = 3, n = 4, k = 5, nnz = 3  ! 假设稀疏矩阵含3个非零元
  integer(c_int) :: col_ptrs(k+1), row_indices(nnz)
  real(c_double) :: values(nnz), B(k, n), C(m, n)
  type(c_ptr) :: sparse_A
  integer(c_int), parameter :: CblasColMajor = 102, CblasNoTrans = 111

  ! 初始化稠密矩阵与结果矩阵
  B = 2.0D0
  C = 0.0D0

  ! 构造CSC格式的稀疏矩阵数据（示例：对角线上的3个1.0）
  col_ptrs = [1,2,3,4,5,6]  ! 列指针（对应Accelerate的CSC格式要求）
  row_indices = [0,1,2]     ! 0-based行索引
  values = [1.0D0, 1.0D0, 1.0D0]

  ! 创建稀疏矩阵对象
  sparse_A = sparse_matrix_create_double(m, k, col_ptrs, row_indices, values, 0)

  ! 执行稀疏-稠密矩阵乘法
  call sparse_matrix_product_dense_double(CblasColMajor, CblasNoTrans, n, 1.0D0, sparse_A, B, k, C, m)

  ! 释放稀疏矩阵资源
  call sparse_matrix_destroy(sparse_A)

  ! 输出结果
  print *, "结果矩阵C:"
  print *, C

end program sparse_dense_demo

关键注意事项 #

• 所有C接口函数必须用bind(C)声明，保证命名匹配
• 稀疏矩阵格式（CSR/CSC）需与函数要求一致，sparse_matrix_product_dense_double默认使用CSC格式
• 类型需严格对应：C的int对应Fortran的c_int，双精度对应c_double，指针用c_ptr
• 编译时必须添加-framework Accelerate链接框架

内容的提问来源于stack exchange，提问作者yarchik