topic Re: Sparse CSR Matrix-Vector Multiplication in Intel® oneAPI Math Kernel Library

Sparse CSR Matrix-Vector Multiplication

dar_io — Tue, 22 Apr 2008 13:10:15 GMT

Hello!

I'm a new user of the MKL, and I have a "big" problem that I don't know how to solve...

I'm trying to multiply a sparse matrix in the CSR format with a vector, and to do this I'm using the mkl_dcsrmv method.

Unfortunately, the result is not correct. All the elements in the resulting vector are zero.
Switching the 4th element of the matdescra parameter (in c++ matdescra[3] = ...) from 'C' to 'F' results in a "clamped" matrix-vector multiplication, where the first row and column are not considered anymore (that's not surprising since doing this the one-based indexing is enabled, as far as I know...).

I was expecting that using 'C' instead of 'F' in the matdescra parameter enable the zero-based indexing but as I said all the elements in the resulting vector become zero.

I'm using C++ and XCode on a Mac Pro 8-Cores. The version of my MKL is 10.0.2.018.

I hope anyone can help me to solve this!
Thank you in advance, best!!

dario

PS: sorry my english:-)

Re: Sparse CSR Matrix-Vector Multiplication

Zhanghong_T_ — Tue, 22 Apr 2008 15:01:54 GMT

Dear Dario,

It seems to be a bug of MKL 10. Try to use MKL 9.1 and you can get a correct result.

Good luck!

Zhanghong Tang

Re: Sparse CSR Matrix-Vector Multiplication

dar_io — Tue, 22 Apr 2008 15:15:33 GMT

Hello,

Thank you very much for your answer!! I will try the MKL 9.1!
Best,

dario

Re: Sparse CSR Matrix-Vector Multiplication

Intel_C_Intel — Thu, 04 Sep 2008 12:37:20 GMT

Dear Dario and Zhanghong,

the mkl_dcsrmv function seems to work with MKL 10 Update 1 (or maybe 9.1), but not with Update 4.

I just updated my libraries and found my C program not working correct anymore. Having matdescra[5]="G F" and added 1 to every indexing entry, everything is fine and the result is equal to the correct result. Having matdescra[5]="G C", the result is a zero-valued vector.

This bug should be fixed asap, or not?

Fabian

Re: Sparse CSR Matrix-Vector Multiplication

Gennady_F_Intel — Mon, 08 Sep 2008 13:04:50 GMT

Hi Fabian,

This is unknown problem for us.Could you get us the test case? We will check the problem with the latest version.

--Gennady

Re: Sparse CSR Matrix-Vector Multiplication

Fabian_K — Thu, 11 Sep 2008 16:12:04 GMT

Hi Gennady,

this is the code for a test case (identity matrix * (0, 1, ... 99)) where I get wrong results:

  // initialize (dim == nonzeros)
  for (int i = 0; i < dim; i++) {
    x = (double)i;
    mkl_c = 1.;
    mkl_f = 0.;
    loop = 0.;
    values = 1.;  // identity matrix
    columns = i;
    pointerB = i;
    pointerE = i+1;
  }

  // Zero-based multiplication
  char transa = 'n';
  char matdescra_c[6] = "G  C ";
  double alpha = 1.;
  double beta = 0.;
  mkl_dcsrmv(&transa, &dim, &dim, α, matdescra_c, values, columns, pointerB, pointerE, x, β, mkl_c);

  // manual multiplication (loop)
  for (int i = 0; i    int fromInd = pointerB;
    int toInd = pointerE;
    for (int j=fromInd; j += x[columns] * values;
  }

  // One-based multiplication
  for (int i = 0; i < dim; i++) {
    columns += 1;
    pointerB += 1;
    pointerE += 1;
  }
  
  char matdescra_f[6] = "G  F ";
  mkl_dcsrmv(&transa, &dim, &dim, α, matdescra_f, values, columns, pointerB, pointerE, x, β, mkl_f);

  // output
  cout << endl << "mkl_C: "; for (int i = 0; i < dim; i++) { cout << mkl_c << " "; }
  cout << endl << "mkl_F: "; for (int i = 0; i < dim; i++) { cout << mkl_f << " "; }
  cout << endl << "loop : "; for (int i = 0; i < dim; i++) { cout << loop << " "; }

I tried it with version 10.0.011 (icpc 10.1 20080602), and the arrays mkl_C, mkl_F and loop contain the correct result (0, 1, ..99).

With version 10.0.4.023 (icpc 10.1 20080801), the array mkl_C is still zero (it is not changed, I tried other initial values). The other arrays contain the the correct result.

I also played with matdescra length and content, also with alpha and beta. This is a minimal example, of course. I use a lot of other MKL routines, and they all seem to work in both versions.

Thank you for your interest,
Fabian

Re: Sparse CSR Matrix-Vector Multiplication

Gennady_F_Intel — Tue, 16 Sep 2008 09:25:25 GMT

Hi Fabian,

This problem has been fixed and the fix is available in our MKL 10.0 update 4 which is released recently. Please download the same and verify the fix.

Thanks, Gennady

Re: Sparse CSR Matrix-Vector Multiplication

Gennady_F_Intel — Tue, 16 Sep 2008 09:28:50 GMT

Fabian,

Sorry for my previous message - I misprinted the MKL version: the problem was fixed in MKL 10.1 beta 2.

Please download the same and verify the fix.

Thanks, Gennady

Re: Sparse CSR Matrix-Vector Multiplication

Fabian_K — Fri, 26 Sep 2008 14:18:56 GMT

Quoting - gfedorov

Fabian,

Sorry for my previous message - I misprinted the MKL version: the problem was fixed in MKL 10.1 beta 2.

Please download the same and verify the fix.

Thanks, Gennady

Thanks.

In the recently published MKL 10.0.5.25 the problem is not fixed. Unfortunately, our admin don't like betas (like in your suggestion)... and me neither. So I use my own (not optimized) multiplication routine. That's okay for me.

Re: Sparse CSR Matrix-Vector Multiplication

Zhanghong_T_ — Tue, 04 May 2010 13:39:26 GMT

Dear all,

It seems that the bug of "mkl_dcsrmv" still exist on Linux version of MKL released together with Intel Fortran compiler 11.1.072 (I tested the intel64 version).

Thanks,
Zhanghong Tang

Re: Sparse CSR Matrix-Vector Multiplication

Gennady_F_Intel — Wed, 05 May 2010 07:45:31 GMT

HelloZhanghong,

1)it will be better if you give us the exact example for reproducing the problem.

2)how do you link the example?

3) CPU type you are working on?

4) the expected result

--Gennady

I am testing the functions of

Adrian_A_ — Thu, 23 Oct 2014 21:19:00 GMT

Hello, I am testing the functions of Intel MKL in a C test-program and I found that I just can't make the spareblas mkl_scsrmm CSR one-based indexing work. I link the example running first mklvar_ia32.bat, and then with VS command prompt and msvs compiler, using nmake libia32 and selecting the function I want to export (which is one I created for this test programs from the examples, because I couldn't find the problem in my original code so I wanted to isolate the function). My CPU is Intel Core i7-3630QM 2,40 Ghz, its a 64 bit processor running Microsoft Windows 8 OS (not 8.1) 64 bits. The PC is a Toshiba laptop (Satellite). Here is the code:

#include <stdio.h>
#include "mkl_types.h"
#include "mkl_spblas.h"

int main() {
    #define M 2        
    #define NNZ 4
        
        MKL_INT        m = M, nnz = NNZ;
        
        float        values[NNZ]     = {2.0,4.0,4.0,2.0};
        MKL_INT        columns[NNZ]  = {1,2,1,2};
        MKL_INT        rowIndex[M+1] = {1,3,5};
        
    #define N 2      
      
        MKL_INT        n = N;
        float         b    = {2.0, 1.0, 5.0, 2.0};
        float         c    = {0.0, 0.0, 0.0, 0.0};
        float         alpha = 1.0, beta = 0.0;
        char          transa, uplo, nonunit;
        char          matdescra[6];
		MKL_INT		  i, j;
        
    transa = 'N';
    
        matdescra[0] = 'S';
        matdescra[1] = 'L';
        matdescra[2] = 'N';
        matdescra[3] = 'F';
        
        mkl_scsrmm(&transa, &m, &n, &m, &alpha, matdescra, values, columns, rowIndex, &(rowIndex[1]), &(b[0][0]), &n,  &beta, &(c[0][0]), &n);

        printf("                             \n");
        printf("   OUTPUT DATA FOR MKL_SCSRMM\n");
        for (i = 0; i < m; i++) {
            for (j = 0; j < n; j++) {
                printf("%7.1f", c);
            };
            printf("\n");
        };
        return 0;
}

The results I get are this:

Zero-based indexing (the right solution):

24.0 10.0

18.0 8.0

One-based indexing:

8.0 10.0

18.0 24.0

Though it seems like it only changes the diagonal elements position, with 3x3 matrix the solution its completly different from the right one. I suspected that it might be something with the input format of the matrix b. I think that there's lack of clarity on the description of array b for the function mkl_scsrmm placed in the MKL reference manual. Thus I changed b format in this example and it worked (I placed elements position in this order: `{2.0, 5.0, 1.0, 2.0}). But I did the same for another 3x3 example I coded and it didn't work so I think it may be just a coincidence. I don't really know what to do with this problem, I have lost three days thinking that there was a problem in my original code (I am a begginer, industrial engineer learning technical software developtment on my own for a final degree proyect, sorry if I miss some important data/information, ask for it and I'll respond as soon as posible), working with mkl core examples in c and trying to figure out what happened. I would like to understand what's going on here.

References:

CSR format

https://software.intel.com/en-us/node/471374

Spare Blas mkl_scsrmm function

https://software.intel.com/sites/products/documentation/doclib/iss/2013/mkl/mklman/hh_goto.htm#GUID-78C55D9B-86FF-4A9F-B5D5-D2F61B9314FC.htm

You are calling a library

mecej4 — Thu, 23 Oct 2014 21:32:00 GMT

You are calling a library routine that was originally designed to be called mainly from Fortran, which uses column-major 2-D arrays, i.e., the columns are contiguous in memory. A C-native 2-D array is kept in row-major format. The result is that when you pass a 2-D array A from C, the library routine interprets the argument as C^T.

But then why the intel c

Adrian_A_ — Fri, 24 Oct 2014 09:53:00 GMT

But then why in the intel c example, which uses zero-based indexing, they don't transpose the matrix b and the result is right?

I don't know which example

mecej4 — Fri, 24 Oct 2014 12:29:26 GMT

I don't know which example you refer to, nor can I compile your example program of #12 and run it to see what the problem is (there are undeclared variables).

As I said before, note that the result array C, which you try to print using the name rhs, is delivered by MKL as a column-major matrix, so your program will print out C^T.

Sorry for the code of #12 I

Adrian_A_ — Fri, 24 Oct 2014 15:01:00 GMT

Sorry for the code of #12 I forgot to copy two instructions and the c name in the printf, now you should be able to compile it. The intel example I am using is located in: "...mkl\examples\examples_core_c\spblasc\source\cspblas_scsr.c". I attached the other example I am using. There are two examples in that code, doing AxB, with the same matrix A in zero-based indexing and in one-based indexing. With zero you don't need to store the matrix B in column major order and everything seems to work. With one based indexing I store the matrix B in column order as you said (ie transposed) and I print both C and C'. None of them are the solution. I don't know what happens. Maby the one based indexing is made only for being called from Fortran...but this is a C interface...

Here is what you said, the paragraph before the title "Differences Between Intel MKL and NIST* Interfaces"

https://software.intel.com/sites/products/documentation/doclib/iss/2013/mkl/mklman/hh_goto.htm#GUID-34C8DB79-0139-46E0-8B53-99F3BEE7B2D4.htm

I read this before coding my program, but as I saw that in the examples coded in c, the native c storage format was used, I thought that the interfaces for C arranged the arrays so that they get to the Fortran code in the proper format. Now I don't know what happens.

The program in #16 works

mecej4 — Fri, 24 Oct 2014 17:50:13 GMT

The program in #16 works correctly, and I think that the problem you have is in interpreting the results. You print C and C^T separately, which is redundant, since visually transposing a 3 X 3 matrix is something that most of us can accomplish.

The first pair of results that your program prints out are (A.B) and (A.B)^T. The second pair of results are (A^T.B^T) and (A^T.B^T) ^T. The last item is equal to B.A but, because B is not symmetric, it is not equal to A.B.

At least as far as your example is concerned, the issue of zero-based or one-based indices is a red herring. That issue affects only sparse matrix representations, and the only sparse matrix in your problem, A, is symmetric. Secondly, you chose square matrices for your example, so that the stride from a row to the next row equals the stride from a column to the next column, so the issue of matdescra[3] being 'C' or 'F' changes nothing.

You may find it beneficial to choose another example where the matrices are not square, so that the aspect of a matrix being symmetric goes away, and compare the results of hand calculations with what MKL gives you.

Well, finally I think I

Adrian_A_ — Sat, 25 Oct 2014 15:37:00 GMT

Well, finally I think I understand what happens. The C interface needs you to store arrays in row-major order when you use zero-based indexing and column-major order when you use one-based indexing, and it returns the arrays in that same order. My mistake was thinking that C interface would accept for both ways of indexing the row-major ordering, and translate it to the Fortran column-major ordering. The format of the sparse matrix just changes in the indexing, from zero to one. So when you call the routine for one and zero-based indexing, changing only matrix A indexing and entering the same matrix B_mxn{a11...a1n, a21....a2n,....,am1....amn} the routine will output the right result for cero-based indexing. For one-based indexing it will misunderstand the matrix B as if it was another matrix with coefficients in different positions, like this:

B = 1   2   3

    4   5   6

arrayB = {1 2 3 4 5 6}

column-order matrix B interpretation, which is used in one-order indexing:

B = 1   3   5

    2   4   6

And the result will also be in column-major order. The last paragraph of #17 really helped me.

Thank you very much, Mr. mecej4 for helping me. I really appreciate you took the time to explain explain this to me.

Regards.

Hi,

Lessig__Christian — Thu, 16 Jan 2020 08:26:11 GMT

Hi,

I have a problem involving sparse matrix-vector multiplication in a loop. In each loop iteration the number of rows that is used in the sparse matrix increases. Is there an efficient way to implement this with mkl_sparse_?_mv()?

Thanks!