Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.
Announcements
The Intel sign-in experience is changing in February to support enhanced security controls. If you sign in, click here for more information.

MKL function dgetrf wrong results

dimdol10
New Contributor I
597 Views

Hello!

I am using Intel oneAPI MKL version 2022.2, and visual studio 2022.

I am currently trying a LU decomposition of a matrix using the function DGETRF provided by MKL. 

 

As I found a wrong result in a simple example, I would like ask you if something is incorrect. The example code is attached in the below. As far as I understood, in this example, the array ipiv should be 

ipiv(1) = 2,

ipiv(2) = 1,

such that 

A_PLU = P*L*U. 

However, the function gives

ipiv(1) = 2,

ipiv(2) = 2.

 

---------------------------------------------------------

PROGRAM TEST_MAT_LU_FAC_INV

DOUBLE PRECISION A(2,2)
DOUBLE PRECISION A_L(2,2), A_U(2,2)
DOUBLE PRECISION A_LU(2,2)
DOUBLE PRECISION A_PLU(2,2), A_PLU_TEST(2,2)
integer info
integer ipiv(2)
DOUBLE PRECISION P(2,2) ! permutation matrix

A = 0.D0
A(1,1) = 1.D0
A(1,2) = 2.D0
A(2,1) = 3.D0
A(2,2) = 4.D0

CALL DGETRF(2,2,A,2,ipiv,info)
A_L = 0.d0
A_L(1,1) = 1.D0 !A(1,1)
A_L(2,1) = A(2,1)
A_L(2,2) = 1.D0 !A(2,2)

A_U = 0.D0
A_U(1,1) = A(1,1)
A_U(1,2) = A(1,2)
A_U(2,2) = A(2,2)

A_LU = 0.D0
DO I = 1,2
DO J = 1,2
DO K = 1,2
A_LU(I,J) = A_LU(I,J) + A_L(I,K)*A_U(K,J)
END DO
END DO
END DO

P = 0.D0
DO I = 1,2
P(I,IPIV(I)) = 1.D0
END DO

A_PLU_TEST = 0.D0
DO I = 1,2
DO J = 1,2
DO K = 1,2
A_PLU_TEST(I,J) = A_PLU_TEST(I,J) + A_L(I,K)*A_U(K,J)
END DO
END DO
END DO

A_PLU = 0.D0
DO I = 1,2
DO J = 1,2
DO K = 1,2
A_PLU(I,J) = A_PLU(I,J) + P(I,K)*A_PLU_TEST(K,J)
END DO
END DO
END DO


END

0 Kudos
5 Replies
ShanmukhS_Intel
Moderator
564 Views

Hi Myung,


Thanks for posting on Intel Communities.


Could you please cross-verify the shared source code? We tried compiling it and we couldn't see any output.


Best Regards,

Shanmukh.SS


dimdol10
New Contributor I
555 Views

Dear Shanmukh,

 

Thanks for the reply. 

Yes, the code print out nothing. But, what I meant is that the array A_PLU is supposed to be the same with the array A. However, if you run the code, you can find that they are different. 

 

A_PLU(1,1) = 1.D0

A_PLU(2,1) = 1.D0

A_PLU(1,2) = 2.D0

A_PLU(2,2) = 2.D0

 

I suspect the output array ipiv from DGETRF is wrong. Could you please take a look at this point?

 

Thank you very much.

 

ShanmukhS_Intel
Moderator
518 Views

Hi Myung,


Please find the below link for more details and usage regarding the DGETRF functionality.


https://www.intel.com/content/www/us/en/develop/documentation/onemkl-developer-reference-fortran/top...


The function describes the input parameters to be used and the proper usage of the routine.


Best Regards,

Shanmukh.SS


ShanmukhS_Intel
Moderator
477 Views

Hi,


A gentle reminder:

We have not heard back from you for a while. Could you please let us know if any updates regarding your issue?


Best Regards,

Shanmukh.SS


ShanmukhS_Intel
Moderator
437 Views

Hi Myung,


We assume that your issue is resolved. If you need any other information, please post a new question as this thread will no longer be monitored by Intel.


Best Regards,

Shanmukh.SS


Reply