Community
cancel
Showing results for 
Search instead for 
Did you mean: 
chubbz
Beginner
77 Views

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

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

0 Kudos
2 Replies
chubbz
Beginner
77 Views

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

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
Employee
77 Views

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.

Reply