使用 MKL 的矢量稀疏矩阵乘法

Question

我有简单的 Fortran 代码，可将 A=[1,1;1,1] 矩阵转换为 CSR 稀疏格式，然后将其乘以 x=(100,200) 作为向量，y=A*x。不幸的是，结果是奇怪的 y=(200,200) 而它应该是 y =(300,300) 向量。谢谢

program main
implicit none
include 'mkl_spblas.fi'
integer :: nzmax, nnz, job( 8 ), m, n, lda, info, irow, k
double precision :: A(2,2)
double precision, allocatable :: Asparse(:)
integer, allocatable :: ia(:), ja(:)
double precision:: x(2)
double precision:: y(2)

A(1,1) = 1.d0
A(1,2) = 1.d0
A(2,1) = 1.d0
A(2,2) = 1.d0
x(1) = 100.d0
x(2) = 200.d0
!! Give an estimate of the number of non-zeros.
nzmax = 4
print *, "nzmax = ", nzmax

m   = size( A, 1 )   !! number of rows
n   = size( A, 2 )   !! number of columns
lda = m              !! leading dimension of A

allocate( Asparse( nzmax ) )
allocate( ja( nzmax ) )   !! <-> columns(:)
allocate( ia( m + 1 ) )   !! <-> rowIndex(:)

job( 1 )   = 0       !! convert dense to sparse A
job( 2:3 ) = 1       !! use 1-based indices
job( 4 )   = 2       !! use the whole A as input
job( 5 )   = nzmax   !! maximum allowed number of non-zeros
job( 6 )   = 1       !! generate Asparse, ia, and ja as output

call mkl_ddnscsr( job, m, n, A, lda, Asparse, ja, ia, info )
if ( info /= 0 ) then
    print *, "insufficient nzmax (stopped at ", info, "row)"; stop
endif

nnz = ia(m+1)
print *, "number of non-zero elements = ", nnz

do irow = 1, m
    !! This loop runs only for rows having nonzero elements.
    do k = ia( irow ), ia( irow + 1 ) - 1
        print "(2i5, f15.8)", irow, ja( k ), Asparse( k )
    enddo
enddo
call mkl_cspblas_dcsrgemv('n', m, Asparse, ia, ja, x, y)
print*, y
end program

Answer 1

那是因为您正在使用 c-style 索引来调用 MKL 子例程。

mkl_cspblas_dcsrgemv 是 zero-based 索引。对于您的程序，您应该使用 mkl_dcsrgemv，即 one-based 索引。

使用 MKL 的矢量稀疏矩阵乘法

Vector- Sparce Matrix multiplication using MKL

fortran

intel-fortran