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

Converting from GSL to MKL LAPACK routines.... confusion

dsl_man
Beginner
1,241 Views
Hi everyone,

I am new to the Intel MKL. I converting some code from the GSL to the MKL LAPACK routines called from C++ on linux. I'm using MKL 10.0.1.014 and icpc 10.1.015.

I want to solve the square system Ax=b.
Currently I set up A like so:

for (int i=0;i for (int j=0;j gsl_matrix_set(A,i,j,some_value) etc....
}
}

So I want to call dgetrf to do the LU factorisation first something like:

dgetrf(&m,&n,&A[0][0],&m,&ipiv,&info);

A is set in the same place and gsl_matrix_set above.

So my first point of confusion is, the result of the LU factorisation is very slightly different (without changing the array ordering) than the one given by the GSL routine, and the value of info is 2 on return which I thought meant that the system cant be solved (?).

I was also under the impression I needed to transpose the matrix for the fortran array model, I tried switching the assigment of A in the above code to A with no success.

I suppose I am missing something here, if anyone could help it would be much appreciated!

Best regards,

Alastair

0 Kudos
1 Reply
Michael_C_Intel4
Employee
1,241 Views

Hi Alastair,

GSL keeps matrices in row-major ordering, but MKL assumes column-major. So, to get the same decomposition you need to explicitly transposeyour matrix (suppose, dimension 4x3 with lda=6):

a11 a12 a13 xxx xxx xxx

a21 a22 a23 xxx xxx xxx

a31 a32 a33 xxx xxx xxx

a41 a42 a43 xxx xxx xxx

which is located in memory as { a11, a12, a13, xxx, xxx, xxx, a21, a22, a23, xxx, xxx, xxx, a31, a32, a33, xxx, xxx, xxx, a41, a42, a43, xxx, xxx, xxx }

into the column-major matrix a', which should be located as { a11, a21, a31, a41, a12, a22, a32, a42, a13, a23, a33, a43, xxx, xxx, xxx,xxx, xxx, xxx, xxx, xxx, xxx,xxx, xxx, xxx }

then apply MKL dgetrf the same as you wrote - in this sample with m=4, n=3.

Be careful about input data. If some_data in your code is a constant, you'll have matrix columns linear dependent, and will have info=2inevitably, since the second column will be equal to the first.

Michael.

0 Kudos
Reply