Errors using mkl_ddiasv

[chubbz]

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
y.resize(m);

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

[chubbz]

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:

0.322399527186761

0.528959810874704

-0.402186761229315

0.870567375886525

-0.153664302600473

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

For my test matrix, I'm using the matrix C fromhttp://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htm

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

Code I used for the above tests:

double values[] = { 666,666,666,4,8,

666,-2,0,2,0,

1,5,4,7,-5,

-1,0,6,0,666,

-3,0,4,666,666};

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]);

[Ying_H_Intel]

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 fromhttp://software.intel.com/sites/products/documentation/hpc/mkl/mklman/index.htmis 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,

where:

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.