Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

Errors using mkl_ddiasv

Ive been running into issues while trying to use the Diagonal matrix solver routine, mkl_ddiasv. My results vary wildly from what I should be getting, when compared to the results from matlab. The only possible issue is that my matrix is stored in row-major order, but transposing the matrix should take care of that. I've also tried converting my matrix to column-major order and it still didnt work. My calling code is below:

char transa = 'C';
int m = 100; //the size of the main diagonal - 100
double alpha = 1.0;
char matdescr[] = { 'D', 'U', 'N', 'C' };
//val is created above
int lval = m; //100
int distance[] = {-10, -1, 0, 1, 10};
int ndiag = 5; //5

mkl_ddiasv(&transa, &m, α, matdescr, &val[0], &lval, distance, &ndiag, &[0], &y[0]);

0 Kudos
2 Replies
Alright, I scaled down to a smaller matrix to test out my implementation of the library,I've tried converting my matrix to CSR format, in order to verify I was using the diagonal format correctly. The matrix appeared to be parsing properly. After confirming this, I then used the dcsrsv routine to see whether it was just an error with ddiasv, but I'm still getting the wrong numbers.

According to MATLAB, the resulting matrix should be:

However the various MKL routines are outuptting:
[0] 1.0000000000000000 double
[1] 0.40000000000000002 double
[2] 0.75000000000000000 double
[3] 0.57142857142857140 double
[4] -1.0000000000000000 double

and my RHS = [1,2,3,4,5]

Code I used for the above tests:
double values[] = { 666,666,666,4,8,
char transa = 'N';
int m = 5;
double alpha = 1.0;
char matdescr[] = { 'D', 'U', 'N', 'C' };
int lval = m; //100
int distance[] = {-3,-1,0,1,2};
int ndiag = 5; //5
double rhs[] = { 1,2,3,4,5 };
std::vector y (5);
int job[] = {1,0,0,0,0,0};
std::vector ascr (13);
std::vector ja (13);
std::vector ia (6);
//mkl_ddiasv(&transa, &m, α, matdescr, values, &lval, distance, &ndiag, rhs, &y[0]);
mkl_dcsrdia(job, &m, &ascr[0], &ja[0], &ia[0], &values[0], &m ,distance, &m, NULL, NULL,NULL, NULL);
mkl_dcsrsv(&transa, &m, α, matdescr, &ascr[0], &ja[0], &ia[0], &ia[1], rhs, &y[0]);
0 Kudos

Hi Chubbz,

As I see, you handled the input parameters and matrixstorage formatproperly. The key problem is that the mkl_ddiasv and mkl_dcsrsvis for triangular sparse matrix.while thematrix C from not a triangular matrx.

See the mkl manual:

The mkl_?csrsv routine solves a system of linear equations with matrix-vector operations for a sparse matrix in the CSR format: y := alpha*inv(A)*x

or y := alpha*inv(A')*x,


alpha is scalar, x and y are vectors, A is a sparse upper or lower triangular matrix with unit or non-unit main diagonal, A' is the transpose of A.

So the output arealways not correct. You may try lapack's function orSparse Solver Routines like pardiso/dss to solve the sparse linear system of equations.

For the function mkl_ddiasv and mkl_dcsrsvthemselves, you may testthe mkl example codes likecspblas_dcsr.c,which be default are in in
Best Regards,
Ying H.

0 Kudos