C语言 Mkl 矩阵乘法,使用MKL的Vector-Sparce矩阵乘法

我有简单的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