Problem with mkl_dcsrmultcsr

[asd__asdqwe]

Hello,

I don't get why in the attached code, if line 36 I set transa to 't', then I get a segmentation fault. The matrix A is symmetric, so there shouldn't be any other changes to do to make this work ? In its Fortran version, dcsr_multiplication.f, there is no problem when I switch trans line 82.

Thanks in advance for any help.

[Gennady_F_Intel]

What mkl version you are using? I can find dcsr_multiplication.f only into $MKLROOT/examples/spblasf/source/ directory. this is version 11.0.

[asd__asdqwe]

Sorry Sir, it turns out this file is not part of MKL, I translated the original Fortran file you were mentioning in C++. I attached it in the first post, do you spot what I am doing wrong ? Thanks for your help, and sorry about that.

[Sergey_P_Intel2]

Hi! We can reproduce the situation for transposed case (transa='t'). It looks like parameter 0 is wrong for the second call to mkl_dcsrmultcsr. After replacing them on any valid pointer the program worked fine without crashes. Could you reproduce the same behavior?

[asd__asdqwe]

Hello, Thanks for your answer. I'm not sure about what you are saying. Here is the diff of something I tried according to your comment: 36c36 < char transa = 'n'; --- > char transa = 't'; 47c47,48 < mkl_dcsrmultcsr(&transa, &job, &sort, &n, &n, &n, --- > char transb = 't'; > mkl_dcsrmultcsr(&transb, &job, &sort, &n, &n, &n, Unfortunately, I still get a segfault, but it still works fine when 't' is replaced by 'n' for both transa and transb.

[Sergey_P_Intel2]

Ok, let me describe the changes in details. I meant that parameter with the value '0' looks wrong for the second call to mkl_dcsrmultcsr: mkl_dcsrmultcsr(&transa, &job, &sort, &n, &n, &n, a, ja, ia, a, ja, ia, c, jc, ic, 0, &ierr); If I replace '0' on any valid pointer, test works fine with transa='t'. E.g. it works if '0' is replaced on &n.

[asd__asdqwe]

Oh ok, it seems indeed to fix the segmentation fault (I have yet to check if the output is correct). Is this the expected behavior though ? The doc I have in front of me clearly states that 'nzmax [..] is used only if request=0' which is not the case in my test program. Another quick question, is it planed to support dcsrmultcsr with arbitrary 'matdescr' (for example symmetric input, unsymmetric output). I'm using this routine to compute A*D where A is symmetric (and so only its lower part is stored) and D is diagonal (but might as well be a general matrix). What I'm doing right now is csradd(csrmultcsr(A, D), csrmultcsr(A', D)) (I don't forget to remove the diagonal terms of csrmultcsr(A', D)), but not only it is making the MKL segfaults (not anymore with your fix though), I don't think this is very efficient.

[asd__asdqwe]

Hello, I think the results are incorrect when transa='t'. If you add the following lines after the two calls to dcsrmultdcsr[cpp]for(unsigned int i = 0; i < n; ++i) { for(unsigned int j = ic - 1; j < ic[i + 1] - 1; ++j) std::cout << i + 1 << " " << jc << " " << c << std::endl; }[/cpp] You will see that when transa='n', the results are correct and the corresponding CSR is upper triangular, while when transa='t', the CSR is not lower triangular, but since A is upper triangular, A' is lower triangular, hence A'*A' should be lower triangular. Do you agree with what I am saying ? Do you get incorrect results too ? Thanks in advance for your help.

[Sergey_P_Intel2]

Hi, In MKL documentation this routine is described in the following way: The mkl_?csrmultcsr routine performs a matrix-matrix operation defined as C := op(A)*B where: A, B, C are the sparse matrices in the CSR format (3-array variation); op(A) is one of op(A) = A, or op(A) =A', or op(A) = conjg(A') . So transa='t' is applied to first matrix A in your example, hence the result B = A' * A should be general matrix but not lower triangular.

[asd__asdqwe]

Hi, Ah ! I knew I forgot about something, sorry about that. Do you have an answer to my previous question ?

qweasd q. wrote:

Another quick question, is it planed to support dcsrmultcsr with arbitrary 'matdescr' (for example symmetric input, unsymmetric output). I'm using this routine to compute A*D where A is symmetric (and so only its lower part is stored) and D is diagonal (but might as well be a general matrix). What I'm doing right now is csradd(csrmultcsr(A, D), csrmultcsr(A', D)) (I don't forget to remove the diagonal terms of csrmultcsr(A', D)), but not only it is making the MKL segfaults (not anymore with your fix though), I don't think this is very efficient.

Thanks for all your help !

[Sergey_P_Intel2]

From my point of view segmentation fault in your example is not expected, it looks like a bug and requires further investigation. Supporting arbitrary matdescra in this functionality is under discussion because it requires implementation of new interfaces. Right now we have no plans to do that in nearest releases.

[Gennady_F_Intel]

hello, pls check the issue with the latest 11.0 update3 and let us know the result.

[asd__asdqwe]

Hello,

It seems Update 3 breaks ARPACK. Because I also use it in my application, I won't update until Update 4. Thanks anyway, I hope everything will be alright in the next update.