使用Intel OneAPI MKL的mkl_dcsrgemv处理大维度稀疏矩阵向量乘积异常问题

问题背景 #

使用Fortran搭配Intel OneAPI编译器ifx及MKL库计算大维度稀疏矩阵与向量的乘积时，用integer(4)描述稀疏矩阵索引数组可正常运行，但将索引数组改为integer(8)后，运行结果异常。

正常运行的代码（integer(4)索引） #

program main 
  implicit none
  integer :: n, ja(9), ia(7)
  real(8) :: a(9), x(6), y(6)
  integer :: i, j
  integer(4), allocatable :: ic(:), jc(:)
  real(8), allocatable :: c(:), xc(:), yc(:)
  integer(8) :: nc

  a  = [1.d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 4.0d0, 3.0d0, 2.0d0, 1.0d0]
  ja = [1,2,3,4,6,5,4,5,6]
  ia = [1,2,3,4,6,7,10]

  x  = [1,2,3,4,5,6]
  call mkl_dcsrgemv('n', size(x), a, ia, ja, x, y)
  print*, y

  if (.true.) then 
  nc = 6     
  allocate(c(9))
  allocate(jc(9)) 
  allocate(ic(nc+1)) 
  allocate(xc(6))  
  allocate(yc(6))  
  c  = [1.d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 4.0d0, 3.0d0, 2.0d0, 1.0d0]
  jc = [1,2,3,4,6,5,4,5,6]
  ic = [1,2,3,4,6,7,10]

  xc = [1,2,3,4,5,6]
  yc = 0.0d0
  call mkl_dcsrgemv('n', nc, c, ic, jc, xc, yc)
  
  print *, "x:", xc
  print *, "y:", yc
  end if 

end program main 

正常运行输出 #

proponent@node3:~/testfortransparse$ make clean
rm -f main.o common_data.mod main
proponent@node3:~/testfortransparse$ make 
ifx -qopenmp -O3 -qmkl=parallel  -c main.f90 
ifx -qopenmp -O3 -qmkl=parallel  -o main main.o 
proponent@node3:~/testfortransparse$ make run 
./main
   1.00000000000000        4.00000000000000        9.00000000000000     
   46.0000000000000        20.0000000000000        28.0000000000000     
 x:   1.00000000000000        2.00000000000000        3.00000000000000     
   4.00000000000000        5.00000000000000        6.00000000000000     
 y:   1.00000000000000        4.00000000000000        9.00000000000000     
   46.0000000000000        20.0000000000000        28.0000000000000     

修改为integer(8)索引后的异常输出 #

修改代码中integer(4), allocatable :: ic(:), jc(:)为integer(8), allocatable :: ic(:), jc(:)后，运行结果异常：

proponent@node3:~/testfortransparse$ make 
ifx -qopenmp -O3 -qmkl=parallel  -c main.f90 
ifx -qopenmp -O3 -qmkl=parallel  -o main main.o 
proponent@node3:~/testfortransparse$ make run 
./main
   1.00000000000000        4.00000000000000        9.00000000000000     
   46.0000000000000        20.0000000000000        28.0000000000000     
 x:   1.00000000000000        2.00000000000000        3.00000000000000     
   4.00000000000000        5.00000000000000        6.00000000000000     
 y:  0.000000000000000E+000   1.00000000000000       0.000000000000000E+000
   1.00000000000000       0.000000000000000E+000   7.00000000000000     

原因分析 #

MKL的mkl_dcsrgemv函数是为32位整数（integer(4)）索引设计的，当传入64位整数（integer(8)）类型的索引数组时，函数会错误地解析内存布局，导致计算逻辑混乱，最终输出异常结果。MKL针对64位整数索引提供了独立的函数接口，不能直接用64位变量调用32位版本的函数。

解决方案 #

使用MKL专门为64位整数索引设计的mkl_dcsrgemv_64函数，同时确保所有相关参数的类型匹配：

修改后的代码 #

program main 
  implicit none
  integer :: n, ja(9), ia(7)
  real(8) :: a(9), x(6), y(6)
  integer :: i, j
  integer(8), allocatable :: ic(:), jc(:)  ! 64位索引数组
  real(8), allocatable :: c(:), xc(:), yc(:)
  integer(8) :: nc  ! 矩阵维度使用64位整数

  a  = [1.d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 4.0d0, 3.0d0, 2.0d0, 1.0d0]
  ja = [1,2,3,4,6,5,4,5,6]
  ia = [1,2,3,4,6,7,10]

  x  = [1,2,3,4,5,6]
  call mkl_dcsrgemv('n', size(x), a, ia, ja, x, y)
  print*, y

  if (.true.) then 
  nc = 6     
  allocate(c(9))
  allocate(jc(9)) 
  allocate(ic(nc+1)) 
  allocate(xc(6))  
  allocate(yc(6))  
  c  = [1.d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 4.0d0, 3.0d0, 2.0d0, 1.0d0]
  jc = [1,2,3,4,6,5,4,5,6]
  ic = [1,2,3,4,6,7,10]

  xc = [1,2,3,4,5,6]
  yc = 0.0d0
  ! 调用64位版本的MKL函数
  call mkl_dcsrgemv_64('n', nc, c, ic, jc, xc, yc)
  
  print *, "x:", xc
  print *, "y:", yc
  end if 

end program main 

修改后的运行结果 #

proponent@node3:~/testfortransparse$ make clean && make && make run
rm -f main.o common_data.mod main
ifx -qopenmp -O3 -qmkl=parallel  -c main.f90 
ifx -qopenmp -O3 -qmkl=parallel  -o main main.o 
./main
   1.00000000000000        4.00000000000000        9.00000000000000     
   46.0000000000000        20.0000000000000        28.0000000000000     
 x:   1.00000000000000        2.00000000000000        3.00000000000000     
   4.00000000000000        5.00000000000000        6.00000000000000     
 y:   1.00000000000000        4.00000000000000        9.00000000000000     
   46.0000000000000        20.0000000000000        28.0000000000000     

关键注意点 #

• 所有涉及稀疏矩阵索引的数组（行指针、列索引）必须统一使用integer(8)类型
• 调用MKL的64位版本函数（后缀为_64）
• 矩阵维度参数（如示例中的nc）也建议使用integer(8)，避免类型不匹配

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