CSR space matrix vector product with stdlib
tof92130 opened this issue · 3 comments
Description
Hello,
1/ With gnu compilers, I am using subroutine spmv from stdlib in order to realize a CSR sparse matrix vector product :
call spmv(matrix=mat%STDL, vec_x=x0, vec_y=x1, alpha=alpha, beta=beta) !> x1 = alpha mat%STDL x0 + beta x0
2/ With Intel compilers, I am using mkl_dcsrmv subroutine for same operation.
3/ I also have a my own procedure for this operation.
If beta=0 we should find x1 = mat%STDL x0
it I the case with mkl_dcsrmv and my own procedure but to obtain this result with procedure spmv from sodlib I must initialize x1(:) to 0d0 before (not good for bandwidth).
Here is my own procedure:
do iRow=1,size(mat%STDL%rowptr)-1
x1(iRow)=beta*x1(iRow)
!$OMP SIMD
do jMat=mat%STDL%rowptr(iRow),mat%STDL%rowptr(iRow+1)-1
x1(iRow)=x1(iRow)+alpha*mat%STDL%data(jMat)*x0(mat%STDL%col(jMat))
enddo
enddo
line 2 initialise x1 to zero if beta=0.
Expected Behaviour
initialize x1 to zero when beta=0
Version of stdlib
Platform and Architecture
Mac/OS ARM and Intel
Additional Information
No response
Hi @tof92130, thanks for signaling this. I'm intrigued though:
At the beginning of the procedure one can find
stdlib/src/stdlib_sparse_spmv.fypp
Lines 148 to 152 in 60d0a76
Which should clean y following:
y = alpha*op(matrix)*x + beta * y
y should indeed be zero if beta == 0. Could you show us the results you are obtaining?
Indeed, it should work...
I investigate my program to see if the problem couldn't come from elsewhere (?)
can you please close this issue.
I solved the problem by initialising x1(:) and x2(:) just after first initialisation. Size of x1(:) and x2(:) are larger than size(mat%STDL%rowptr)-1 because I also include place for MPI exchange and I probably introduced bad values because of that.
intel compilers initialize to zero after allocation, gnu compilers do not.