Problem with sparse BLAS Level 2 function mkl_scsrgemv()

[Wu__Wentao]

Hi!

When programming with mkl sparse BLAS level 2 fucntion to implement a matrix-vector mulitplication, I met a trouble in function mlk_?csrgemv(). A simplified and clear code is given below, in which the same problem happened.

!-----------------------------------------------

program main

implicit none

real, dimension(12) :: val = (/-1, 1,-1, &

1,-1, 1, &

1,-1, 1, &

-1, 1,-1 /)

integer, dimension(12) :: column = (/1,2,3, 4,5,6, 7,8,9, 1,4,7/)

integer, dimension(5) :: rowIndex= (/1,4,7, 10,13/)

real, dimension(9) :: v = 1.0

real, dimension(9) :: res = 0

character :: transa = 'n'

integer :: row_num = 9

call mkl_scsrgemv(transa, row_num, val, rowIndex, column, v, res)

write(*,*) "result = "

write(*,*) res

stop

end program

!------------------------------------------

Error info is: Unhandled exception at 0x0117d008 in SparseBLAS.exe: 0xC0000005: Access violation reading location 0xff25e07c.

It shoud be pointed out that the matrix is not a square matrix of the dimension m*m.

But another non square sparse matrix,

val = (/-1,1,1,-1,1,-1,-1,1/)

column = (/1,2,3,4,5,6,2,4/)

rowIndex = (/1,3,5,7,9/)

worked successfully.

My question is, what is wrong with the code? Can the nonsquare matrix vector multiplication be implemented with CSR(3 array variation) data structure? And How?

(Platform: Intel Pentium T4200, VS2008, MKL 10.3)

Expecting for your help and Thanks a lot.

[Victor_Gladkikh]

Hi David

As I understand your first matrix has 4 rows and 9 columns with next scheme:
(x x x - - - - - -)
(- - - x x x- - -)
(- - -- - - x x x)
(x- - x - - x - -)


In this cases youshould use row_num =4 instead of row_num = 9.

Victor

[Victor_Gladkikh · Accepted solution]

Or you can keep row_num = 9 and extend RowIndex: integer, dimension(10) ::rowindex=(/1,4,7,10,13,13,13,13,13,13/)

Victor

[mecej4]

What is the size of your matrix? Could you show the full matrix in standard mathematical notation? Then we could check whether you are describing it in compact row-major storage format.

See also the examples that came with MKL or the ones in the MKL documentation.

[Wu__Wentao]

Victor, Thank you so much!

This is the right answer I am expecting. I have tried and it worked successfully.

But setting the row number to be 4 seems not a perfect idea, though it returns right answer if transa = 'N'. The result is incorrect without any warning and error reports from compiler when transa = 'T'. However, the extended rowIndex works no matter transa = 'N' or 'T'.

Best regards,

David.

[Sergey_P_Intel2]

Hi David,

In MKL manual this functionality is described as working with square matrices:

The mkl_?csrgemv routine performs a matrix-vector operation defined as

y := A*x

or

y := A'*x,

where: x and y are vectors, A is an m-by-m sparse square matrix in the CSR format (3-array variation), A' is the transpose of A.

So for rectangular matrices it works improperly in some cases.

Regards,

Sergey

[Victor_Gladkikh]

Anyway the second solution have to work because we extend rectangular matrixto square matrix.

Victor

[Wu__Wentao]

Quoting Victor Gladkikh (Intel)

Anyway the second solution have to work because we extend rectangular matrixto square matrix.

Victor

Victor, Thank you very much again.

While, When you said we extend the rectanglar matrix to a square one, I get a new idea:

If we have a exact matrix,

(/ x x x x

x 1 2 x

x x x x

3 x 4 x

x x x x

x 5 6 x /)

then it can be presented in a sparse form with CSR (3 array variation)

value = (/1, 2, 3, 4, 5, 6/)

column = (/2, 3, 1, 3, 2, 3/)

rowIndex = (/1, 1, 3, 3, 5, 5, 7/).

That is to say, if there is a trival row(all zero elements row) in the matrix, a duplicated row index can be used to represent this row. I have tried some examples, and it does work!

But according to the MKL reference manual, "The 3-array variation of the CSR format has a restriction: all non-zero elements are stored continuously, that is the set of non-zero elements in the rowJgoes just after the set of non-zero elements in the rowJ-1."

What do you think of it, does my duplicated row index method always work? Or it just fits for the lucky simple cases.

Regrads,

David.

[mecej4]

Your proposal, it appears to me, goes against the principal reason for the existence of sparse storage schemes: to avoid calculations with zero elements.

Once you allow more than the minimal number of zeros (this number is non-zero for banded matrices, zero in packed storage schemes), the efficiency of the sparse matrix method deteriorates.